feat(forecast): project irradiance onto each PV array's plane

.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.

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)

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)
+}

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
+}