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/config/config.go b/go/internal/config/config.go index 074d84e8..21f5ed0d 100644 --- a/go/internal/config/config.go +++ b/go/internal/config/config.go @@ -978,9 +978,10 @@ type Weather struct { // predictions when each plane is described separately than when // everything is averaged into a single tilt/azimuth. // - // When set, PVArrays overrides the legacy single-array fields. - // Providers that can't use site geometry (met_no, open_meteo) - // ignore this entirely and just use PVRatedW. + // When set, PVArrays overrides the legacy single-array fields for + // geometry-aware providers. GHI providers such as open_meteo project + // radiation onto these planes; incomplete entries are ignored and the + // provider falls back to its flat estimate. PVArrays []PVArray `yaml:"pv_arrays,omitempty" json:"pv_arrays,omitempty"` // HeatingWPerDegC adds load proportional to max(18°C − outdoor_temp, 0). @@ -995,10 +996,27 @@ type Weather struct { // + east roof + garage) with different tilt/azimuth. The sum of all // KWp values should match the total PV nameplate at the site. type PVArray struct { - Name string `yaml:"name,omitempty" json:"name,omitempty"` - KWp float64 `yaml:"kwp" json:"kwp"` - TiltDeg float64 `yaml:"tilt_deg" json:"tilt_deg"` - AzimuthDeg float64 `yaml:"azimuth_deg" json:"azimuth_deg"` + Name string `yaml:"name,omitempty" json:"name,omitempty"` + KWp float64 `yaml:"kwp" json:"kwp"` + TiltDeg *float64 `yaml:"tilt_deg" json:"tilt_deg"` + AzimuthDeg *float64 `yaml:"azimuth_deg" json:"azimuth_deg"` +} + +// CompleteGeometry returns one usable PV plane. Tilt and azimuth are +// pointers so an omitted field cannot be confused with a valid 0° value. +// Invalid or partial entries are intentionally not fatal: callers use the +// flat forecast path when no complete plane remains. +func (a PVArray) CompleteGeometry() (tiltDeg, azimuthDeg, kWp float64, ok bool) { + if a.KWp <= 0 || math.IsNaN(a.KWp) || math.IsInf(a.KWp, 0) || + a.TiltDeg == nil || a.AzimuthDeg == nil { + return 0, 0, 0, false + } + tiltDeg, azimuthDeg = *a.TiltDeg, *a.AzimuthDeg + if math.IsNaN(tiltDeg) || math.IsInf(tiltDeg, 0) || tiltDeg < 0 || tiltDeg > 90 || + math.IsNaN(azimuthDeg) || math.IsInf(azimuthDeg, 0) || azimuthDeg < 0 || azimuthDeg > 360 { + return 0, 0, 0, false + } + return tiltDeg, azimuthDeg, a.KWp, true } // Battery is per-battery overrides (keyed by driver name in the top-level map). diff --git a/go/internal/config/config_test.go b/go/internal/config/config_test.go index 88f3e4f0..eea190bd 100644 --- a/go/internal/config/config_test.go +++ b/go/internal/config/config_test.go @@ -574,6 +574,42 @@ batteries: } } +func TestPVArrayGeometryDistinguishesMissingFromZero(t *testing.T) { + yaml := minimalYAML + ` +weather: + provider: open_meteo + latitude: 59.3293 + longitude: 18.0686 + pv_arrays: + - name: partial + kwp: 10 + tilt_deg: 35 + - name: north flat + kwp: 5 + tilt_deg: 0 + azimuth_deg: 0 +` + c, err := Parse([]byte(yaml), "/tmp") + if err != nil { + t.Fatal(err) + } + if c.Weather == nil || len(c.Weather.PVArrays) != 2 { + t.Fatalf("weather arrays missing: %+v", c.Weather) + } + partial := c.Weather.PVArrays[0] + if partial.AzimuthDeg != nil { + t.Fatalf("omitted azimuth should remain nil, got %v", *partial.AzimuthDeg) + } + if _, _, _, ok := partial.CompleteGeometry(); ok { + t.Fatal("partial geometry must not be treated as a north-facing array") + } + northFlat := c.Weather.PVArrays[1] + tilt, azimuth, kwp, ok := northFlat.CompleteGeometry() + if !ok || tilt != 0 || azimuth != 0 || kwp != 5 { + t.Fatalf("explicit zero geometry should remain valid: tilt=%v azimuth=%v kwp=%v ok=%v", tilt, azimuth, kwp, ok) + } +} + func TestSiteMeterDriverReturnsName(t *testing.T) { c, err := Parse([]byte(minimalYAML), ".") if err != nil { diff --git a/go/internal/config/restart_required_test.go b/go/internal/config/restart_required_test.go index 2db78e1e..08cbe823 100644 --- a/go/internal/config/restart_required_test.go +++ b/go/internal/config/restart_required_test.go @@ -89,8 +89,9 @@ func TestRestartRequiredFor_BootSections(t *testing.T) { c.Weather = &Weather{Provider: "open_meteo", Latitude: 59, Longitude: 18} }, "weather"}, {"weather pv_arrays added", func(c *Config) { + tilt, azimuth := 30.0, 180.0 c.Weather = &Weather{Provider: "met_no", Latitude: 59, Longitude: 18, - PVArrays: []PVArray{{KWp: 5, TiltDeg: 30, AzimuthDeg: 180}}} + PVArrays: []PVArray{{KWp: 5, TiltDeg: &tilt, AzimuthDeg: &azimuth}}} }, "weather"}, {"weather heating coefficient", func(c *Config) { c.Weather.HeatingWPerDegC = 250 diff --git a/go/internal/forecast/forecast.go b/go/internal/forecast/forecast.go index 08a55efb..22da4c5b 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{} @@ -273,9 +280,9 @@ func FromConfig(cfg *config.Weather, ratedPVW float64, st *state.Store, userAgen // their geometry just because the config model grew. var arrays []Array for _, a := range cfg.PVArrays { - arrays = append(arrays, Array{ - TiltDeg: a.TiltDeg, AzimuthDeg: a.AzimuthDeg, KWp: a.KWp, - }) + if converted, ok := arrayFromConfig(a); ok { + arrays = append(arrays, converted) + } } if len(arrays) == 0 { arrays = append(arrays, Array{ @@ -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 converted, ok := arrayFromConfig(a); ok { + arrays = append(arrays, converted) + } + } return &Service{ Provider: p, Store: st, Lat: cfg.Latitude, Lon: cfg.Longitude, RatedPVW: ratedPVW, + Arrays: arrays, stop: make(chan struct{}), done: make(chan struct{}), } @@ -333,6 +351,18 @@ func (s *Service) fetchAndStore(ctx context.Context) { nowMs := time.Now().UnixMilli() points := make([]state.ForecastPoint, 0, len(rows)) for _, r := range rows { + // A negative irradiance is not physical; retain the row with a + // zero signal. Non-finite values are invalid provider data and must + // not reach SQLite, where they can become NULL silently. + var solarWm2 *float64 + if r.SolarWm2 != nil { + ghi, ok := normalizeIrradiance(*r.SolarWm2) + if !ok { + slog.Warn("forecast row skipped", "reason", "non-finite irradiance", "provider", s.Provider.Name(), "slot", r.HourStart) + continue + } + solarWm2 = &ghi + } // Pick the most direct PV signal the provider gave us. Forecast.Solar // returns site-calibrated watts directly; Open-Meteo returns shortwave // radiation we turn into watts via rated × W/m²/1000; met.no only has @@ -341,18 +371,30 @@ func (s *Service) fetchAndStore(ctx context.Context) { switch { case r.PVWEstimated != nil: pvW = *r.PVWEstimated - case r.SolarWm2 != nil && s.RatedPVW > 0: - pvW = s.RatedPVW * (*r.SolarWm2) / 1000.0 + case solarWm2 != nil: + var ok bool + pvW, ok = pvWFromGHI(s.Lat, s.Lon, r.HourStart, *solarWm2, s.RatedPVW, s.Arrays) + if !ok { + slog.Warn("forecast row skipped", "reason", "non-finite irradiance", "provider", s.Provider.Name(), "slot", r.HourStart) + continue + } default: pvW = EstimatePVW(s.Lat, s.Lon, r.HourStart, r.CloudCoverPct, s.RatedPVW) } + if math.IsNaN(pvW) || math.IsInf(pvW, 0) { + slog.Warn("forecast row skipped", "reason", "non-finite PV estimate", "provider", s.Provider.Name(), "slot", r.HourStart) + continue + } + if pvW < 0 { + pvW = 0 + } pvPtr := &pvW points = append(points, state.ForecastPoint{ SlotTsMs: r.HourStart.UnixMilli(), SlotLenMin: 60, CloudCoverPct: r.CloudCoverPct, TempC: r.TempC, - SolarWm2: r.SolarWm2, + SolarWm2: solarWm2, PVWEstimated: pvPtr, Source: s.Provider.Name(), FetchedAtMs: nowMs, @@ -365,6 +407,61 @@ func (s *Service) fetchAndStore(ctx context.Context) { slog.Info("forecast fetched", "count", len(points), "provider", s.Provider.Name()) } +func arrayFromConfig(a config.PVArray) (Array, bool) { + tiltDeg, azimuthDeg, kWp, ok := a.CompleteGeometry() + if !ok { + return Array{}, false + } + return Array{TiltDeg: tiltDeg, AzimuthDeg: azimuthDeg, KWp: kWp}, true +} + +func normalizeIrradiance(ghiWm2 float64) (float64, bool) { + if math.IsNaN(ghiWm2) || math.IsInf(ghiWm2, 0) { + return 0, false + } + if ghiWm2 < 0 { + return 0, true + } + return ghiWm2, true +} + +func pvWFromGHI(lat, lon float64, t time.Time, ghiWm2, ratedPVW float64, arrays []Array) (float64, bool) { + ghiWm2, ok := normalizeIrradiance(ghiWm2) + if !ok { + return 0, false + } + if len(arrays) > 0 { + return poaPVWattsFromGHI(lat, lon, t, ghiWm2, arrays), true + } + if ratedPVW <= 0 { + return 0, true + } + return ratedPVW * ghiWm2 / 1000.0, true +} + +// 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 { + ghiWm2, ok := normalizeIrradiance(ghiWm2) + if !ok { + return 0 + } + 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..0517dea7 100644 --- a/go/internal/forecast/forecast_test.go +++ b/go/internal/forecast/forecast_test.go @@ -14,6 +14,25 @@ import ( "github.com/srcfl/ftw/go/internal/state" ) +type staticForecastProvider struct { + rows []RawForecast +} + +func (p staticForecastProvider) Name() string { return "static" } + +func (p staticForecastProvider) Fetch(context.Context, float64, float64) ([]RawForecast, error) { + return p.rows, nil +} + +func testPVArray(name string, kwp, tiltDeg, azimuthDeg float64) config.PVArray { + return config.PVArray{ + Name: name, + KWp: kwp, + TiltDeg: &tiltDeg, + AzimuthDeg: &azimuthDeg, + } +} + // ---- Clear-sky model sanity ---- func TestClearSkyIsZeroAtMidnight(t *testing.T) { @@ -247,3 +266,186 @@ 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{ + testPVArray("south", 6, 35, 180), + testPVArray("east", 4, 30, 90), + testPVArray("empty", 0, 10, 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) + } +} + +func TestFromConfigSkipsPartialArrayGeometry(t *testing.T) { + st, err := state.Open(filepath.Join(t.TempDir(), "t.db")) + if err != nil { + t.Fatal(err) + } + defer st.Close() + tilt := 35.0 + cfg := &config.Weather{ + Provider: "open_meteo", Latitude: 59.3293, Longitude: 18.0686, + PVArrays: []config.PVArray{ + {Name: "missing azimuth", KWp: 10, TiltDeg: &tilt}, + testPVArray("Stockholm south", 6, 35, 180), + }, + } + s := FromConfig(cfg, 16000, st, "ua") + if s == nil { + t.Fatal("expected service") + } + if len(s.Arrays) != 1 { + t.Fatalf("expected only complete Stockholm geometry, got %d arrays: %+v", len(s.Arrays), s.Arrays) + } + if s.Arrays[0].AzimuthDeg != 180 || s.Arrays[0].KWp != 6 { + t.Fatalf("unexpected complete geometry: %+v", s.Arrays[0]) + } +} + +// ---- 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) + } +} + +func TestServiceGHIPhysicalBounds(t *testing.T) { + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + cases := []struct { + name string + ghi float64 + valid bool + wantSolar float64 + }{ + {name: "negative", ghi: -100, valid: true, wantSolar: 0}, + {name: "zero", ghi: 0, valid: true, wantSolar: 0}, + {name: "nan", ghi: math.NaN(), valid: false}, + {name: "positive infinity", ghi: math.Inf(1), valid: false}, + } + for _, withArrays := range []bool{false, true} { + path := "without arrays" + if withArrays { + path = "with arrays" + } + for _, tc := range cases { + t.Run(path+"/"+tc.name, func(t *testing.T) { + st, err := state.Open(filepath.Join(t.TempDir(), "state.db")) + if err != nil { + t.Fatal(err) + } + defer st.Close() + ghi := tc.ghi + s := &Service{ + Provider: staticForecastProvider{rows: []RawForecast{{HourStart: tt, SolarWm2: &ghi}}}, + Store: st, + Lat: 59.3293, + Lon: 18.0686, + RatedPVW: 10000, + Arrays: nil, + } + if withArrays { + s.Arrays = []Array{{TiltDeg: 35, AzimuthDeg: 180, KWp: 10}} + } + s.fetchAndStore(context.Background()) + + rows, err := st.LoadForecasts(tt.UnixMilli(), tt.Add(time.Hour).UnixMilli()) + if err != nil { + t.Fatal(err) + } + if !tc.valid { + if len(rows) != 0 { + t.Fatalf("non-finite irradiance should omit the row, got %+v", rows) + } + return + } + if len(rows) != 1 || rows[0].PVWEstimated == nil || rows[0].SolarWm2 == nil { + t.Fatalf("expected one finite forecast row, got %+v", rows) + } + if got := *rows[0].PVWEstimated; got != 0 { + t.Errorf("PV estimate from %s irradiance = %.2f, want 0", tc.name, got) + } + if got := *rows[0].SolarWm2; got != tc.wantSolar { + t.Errorf("stored irradiance = %.2f, want %.2f", got, tc.wantSolar) + } + }) + } + } +} + +// 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/forecast/open_meteo.go b/go/internal/forecast/open_meteo.go index 5eba47e0..ba9d9ce2 100644 --- a/go/internal/forecast/open_meteo.go +++ b/go/internal/forecast/open_meteo.go @@ -17,10 +17,11 @@ import ( // than cloud_area_fraction alone because "49% cloud in the sky" can // mean anything from full sun to overcast depending on which part of // the sky the Sun is behind — radiation measures what actually reaches -// a horizontal surface. The downstream PV derivation simplifies to -// `rated × (W/m² / 1000)` with a single calibration coefficient per -// site (orientation + soiling), not the brittle `(1 - cloud)^1.5` -// curve used when only cloud fraction is available. +// a horizontal surface. The downstream PV derivation projects this +// irradiance onto each complete `weather.pv_arrays` plane when present, +// and uses the safe flat `rated × (W/m² / 1000)` estimate otherwise. +// Forecast.Solar remains site-calibrated server-side; Open-Meteo uses the +// configured arrays directly rather than ignoring them. type OpenMeteoProvider struct { Client *http.Client BaseURL string diff --git a/go/internal/sunpos/sunpos.go b/go/internal/sunpos/sunpos.go index 8ce95cc3..fd5f81a5 100644 --- a/go/internal/sunpos/sunpos.go +++ b/go/internal/sunpos/sunpos.go @@ -125,35 +125,120 @@ 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 { + if math.IsNaN(ghi) || math.IsInf(ghi, 0) || ghi <= 0 { + return 0 + } + sun := At(t, lat, lon) + if sun.ZenithDeg >= 90 { + 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..392db3ad 100644 --- a/go/internal/sunpos/sunpos_test.go +++ b/go/internal/sunpos/sunpos_test.go @@ -66,3 +66,111 @@ 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) + } +} + +func TestPOAFromGHINonFiniteOrNonPositiveGHIIsZero(t *testing.T) { + when := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + tests := []struct { + name string + ghi float64 + }{ + {name: "NaN", ghi: math.NaN()}, + {name: "positive infinity", ghi: math.Inf(1)}, + {name: "negative infinity", ghi: math.Inf(-1)}, + {name: "negative", ghi: -1}, + {name: "zero", ghi: 0}, + } + for _, tc := range tests { + t.Run(tc.name, func(t *testing.T) { + if got := POAFromGHI(when, 59.33, 18.07, tc.ghi, 35, 180); got != 0 { + t.Errorf("POAFromGHI(%v) = %v, want 0", tc.ghi, got) + } + }) + } +} + +// 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) + } +} diff --git a/web/settings/tabs/weather.js b/web/settings/tabs/weather.js index 7d90365d..85af1ad5 100644 --- a/web/settings/tabs/weather.js +++ b/web/settings/tabs/weather.js @@ -181,8 +181,8 @@ field("API key (OpenWeather only)", "weather.api_key", "text", "") + '' + '
PV arrays ' + help( - 'Optional. If set, forecast_solar uses these per-plane values to produce a site-calibrated forecast. ' + - 'Leave empty to let the model learn your orientation from telemetry — predictions are fine after a few varied days.') + '' + + 'Optional. Open-Meteo uses these per-plane values to project shortwave radiation onto each array. ' + + 'Forecast.Solar uses them for its site-calibrated forecast. Leave empty for the safe flat estimate or the provider default.') + '' + '
' + '' + '

' +