diff --git a/.changeset/strang-poa-pv-forecast.md b/.changeset/strang-poa-pv-forecast.md new file mode 100644 index 00000000..687c8a7e --- /dev/null +++ b/.changeset/strang-poa-pv-forecast.md @@ -0,0 +1,15 @@ +--- +"ftw": minor +--- + +PV forecasts are now orientation-aware. When per-plane geometry is configured +(the Weather tab's PV arrays: tilt/azimuth/kWp), a radiation-bearing forecast +provider's global horizontal irradiance is projected onto each panel plane via +the physics `sunpos` model and summed, instead of the previous flat +`rated × (W/m² / 1000)` estimate that ignored panel orientation. Sites with no +arrays configured keep the existing behaviour, and providers that already return +site-calibrated watts (Forecast.Solar) are left untouched. + +Providers that publish only global horizontal irradiance get an Erbs correlation +to split it into direct and diffuse components before projection, so a +south-facing 35° roof and a flat one no longer receive the same forecast. diff --git a/go/internal/forecast/forecast.go b/go/internal/forecast/forecast.go index 08a55efb..6222871b 100644 --- a/go/internal/forecast/forecast.go +++ b/go/internal/forecast/forecast.go @@ -29,6 +29,7 @@ import ( "github.com/srcfl/ftw/go/internal/config" "github.com/srcfl/ftw/go/internal/state" + "github.com/srcfl/ftw/go/internal/sunpos" ) // Provider is implemented by each weather source. @@ -243,10 +244,16 @@ func EstimatePVW(lat, lon float64, t time.Time, cloudPct *float64, ratedW float6 // Service wraps a provider + store + scheduler for forecasts. type Service struct { - Provider Provider - Store *state.Store - Lat, Lon float64 - RatedPVW float64 // total rated PV across all arrays (used for estimate) + Provider Provider + Store *state.Store + Lat, Lon float64 + RatedPVW float64 // total rated PV across all arrays (used for estimate) + + // Arrays holds per-plane geometry (tilt/azimuth/kWp) mirrored from the + // weather config. When set, a radiation-bearing provider's horizontal + // GHI is projected onto each plane via sunpos and summed, instead of the + // orientation-blind flat rated×(W/m²/1000) estimate. Empty → flat estimate. + Arrays []Array stop chan struct{} done chan struct{} @@ -286,10 +293,21 @@ func FromConfig(cfg *config.Weather, ratedPVW float64, st *state.Store, userAgen default: return nil } + // Mirror per-plane geometry for the POA path. Shared across all + // providers: forecast_solar already applies geometry server-side (so + // this stays unused there), but open_meteo / STRÅNG-style GHI providers + // use it to project irradiance onto each plane in fetchAndStore. + var arrays []Array + for _, a := range cfg.PVArrays { + if a.KWp > 0 { + arrays = append(arrays, Array{TiltDeg: a.TiltDeg, AzimuthDeg: a.AzimuthDeg, KWp: a.KWp}) + } + } return &Service{ Provider: p, Store: st, Lat: cfg.Latitude, Lon: cfg.Longitude, RatedPVW: ratedPVW, + Arrays: arrays, stop: make(chan struct{}), done: make(chan struct{}), } @@ -341,7 +359,11 @@ func (s *Service) fetchAndStore(ctx context.Context) { switch { case r.PVWEstimated != nil: pvW = *r.PVWEstimated + case r.SolarWm2 != nil && len(s.Arrays) > 0: + // Orientation-aware: project GHI onto each configured plane and sum. + pvW = poaPVWattsFromGHI(s.Lat, s.Lon, r.HourStart, *r.SolarWm2, s.Arrays) case r.SolarWm2 != nil && s.RatedPVW > 0: + // No geometry configured — flat, orientation-blind estimate. pvW = s.RatedPVW * (*r.SolarWm2) / 1000.0 default: pvW = EstimatePVW(s.Lat, s.Lon, r.HourStart, r.CloudCoverPct, s.RatedPVW) @@ -365,6 +387,25 @@ func (s *Service) fetchAndStore(ctx context.Context) { slog.Info("forecast fetched", "count", len(points), "provider", s.Provider.Name()) } +// poaPVWattsFromGHI converts a global-horizontal irradiance (W/m², positive) +// into expected DC PV output (W, positive) by projecting it onto each +// configured array's plane via sunpos and scaling by nameplate. This is the +// orientation-aware replacement for the flat rated×(W/m²/1000) estimate; it +// is used whenever the provider supplies GHI and the site has per-plane +// geometry. Returns 0 when the sun is down or no arrays produce output. +func poaPVWattsFromGHI(lat, lon float64, t time.Time, ghiWm2 float64, arrays []Array) float64 { + var total float64 + for _, a := range arrays { + if a.KWp <= 0 { + continue + } + poa := sunpos.POAFromGHI(t, lat, lon, ghiWm2, a.TiltDeg, a.AzimuthDeg) + // kWp×1000 = nameplate W at STC (1000 W/m²); scale by POA/1000. + total += a.KWp * 1000.0 * (poa / 1000.0) + } + return total +} + // Load returns forecasts in [sinceMs, untilMs]. func (s *Service) Load(sinceMs, untilMs int64) ([]state.ForecastPoint, error) { return s.Store.LoadForecasts(sinceMs, untilMs) diff --git a/go/internal/forecast/forecast_test.go b/go/internal/forecast/forecast_test.go index 73450fa0..e0cb5526 100644 --- a/go/internal/forecast/forecast_test.go +++ b/go/internal/forecast/forecast_test.go @@ -247,3 +247,97 @@ func TestFromConfigBuildsMetNo(t *testing.T) { if s.Lat != 59 { t.Errorf("lat: %f", s.Lat) } if s.RatedPVW != 10000 { t.Errorf("rated: %f", s.RatedPVW) } } + +func TestFromConfigPopulatesArrays(t *testing.T) { + st, _ := state.Open(filepath.Join(t.TempDir(), "t.db")) + defer st.Close() + cfg := &config.Weather{ + Provider: "open_meteo", Latitude: 59, Longitude: 18, + PVArrays: []config.PVArray{ + {Name: "south", KWp: 6, TiltDeg: 35, AzimuthDeg: 180}, + {Name: "east", KWp: 4, TiltDeg: 30, AzimuthDeg: 90}, + {Name: "empty", KWp: 0, TiltDeg: 10, AzimuthDeg: 200}, // skipped (kWp 0) + }, + } + s := FromConfig(cfg, 10000, st, "ua") + if s == nil { t.Fatal("expected service") } + if len(s.Arrays) != 2 { + t.Fatalf("expected 2 arrays (kWp>0 only), got %d", len(s.Arrays)) + } + if s.Arrays[0].KWp != 6 || s.Arrays[1].AzimuthDeg != 90 { + t.Errorf("array geometry mismatch: %+v", s.Arrays) + } +} + +// ---- POA-per-array (orientation-aware) estimate ---- + +func TestPOAPVWattsSumsArrays(t *testing.T) { + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + one := poaPVWattsFromGHI(59.3293, 18.0686, tt, 700, []Array{{TiltDeg: 35, AzimuthDeg: 180, KWp: 5}}) + two := poaPVWattsFromGHI(59.3293, 18.0686, tt, 700, []Array{ + {TiltDeg: 35, AzimuthDeg: 180, KWp: 5}, + {TiltDeg: 35, AzimuthDeg: 180, KWp: 5}, + }) + if one <= 0 { + t.Fatalf("expected positive POA watts, got %.1f", one) + } + if math.Abs(two-2*one) > 1e-6 { + t.Errorf("two identical arrays should double output: one=%.2f two=%.2f", one, two) + } +} + +func TestPOAPVWattsZeroAtNight(t *testing.T) { + tt := time.Date(2026, 12, 21, 23, 0, 0, 0, time.UTC) + w := poaPVWattsFromGHI(59.3293, 18.0686, tt, 500, []Array{{TiltDeg: 35, AzimuthDeg: 180, KWp: 10}}) + if w != 0 { + t.Errorf("night POA watts should be 0, got %.2f", w) + } +} + +// End-to-end: with per-plane geometry, a GHI-bearing provider's stored PV +// estimate comes from the POA path and differs from the flat rated×GHI/1000. +func TestServicePOAPathDiffersFromFlat(t *testing.T) { + handler := http.HandlerFunc(func(w http.ResponseWriter, r *http.Request) { + resp := map[string]any{ + "hourly": map[string]any{ + "time": []string{"2026-06-21T11:00"}, + "shortwave_radiation": []float64{700}, + "cloud_cover": []float64{5}, + "temperature_2m": []float64{20}, + }, + } + _ = json.NewEncoder(w).Encode(resp) + }) + srv := httptest.NewServer(handler) + defer srv.Close() + + st, _ := state.Open(filepath.Join(t.TempDir(), "state.db")) + defer st.Close() + + p := NewOpenMeteo() + p.BaseURL = srv.URL + s := &Service{ + Provider: p, Store: st, Lat: 59.3293, Lon: 18.0686, RatedPVW: 10000, + Arrays: []Array{{TiltDeg: 35, AzimuthDeg: 180, KWp: 10}}, + } + s.fetchAndStore(context.Background()) + + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + rows, err := st.LoadForecasts(tt.UnixMilli(), tt.Add(time.Hour).UnixMilli()) + if err != nil { + t.Fatal(err) + } + if len(rows) != 1 || rows[0].PVWEstimated == nil { + t.Fatalf("expected 1 forecast with PV estimate, got %+v", rows) + } + got := *rows[0].PVWEstimated + flat := 10000 * 700.0 / 1000.0 // orientation-blind estimate = 7000 W + want := poaPVWattsFromGHI(59.3293, 18.0686, tt, 700, s.Arrays) + if math.Abs(got-want) > 1.0 { + t.Errorf("service should use POA path: got %.1f want %.1f", got, want) + } + if math.Abs(got-flat) < 1.0 { + t.Errorf("POA estimate should differ from flat %.0f, got %.1f", flat, got) + } + t.Logf("POA-per-array estimate %.0fW vs flat %.0fW", got, flat) +} diff --git a/go/internal/sunpos/sunpos.go b/go/internal/sunpos/sunpos.go index 8ce95cc3..787e7c27 100644 --- a/go/internal/sunpos/sunpos.go +++ b/go/internal/sunpos/sunpos.go @@ -125,35 +125,117 @@ func ClearSkyW(t time.Time, lat, lon float64) float64 { } // POA estimates plane-of-array irradiance for one tilted panel using the -// isotropic-sky model. Splits clear-sky horizontal irradiance into beam -// (DNI) and diffuse (DHI) components via a simple Erbs correlation, then -// projects each onto the panel. +// isotropic-sky model, driven by the package's own clear-sky prior. It splits +// that clear-sky horizontal irradiance into beam (DNI) and diffuse (DHI) +// components with a fixed 20% diffuse fraction, then projects each onto the +// panel. Used as the prior signal for the PV twin when no measured irradiance +// is available. // // Returns W/m² on the panel surface; clamped to ≥ 0. +// +// When a data source supplies measured irradiance, prefer the two variants +// below: POAFromComponents (GHI + DHI both known, e.g. SMHI STRÅNG params +// 117 + 122) or POAFromGHI (only GHI known, e.g. Open-Meteo shortwave). func POA(t time.Time, lat, lon, panelTiltDeg, panelAzDeg float64) float64 { sun := At(t, lat, lon) - if sun.ZenithDeg >= 90 { - return 0 - } ghi := ClearSkyW(t, lat, lon) - if ghi <= 0 { + // No measured diffuse component from the clear-sky prior, so keep the + // historical fixed 20% diffuse fraction for this variant. + return POAFromComponents(sun, ghi, 0.2*ghi, panelTiltDeg, panelAzDeg) +} + +// POAFromComponents projects measured global (GHI) and diffuse (DHI) +// horizontal irradiance onto a tilted panel using the isotropic-sky model, +// reusing AOI for the beam projection. All irradiances in W/m²; returns the +// plane-of-array irradiance in W/m², clamped ≥ 0. +// +// Use this when a source gives both GHI and DHI directly (e.g. SMHI STRÅNG +// parameters 117 + 122). When only GHI is available use POAFromGHI, which +// estimates the diffuse split via the Erbs correlation first. +func POAFromComponents(sun Position, ghi, dhi, panelTiltDeg, panelAzDeg float64) float64 { + if sun.ZenithDeg >= 90 || ghi <= 0 { return 0 } - // Erbs et al. (1982) clearness-based diffuse fraction. We don't have - // real GHI measurements, so kt comes from our own clear-sky model → - // always ~0.7-0.75 → diffuse ratio ~0.2. Conservative. - kd := 0.2 - dhi := ghi * kd - dni := (ghi - dhi) / math.Cos(sun.ZenithDeg*math.Pi/180) - + if dhi < 0 { + dhi = 0 + } + if dhi > ghi { + dhi = ghi + } + tiltR := panelTiltDeg * math.Pi / 180 + diffusePOA := dhi * (1 + math.Cos(tiltR)) / 2 // isotropic sky dome aoi := AOI(sun, panelTiltDeg, panelAzDeg) if aoi > 90 { - // Sun behind panel — only diffuse counts. - return dhi * (1 + math.Cos(panelTiltDeg*math.Pi/180)) / 2 + // Sun behind the panel — only diffuse reaches the surface. + return diffusePOA } + cosZ := math.Cos(sun.ZenithDeg * math.Pi / 180) + if cosZ < 0.01 { + // Sun on the horizon: beam projection is ill-conditioned + // (divide-by-~0) and diffuse dominates anyway. + return diffusePOA + } + dni := (ghi - dhi) / cosZ beamPOA := dni * math.Cos(aoi*math.Pi/180) - diffusePOA := dhi * (1 + math.Cos(panelTiltDeg*math.Pi/180)) / 2 out := beamPOA + diffusePOA - if out < 0 { out = 0 } + if out < 0 { + out = 0 + } return out } + +// POAFromGHI projects a measured/forecast global horizontal irradiance (GHI, +// W/m²) onto a tilted panel when no diffuse component is available. It +// estimates the diffuse fraction from the hourly clearness index via the Erbs +// et al. (1982) correlation, then delegates to POAFromComponents. +// +// Use this for radiation providers that expose shortwave/GHI but not diffuse +// (e.g. Open-Meteo shortwave_radiation, or SMHI STRÅNG global-only windows). +func POAFromGHI(t time.Time, lat, lon, ghi, panelTiltDeg, panelAzDeg float64) float64 { + sun := At(t, lat, lon) + if sun.ZenithDeg >= 90 || ghi <= 0 { + return 0 + } + cosZ := math.Cos(sun.ZenithDeg * math.Pi / 180) + i0h := extraterrestrialHorizontalW(t, cosZ) + kt := 0.0 + if i0h > 0 { + kt = ghi / i0h + } + dhi := ghi * ErbsDiffuseFraction(kt) + return POAFromComponents(sun, ghi, dhi, panelTiltDeg, panelAzDeg) +} + +// ErbsDiffuseFraction returns the diffuse fraction (DHI/GHI) for an hourly +// clearness index kt, per Erbs, Klein & Duffie (1982). Result is in +// [0.165, 1]: overcast skies (low kt) are almost entirely diffuse, clear +// skies (high kt) settle near 16.5% diffuse. +func ErbsDiffuseFraction(kt float64) float64 { + switch { + case kt <= 0: + return 1 + case kt <= 0.22: + return 1 - 0.09*kt + case kt <= 0.80: + return 0.9511 - 0.1604*kt + 4.388*kt*kt - 16.638*kt*kt*kt + 12.336*kt*kt*kt*kt + default: + return 0.165 + } +} + +// extraterrestrialHorizontalW returns top-of-atmosphere irradiance on a +// horizontal surface (W/m²) at time t for a solar cosine-zenith cosZ. Used as +// the denominator of the clearness index. Returns 0 when the sun is at/below +// the horizon. +func extraterrestrialHorizontalW(t time.Time, cosZ float64) float64 { + if cosZ <= 0 { + return 0 + } + doy := float64(t.UTC().YearDay()) + gamma := 2 * math.Pi * (doy - 1) / 365 + e0 := 1.000110 + + 0.034221*math.Cos(gamma) + 0.001280*math.Sin(gamma) + + 0.000719*math.Cos(2*gamma) + 0.000077*math.Sin(2*gamma) + const i0 = 1361.0 // solar constant W/m² + return i0 * e0 * cosZ +} diff --git a/go/internal/sunpos/sunpos_test.go b/go/internal/sunpos/sunpos_test.go index efe085fa..fb48f8d5 100644 --- a/go/internal/sunpos/sunpos_test.go +++ b/go/internal/sunpos/sunpos_test.go @@ -66,3 +66,90 @@ func TestPOAFlatEqualsGHI(t *testing.T) { t.Errorf("flat POA should equal GHI: ghi=%.1f poa=%.1f", ghi, poa) } } + +// A flat panel receives all of GHI regardless of the diffuse split, because +// beam-on-horizontal + diffuse-on-horizontal reconstructs GHI exactly. +func TestPOAFromComponentsFlatEqualsGHI(t *testing.T) { + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + sun := At(tt, 59.33, 18.07) + const ghi = 700.0 + for _, dhi := range []float64{0, 140, 350, 700} { + poa := POAFromComponents(sun, ghi, dhi, 0, 180) + if math.Abs(poa-ghi) > 0.5 { + t.Errorf("flat POA should equal GHI for dhi=%.0f: got %.2f want %.0f", dhi, poa, ghi) + } + } +} + +// A south-tilted panel at Stockholm winter noon (low sun) collects more than a +// flat one — the whole point of projecting onto the plane. +func TestPOAFromComponentsSouthTiltBeatsFlatInWinter(t *testing.T) { + tt := time.Date(2026, 12, 21, 11, 0, 0, 0, time.UTC) // ~noon local, low sun + sun := At(tt, 59.33, 18.07) + if sun.ZenithDeg >= 90 { + t.Skip("sun below horizon") + } + const ghi, dhi = 200.0, 60.0 + flat := POAFromComponents(sun, ghi, dhi, 0, 180) + tilt := POAFromComponents(sun, ghi, dhi, 45, 180) // 45° south + if tilt <= flat { + t.Errorf("south-tilted winter POA (%.1f) should beat flat (%.1f)", tilt, flat) + } +} + +// Sun behind the panel → only the diffuse dome contributes (no negative beam). +func TestPOAFromComponentsSunBehindPanelIsDiffuseOnly(t *testing.T) { + tt := time.Date(2026, 6, 21, 5, 0, 0, 0, time.UTC) // morning, sun in the east + sun := At(tt, 59.33, 18.07) + if sun.ZenithDeg >= 90 { + t.Skip("sun below horizon") + } + const ghi, dhi = 300.0, 90.0 + // A steep west-facing wall can't see the eastern morning sun's beam. + poa := POAFromComponents(sun, ghi, dhi, 90, 270) + diffuseOnly := dhi * (1 + math.Cos(90*math.Pi/180)) / 2 + if math.Abs(poa-diffuseOnly) > 0.5 { + t.Errorf("sun-behind panel should be diffuse-only %.2f, got %.2f", diffuseOnly, poa) + } +} + +// Erbs diffuse fraction: overcast → ~all diffuse, clear → floor at 0.165, +// monotonically non-increasing across the mid range. +func TestErbsDiffuseFraction(t *testing.T) { + if f := ErbsDiffuseFraction(0); f != 1 { + t.Errorf("kt=0 → fraction 1, got %.3f", f) + } + if f := ErbsDiffuseFraction(1.0); math.Abs(f-0.165) > 1e-9 { + t.Errorf("kt>0.8 → 0.165, got %.3f", f) + } + clear := ErbsDiffuseFraction(0.75) + murky := ErbsDiffuseFraction(0.30) + if !(clear < murky) { + t.Errorf("clearer sky should have less diffuse: clear=%.3f murky=%.3f", clear, murky) + } + if clear < 0.165 || clear > 1 { + t.Errorf("fraction out of [0.165,1]: %.3f", clear) + } +} + +// POAFromGHI (Erbs split) on a flat panel still reconstructs ~GHI. +func TestPOAFromGHIFlatApproxGHI(t *testing.T) { + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + const ghi = 600.0 + poa := POAFromGHI(tt, 59.33, 18.07, ghi, 0, 180) + if math.Abs(poa-ghi) > 0.5 { + t.Errorf("flat POAFromGHI should equal GHI: got %.2f want %.0f", poa, ghi) + } +} + +// Night → zero from both measured-irradiance variants regardless of input. +func TestPOAVariantsZeroAtNight(t *testing.T) { + tt := time.Date(2026, 12, 21, 23, 0, 0, 0, time.UTC) + sun := At(tt, 59.33, 18.07) + if p := POAFromComponents(sun, 500, 100, 35, 180); p != 0 { + t.Errorf("night POAFromComponents should be 0, got %.2f", p) + } + if p := POAFromGHI(tt, 59.33, 18.07, 500, 35, 180); p != 0 { + t.Errorf("night POAFromGHI should be 0, got %.2f", p) + } +}