From b7fd7046ffe814515de99b9d1d175134150f299d Mon Sep 17 00:00:00 2001 From: "a.antonov" Date: Sat, 30 May 2026 06:38:38 +0900 Subject: [PATCH] feat: move stuff to numerics --- internal/engine/models.go | 24 ++++--------- internal/engine/registry.go | 30 ++++++++++------- internal/numerics/atmosphere.go | 13 +++++++ internal/numerics/motion_test.go | 58 ++++++++++++++++++++++++++++++++ internal/numerics/vec.go | 23 +++++++++++++ 5 files changed, 119 insertions(+), 29 deletions(-) create mode 100644 internal/numerics/motion_test.go diff --git a/internal/engine/models.go b/internal/engine/models.go index cdfdf1a..9a738d8 100644 --- a/internal/engine/models.go +++ b/internal/engine/models.go @@ -1,7 +1,6 @@ package engine import ( - "math" "sort" "predictor-refactored/internal/numerics" @@ -19,10 +18,7 @@ func Sum(models ...Model) Model { return func(t float64, s State) State { var sum State for _, m := range models { - d := m(t, s) - sum.Lat += d.Lat - sum.Lng += d.Lng - sum.Altitude += d.Altitude + sum = numerics.AddGeo(sum, m(t, s)) } return sum } @@ -44,9 +40,8 @@ func ConstantRate(rate float64) Model { // // using the NASA atmosphere model for rho. Equivalent to Tawhiri's drag_descent. func ParachuteDescent(seaLevelRate float64) Model { - k := seaLevelRate * 1.1045 return func(_ float64, s State) State { - return State{Altitude: -k / math.Sqrt(numerics.NasaDensity(s.Altitude))} + return State{Altitude: numerics.DragTerminalVelocity(seaLevelRate, s.Altitude)} } } @@ -79,15 +74,13 @@ func Piecewise(segments []RateSegment) Model { } // WindTransport returns a model that moves laterally at the wind velocity -// sampled from field. Vertical component is zero. Wind components in m/s -// are converted to deg/s on Earth's surface using R = 6371009 m. +// sampled from field. The vertical component is zero. Sampling and the +// non-fatal "above_model" event live here (orchestration); the m/s → deg/s +// conversion is numerics.WindToGeoRate. // // If events is non-nil, an "above_model" event is emitted whenever the // wind field reports altitude above the highest pressure level. func WindTransport(field weather.WindField, events *EventSink) Model { - const earthR = 6371009.0 - const piOver180 = math.Pi / 180.0 - const degPerRad = 180.0 / math.Pi return func(t float64, s State) State { sample, err := field.Wind(t, s.Lat, s.Lng, s.Altitude) if err != nil { @@ -97,10 +90,7 @@ func WindTransport(field weather.WindField, events *EventSink) Model { events.Emit("above_model", t, s, "altitude exceeded the highest pressure level of the wind dataset; samples extrapolated") } - r := earthR + s.Altitude - return State{ - Lat: degPerRad * sample.V / r, - Lng: degPerRad * sample.U / (r * math.Cos(s.Lat*piOver180)), - } + dLat, dLng := numerics.WindToGeoRate(sample.U, sample.V, s.Lat, s.Altitude) + return State{Lat: dLat, Lng: dLng} } } diff --git a/internal/engine/registry.go b/internal/engine/registry.go index db6d53f..27a7a9f 100644 --- a/internal/engine/registry.go +++ b/internal/engine/registry.go @@ -241,25 +241,31 @@ func buildPiecewise(spec ModelSpec, deps BuildDeps) (BuiltModel, error) { }, nil } -// resolveSegments converts spec segments to engine.RateSegment using the -// stage context to resolve relative references. +// resolveSegments converts spec segments to engine.RateSegment, turning each +// segment's reference-relative Until into an absolute UNIX time. References +// are validated by buildPiecewise, so an unrecognised one here is treated as +// absolute rather than re-erroring. func resolveSegments(in []PiecewiseSegmentSpec, ctx StageContext) []RateSegment { out := make([]RateSegment, 0, len(in)) for _, s := range in { - var until float64 - switch s.Reference { - case "", "absolute": - until = s.Until - case "profile_start": - until = ctx.ProfileStart + s.Until - case "propagator_start": - until = ctx.PropagatorStart + s.Until - } - out = append(out, RateSegment{Until: until, Rate: s.Rate}) + out = append(out, RateSegment{Until: segmentBase(s.Reference, ctx) + s.Until, Rate: s.Rate}) } return out } +// segmentBase returns the absolute time a piecewise segment's Until is +// measured from, per its reference. +func segmentBase(reference string, ctx StageContext) float64 { + switch reference { + case "profile_start": + return ctx.ProfileStart + case "propagator_start": + return ctx.PropagatorStart + default: // "", "absolute" + return 0 + } +} + // maybeAddWind sums a WindTransport model into base when the spec asks for it. func maybeAddWind(base Model, includeWind bool, deps BuildDeps) Model { if !includeWind { diff --git a/internal/numerics/atmosphere.go b/internal/numerics/atmosphere.go index fa6bb31..9002360 100644 --- a/internal/numerics/atmosphere.go +++ b/internal/numerics/atmosphere.go @@ -23,3 +23,16 @@ func NasaDensity(alt float64) float64 { } return pressure / (0.2869 * (temp + 273.1)) } + +// DragTerminalVelocity returns the vertical velocity (m/s, negative = downward) +// of a parachute descent at the given altitude. seaLevelRate is the descent +// speed at sea level (positive m/s); the rate grows with altitude as the +// thinner air provides less drag: +// +// v = -k / sqrt(rho(alt)), k = seaLevelRate * 1.1045 +// +// Matches Tawhiri's drag_descent. +func DragTerminalVelocity(seaLevelRate, alt float64) float64 { + k := seaLevelRate * 1.1045 + return -k / math.Sqrt(NasaDensity(alt)) +} diff --git a/internal/numerics/motion_test.go b/internal/numerics/motion_test.go new file mode 100644 index 0000000..3bc8051 --- /dev/null +++ b/internal/numerics/motion_test.go @@ -0,0 +1,58 @@ +package numerics + +import ( + "math" + "testing" +) + +func TestAddGeo(t *testing.T) { + // Rates sum component-wise with no longitude wrapping. + got := AddGeo(GeoVec{Lat: 1, Lng: 350, Altitude: 2}, GeoVec{Lat: 3, Lng: 20, Altitude: 4}) + want := GeoVec{Lat: 4, Lng: 370, Altitude: 6} + if got != want { + t.Errorf("AddGeo = %+v, want %+v (no wrap on rates)", got, want) + } +} + +func TestWindToGeoRate(t *testing.T) { + // Pure eastward 10 m/s at the equator, sea level. + dLat, dLng := WindToGeoRate(10, 0, 0, 0) + wantLng := (180.0 / math.Pi) * 10.0 / EarthRadius + if math.Abs(dLat) > 1e-15 { + t.Errorf("dLat = %v, want 0", dLat) + } + if math.Abs(dLng-wantLng) > 1e-15 { + t.Errorf("dLng = %v, want %v", dLng, wantLng) + } + + // Northward 5 m/s at 60°N: dLat independent of longitude scaling. + dLat, _ = WindToGeoRate(0, 5, 60, 0) + wantLat := (180.0 / math.Pi) * 5.0 / EarthRadius + if math.Abs(dLat-wantLat) > 1e-15 { + t.Errorf("dLat at 60N = %v, want %v", dLat, wantLat) + } + + // cos(lat) factor makes eastward motion span more degrees nearer the poles. + _, dLngEq := WindToGeoRate(10, 0, 0, 0) + _, dLng60 := WindToGeoRate(10, 0, 60, 0) + if dLng60 <= dLngEq { + t.Errorf("eastward deg/s should grow with latitude: eq=%v 60N=%v", dLngEq, dLng60) + } +} + +func TestDragTerminalVelocity(t *testing.T) { + // Descent is downward (negative) and faster (more negative) at altitude + // where the air is thinner. + sea := DragTerminalVelocity(5, 0) + high := DragTerminalVelocity(5, 20000) + if sea >= 0 { + t.Errorf("sea-level rate = %v, want negative (downward)", sea) + } + if high >= sea { + t.Errorf("expected faster descent at altitude: sea=%v high=%v", sea, high) + } + // Sanity: at sea level rho≈1.225, so v ≈ -5*1.1045/sqrt(1.225) ≈ -4.99 m/s. + if math.Abs(sea-(-5*1.1045/math.Sqrt(NasaDensity(0)))) > 1e-12 { + t.Errorf("sea-level formula mismatch: %v", sea) + } +} diff --git a/internal/numerics/vec.go b/internal/numerics/vec.go index dc120b8..08d22a5 100644 --- a/internal/numerics/vec.go +++ b/internal/numerics/vec.go @@ -64,3 +64,26 @@ func LngLerp(a, b, l float64) float64 { func Lerp(a, b, l float64) float64 { return (1-l)*a + l*b } + +// AddGeo returns the component-wise sum a+b without longitude wrapping. Use it +// to combine derivative (rate) vectors — rates accumulate linearly, unlike +// positions, which wrap via GeoAdd. +func AddGeo(a, b GeoVec) GeoVec { + return GeoVec{Lat: a.Lat + b.Lat, Lng: a.Lng + b.Lng, Altitude: a.Altitude + b.Altitude} +} + +// EarthRadius is the spherical Earth radius (metres) used for horizontal +// motion, matching the reference Tawhiri implementation. +const EarthRadius = 6371009.0 + +// WindToGeoRate converts eastward (u) and northward (v) wind in m/s at the +// given latitude (deg) and altitude (m) into the geographic rate in deg/s on a +// spherical Earth. The returned dLng diverges near the poles as cos(lat) → 0. +func WindToGeoRate(u, v, lat, alt float64) (dLat, dLng float64) { + const degPerRad = 180.0 / math.Pi + const piOver180 = math.Pi / 180.0 + r := EarthRadius + alt + dLat = degPerRad * v / r + dLng = degPerRad * u / (r * math.Cos(lat*piOver180)) + return dLat, dLng +}