Add day/night estimation for NWS observations based on solar elevation and update observation model.

internal/normalizers/nws/day_night.go

@@ -0,0 +1,74 @@
+// FILE: internal/normalizers/nws/day_night.go
+package nws
+
+import (
+	"math"
+	"time"
+)
+
+const (
+	degToRad = math.Pi / 180.0
+	radToDeg = 180.0 / math.Pi
+)
+
+func observationLatLon(coords []float64) (lat *float64, lon *float64) {
+	if len(coords) < 2 {
+		return nil, nil
+	}
+	latVal := coords[1]
+	lonVal := coords[0]
+	return &latVal, &lonVal
+}
+
+// isDayFromLatLonTime estimates day/night from solar elevation.
+// Uses a refraction-adjusted cutoff (-0.833 degrees).
+func isDayFromLatLonTime(lat, lon float64, ts time.Time) *bool {
+	if ts.IsZero() || math.IsNaN(lat) || math.IsNaN(lon) || math.IsInf(lat, 0) || math.IsInf(lon, 0) {
+		return nil
+	}
+	if lat < -90.0 || lat > 90.0 || lon < -180.0 || lon > 180.0 {
+		return nil
+	}
+
+	t := ts.UTC()
+	day := float64(t.YearDay())
+	hour := float64(t.Hour())
+	min := float64(t.Minute())
+	sec := float64(t.Second()) + float64(t.Nanosecond())/1e9
+	utcHours := hour + min/60.0 + sec/3600.0
+
+	gamma := 2.0 * math.Pi / 365.0 * (day - 1.0 + (utcHours-12.0)/24.0)
+
+	eqtime := 229.18 * (0.000075 + 0.001868*math.Cos(gamma) - 0.032077*math.Sin(gamma) -
+		0.014615*math.Cos(2.0*gamma) - 0.040849*math.Sin(2.0*gamma))
+
+	decl := 0.006918 - 0.399912*math.Cos(gamma) + 0.070257*math.Sin(gamma) -
+		0.006758*math.Cos(2.0*gamma) + 0.000907*math.Sin(2.0*gamma) -
+		0.002697*math.Cos(3.0*gamma) + 0.00148*math.Sin(3.0*gamma)
+
+	timeOffset := eqtime + 4.0*lon
+	trueSolarMinutes := hour*60.0 + min + sec/60.0 + timeOffset
+	for trueSolarMinutes < 0 {
+		trueSolarMinutes += 1440.0
+	}
+	for trueSolarMinutes >= 1440.0 {
+		trueSolarMinutes -= 1440.0
+	}
+
+	hourAngleDeg := trueSolarMinutes/4.0 - 180.0
+	ha := hourAngleDeg * degToRad
+	latRad := lat * degToRad
+
+	cosZenith := math.Sin(latRad)*math.Sin(decl) + math.Cos(latRad)*math.Cos(decl)*math.Cos(ha)
+	if cosZenith > 1.0 {
+		cosZenith = 1.0
+	} else if cosZenith < -1.0 {
+		cosZenith = -1.0
+	}
+
+	zenith := math.Acos(cosZenith)
+	elevation := 90.0 - zenith*radToDeg
+
+	isDay := elevation > -0.833
+	return &isDay
+}