roam-location/weather/condition.go
Adriano Caloiaro 567845f88e
fix(weather): apply Beer-Lambert extinction to clear-sky model
The simple atmosphericCeiling*sin(elevation) model overestimates clear-sky
radiation at low sun angles because it ignores the longer atmospheric path
light travels near the horizon. At 5° elevation the path is ~11x vertical,
making the real clear-sky ceiling ~26 W/m² rather than the 83 W/m² the
naive model predicts — causing sunny low-angle readings to be misclassified
as mostly cloudy.

Replace with the Beer-Lambert correction:
  expectedClearSky = 950 * exp(-opticalDepth / sin(elevation)) * sin(elevation)

opticalDepth=0.1 (clear atmosphere) produces accurate values across all
angles: ~26 W/m² at 5°, ~93 W/m² at 10°, ~715 W/m² at solar noon.
2026-04-12 07:39:31 -06:00

109 lines
3.7 KiB
Go

package weather
import (
"math"
"time"
)
// Condition is an inferred weather condition with a human label and emoji.
type Condition struct {
Label string `json:"label"`
Emoji string `json:"emoji"`
}
// atmosphericCeiling is the extraterrestrial solar irradiance at sea level (W/m²)
// with the sun directly overhead.
const atmosphericCeiling = 950.0
// opticalDepth is the vertical atmospheric optical depth used in the Beer-Lambert
// extinction model. At low sun angles the light path through the atmosphere is
// 1/sin(elevation) times longer than vertical, so the effective transmittance
// drops sharply. A value of 0.1 represents a clear sky and produces realistic
// clear-sky readings across all elevation angles (e.g. ~26 W/m² at 5°, ~715 W/m²
// at solar noon for mid-latitudes in spring).
const opticalDepth = 0.1
// solarPosition returns the solar elevation angle (radians) and hour angle
// (degrees) for the given latitude, longitude, and UTC time.
//
// Elevation > 0 means the sun is above the horizon.
// Hour angle < 0 means the sun is east of the meridian (morning);
// > 0 means west of the meridian (afternoon).
func solarPosition(lat, lon float64, t time.Time) (elevationRad, hourAngleDeg float64) {
decRad := 23.45 * math.Pi / 180 * math.Sin(2*math.Pi/365*float64(t.YearDay()-81))
utcHours := float64(t.Hour()) + float64(t.Minute())/60.0
// math.Remainder normalises to [-180, 180], handling midnight crossings cleanly.
hourAngleDeg = math.Remainder((utcHours-(12.0-lon/15.0))*15, 360)
hourAngleRad := hourAngleDeg * math.Pi / 180
latRad := lat * math.Pi / 180
sinElev := math.Sin(latRad)*math.Sin(decRad) + math.Cos(latRad)*math.Cos(decRad)*math.Cos(hourAngleRad)
elevationRad = math.Asin(math.Max(-1, math.Min(1, sinElev)))
return
}
// Infer derives a weather condition from sensor readings and position.
// Solar radiation is normalised against the clear-sky expectation for the
// actual solar elevation angle, so conditions are accurate at any time of day
// without any hardcoded sunrise or sunset times.
func Infer(d Data, lat, lon float64) *Condition {
solar := 0.0
if d.SolarRadiation != nil {
solar = *d.SolarRadiation
}
humidity := 0.0
if d.Humidity != nil {
humidity = *d.Humidity
}
// Rain takes priority.
if d.HourlyRainIn != nil && *d.HourlyRainIn > 0.02 {
if *d.HourlyRainIn > 0.1 {
return &Condition{"Heavy rain", "⛈️"}
}
return &Condition{"Rainy", "🌧️"}
}
// Near-rain: very high humidity with no solar.
if humidity >= 95 && solar < 10 {
return &Condition{"Foggy", "🌫️"}
}
elevationRad, hourAngleDeg := solarPosition(lat, lon, time.Now().UTC())
isMorning := hourAngleDeg < 0 // sun east of meridian
if elevationRad <= 0 {
// Sun is below the horizon.
if solar > 0 {
// Sensor still catching scattered twilight.
if isMorning {
return &Condition{"Dawn", "🌅"}
}
return &Condition{"Dusk", "🌇"}
}
if humidity > 80 {
return &Condition{"Cloudy night", "☁️"}
}
return &Condition{"Clear night", "🌙"}
}
// Sun is above the horizon. Normalise against the clear-sky radiation
// expected at this elevation angle, corrected for atmospheric extinction
// (Beer-Lambert): at low angles the light path is ~1/sin(elev) times
// longer than vertical, so the clear-sky ceiling drops sharply.
expectedClearSky := atmosphericCeiling * math.Exp(-opticalDepth/math.Sin(elevationRad)) * math.Sin(elevationRad)
if solar < 10 {
return &Condition{"Overcast", "☁️"}
}
normalized := math.Min(solar/expectedClearSky, 1.0)
if normalized > 0.85 {
return &Condition{"Sunny", "🌞"}
}
if normalized > 0.55 {
return &Condition{"Partly cloudy", "⛅"}
}
if normalized > 0.15 {
return &Condition{"Mostly cloudy", "🌥️"}
}
return &Condition{"Overcast", "☁️"}
}