diff --git a/docs/numerics.tex b/docs/numerics.tex
index 16aaacd..c926a5f 100644
--- a/docs/numerics.tex
+++ b/docs/numerics.tex
@@ -74,25 +74,66 @@ The contribution at time $t$ is
 \]
 
 \paragraph{Wind transport.} The horizontal contribution from sampling the
-loaded wind field $W$:
+loaded wind field $W$ is the wind itself:
 \[
-  \mathbf{F}_{\text{wind}}(t, \mathbf{s}) = \Bigl(
-    \frac{180}{\pi}\,\frac{v}{R + h},\;\;
-    \frac{180}{\pi}\,\frac{u}{(R + h)\cos\bigl(\varphi\,\pi/180\bigr)},\;\;
-    0
-  \Bigr),
+  \mathbf{F}_{\text{wind}}(t, \mathbf{s}) = (u,\; v,\; 0),
+  \qquad (u, v) = W(t, \varphi, \lambda, h),
 \]
-where $(u, v) = W(t, \varphi, \lambda, h)$ are the eastward and northward
-wind components in metres per second, and $R = 6{,}371{,}009$~m is the
-spherical Earth radius. The implementation lives in
+in metres per second east and north. No conversion is performed.
+
+Earlier revisions converted this to degrees per second, which introduced a
+$1/\cos\varphi$ factor in longitude. That factor diverges at the poles: for
+$u = 10$~m/s it grows from $8.95\times10^{-5}$~deg/s at the equator to
+$0.513$~deg/s at $\varphi = 89.99^\circ$ and $1.46\times10^{12}$~deg/s at
+$\varphi = 90^\circ$ (large but finite, since $\cos(\pi/2)$ evaluates to
+$6.12\times10^{-17}$ in double precision rather than to zero). Against a fixed
+step this made the integrator meaningless near the poles. The factor is now
+absent from the formulation rather than guarded against; see
+section~\ref{sec:geostep}. The implementation lives in
 \verb|engine.WindTransport| (\verb|engine/models.go|).
 
-\paragraph{Coordinate system.} The model is a spherical Earth in
-plate-carrée (latitude/longitude/altitude) coordinates. This matches the
-reference Tawhiri predictor exactly and is necessary for bit-identical
-back-to-back testing. A WGS84/ECEF variant is planned but deferred: it
-would require converting U/V wind components from the GFS sphere model
-to the ellipsoid, which is not a trivial coordinate transform.
+\paragraph{Coordinate system.} The model is a spherical Earth. State is
+still carried as $(\varphi, \lambda, h)$ in degrees and metres, because
+constraints, path recording and the REST API all speak latitude and
+longitude --- but motion is \emph{not} integrated in those coordinates.
+Displacement is applied by rotating the position vector along a great
+circle (section~\ref{sec:geostep}), so no longitude derivative is ever
+formed and there is no coordinate singularity at the poles. Latitude also
+cannot leave $[-90, 90]$, since it is read back from a unit vector instead
+of accumulated.
+
+This is a deliberate departure from the reference Tawhiri predictor, which
+integrates in plate-carrée coordinates. Outputs are therefore no longer
+bit-identical to it. The measured cost is small: on GFS
+\verb|2026-08-03T00:00:00Z|, launch $89^\circ$N $68^\circ$E, burst
+25~km --- a 114~km flight --- the two integrators differ by at most
+$1.608$~m along the track and $0.058$~m at the landing point. That is the
+order of the truncation error, and far below the representation error of
+the wind data itself. Back-to-back testing against Tawhiri is retained as
+an agreement-within-tolerance check rather than an equality check, run with
+\verb|cmd/compare-tawhiri| and its \verb|-align-dataset| flag (on by default),
+which asks the hosted service to use the local predictor's GFS run --- without
+it the two sides silently compare different weather. Note the hosted service
+retains only recent runs, so alignment fails once the local dataset ages out.
+
+Measured on GFS \verb|2026-08-03T06:00:00Z|, burst 25~km:
+\begin{center}
+\begin{tabular}{lrrrr}
+launch & burst $\Delta$ & landing $\Delta$ & apex alt $\Delta$ & land alt $\Delta$ \\
+\hline
+$52.2^\circ$N $0.1^\circ$E & 1~m & 60~m  & 0~m & 47~m \\
+$89^\circ$N $68^\circ$E    & 2~m & 0~m   & 0~m & 0~m  \\
+\end{tabular}
+\end{center}
+The polar launch agrees exactly. The 60~m at mid-latitude is not integrator
+error: it follows from the 47~m difference in termination altitude, because the
+reference has the ruaumoko elevation dataset and terminates on terrain while
+this deployment has none and terminates at sea level. At $89^\circ$N the
+surface is sea ice, so sea level is the terrain and the difference vanishes.
+
+A WGS84 ellipsoid variant remains deferred: it would require converting
+U/V wind components from the GFS sphere model to the ellipsoid, which is
+not a trivial coordinate transform.
 
 % =========================================================================
 \section{Profiles and propagators}
@@ -181,12 +222,29 @@ $x_i = \ell + i \cdot s$ for $i = 0, 1, \ldots, N - 1$, parameterised by
 the left edge $\ell$, the step $s > 0$, and the point count $N$.
 
 Given a query $v$, the \emph{bracket} is the pair $(i_0, i_1)$ with
-$x_{i_0} \le v < x_{i_1}$ and the dimensionless position
+$x_{i_0} \le v \le x_{i_1}$ and the dimensionless position
 \[
-  f = \frac{v - x_{i_0}}{s} \in [0, 1).
+  f = \frac{v - x_{i_0}}{s} \in [0, 1].
 \]
 Implemented as \verb|Axis.Locate| in \verb|internal/numerics/grid.go|.
 
+\paragraph{Both ends are closed.} The accepted range is
+$[\ell, \ell + (N-1)s]$, and the upper end resolves to the last cell at
+$f = 1$ rather than opening a cell with no neighbour above it. This matters
+on the latitude axis, where $\ell + (N-1)s = 90^\circ$ is the north pole:
+that row carries real data --- NCEP resolves the GFS pole row per longitude
+--- so it is a usable grid row like any other. While the upper end was open,
+sampling exactly $90^\circ$ returned an error; the wind model discarded that
+error and returned a zero rate, so a prediction launched at the pole froze in
+place, and the wind-field endpoint reported a row of calm where a 29~m/s flow
+was blowing. The same argument applies to the last forecast hour and the
+topmost pressure level.
+
+Bound-checking is done on $p = (v - \ell)/s$ before truncation. Testing the
+truncated index instead admitted values just below $\ell$, because Go's
+\verb|int()| truncates toward zero: $p = -0{.}002$ became index $0$ and
+extrapolated off the end of the axis.
+
 \paragraph{Wrapping axes.} For periodic axes (e.g.\ longitude), the
 sequence is extended by the convention $x_N = x_0$ so a value approaching
 $x_N$ from below brackets $(N{-}1, 0)$ with fraction
@@ -194,7 +252,8 @@ $f = (v - x_{N-1})/s$.
 
 \paragraph{Worked example.} Latitude axis with $\ell = -90$, $s = 0{.}5$,
 $N = 361$. Query $v = -89{.}75$ yields $p = 0{.}5$, so $i_0 = 0$,
-$i_1 = 1$, $f = 0{.}5$.
+$i_1 = 1$, $f = 0{.}5$. Query $v = 90$ yields $p = 360$, which clamps to
+$i_0 = 359$, $i_1 = 360$, $f = 1$ --- the north pole row at full weight.
 
 \subsection{Multilinear interpolation}
 
@@ -238,8 +297,56 @@ and step $\Delta t$, \verb|RK4Step| applies
 \end{aligned}
 \]
 Reverse-time integration uses $\Delta t < 0$ unchanged; the implementation
-contains no branch on the sign of $\Delta t$. Domain-specific vector
-arithmetic (longitude wrap) is injected via \verb|VecAdd|.
+contains no branch on the sign of $\Delta t$.
+
+\paragraph{Stage combination on the sphere.} The additions above are not
+performed in $(\varphi, \lambda, h)$. Each stage rate is a velocity in the
+local horizontal frame at \emph{its own} evaluation point, and that frame
+rotates from one stage to the next --- near a pole, fast enough that
+averaging east/north components directly would reintroduce the error this
+formulation exists to remove. The stages are therefore converted to
+earth-centred vectors, combined with the usual $\tfrac16(1,2,2,1)$
+weights, read back in the frame at the starting point (which also discards
+the small radial component the averaging introduces), and applied as a
+single great-circle step:
+\[
+  \bar{\mathbf V} = \sum_i w_i \bigl(k_i^{E}\,\hat{\mathbf e}_i + k_i^{N}\,\hat{\mathbf n}_i\bigr),
+  \qquad
+  y(t + \Delta t) = \mathrm{GeoStep}\bigl(y,\; \bar{\mathbf V},\; \Delta t\bigr).
+\]
+
+\subsection{Great-circle stepping}
+\label{sec:geostep}
+
+\verb|GeoStep| advances a position by a horizontal velocity
+$(u, v)$ and a vertical rate $w$ over $\Delta t$. With
+$\hat{\mathbf r}, \hat{\mathbf e}, \hat{\mathbf n}$ the outward radial,
+east and north unit vectors at $(\varphi, \lambda)$, speed
+$V = \sqrt{u^2 + v^2}$ and unit travel direction
+$\hat{\mathbf t} = (u\,\hat{\mathbf e} + v\,\hat{\mathbf n})/V$:
+\[
+  \delta = \frac{V \Delta t}{R + h},
+  \qquad
+  \hat{\mathbf r}' = \hat{\mathbf r}\cos\delta + \hat{\mathbf t}\sin\delta,
+  \qquad
+  h' = h + w\,\Delta t,
+\]
+and $(\varphi', \lambda')$ are read back from $\hat{\mathbf r}'$ via
+$\varphi' = \arcsin r'_z$, $\lambda' = \operatorname{atan2}(r'_y, r'_x)$.
+
+The step is exact for a constant rate. Because $\hat{\mathbf t}$ is
+orthogonal to $\hat{\mathbf r}$ by construction, $\hat{\mathbf r}'$ stays
+on the unit sphere without renormalisation. All three basis vectors are
+unit length at every latitude including the poles; what happens at a pole
+is not a degeneracy but a genuine ambiguity, since east and north there
+depend on which meridian $\lambda$ names. That matches the data --- NCEP
+resolves the GFS pole row per longitude for exactly this reason, so any
+choice yields the same physical vector.
+
+A step that reaches a pole simply continues down the far side and
+$\lambda$ picks up $180^\circ$ on its own, with no special case and no
+threshold latitude. This is asserted directly in
+\verb|numerics/spherical_test.go|.
 
 \subsection{Termination refinement}
 
@@ -335,13 +442,15 @@ arbitrary State types via generics in numerics; the engine could lift
 its State to $(\mathbf{s}, \mathbf{v}_p)$ for a future mass-aware
 propagator without breaking the existing models.
 
-\paragraph{Coordinate system upgrades.} Migrating to WGS84/ECEF would
-remove the cosine factor in the horizontal wind transport equation and
-make distances metric directly. GFS itself uses a spherical Earth; the
-wind components are not directly portable. A clean implementation
-provides a coordinate-system parameter on the profile request; for now,
-the spherical model is used uniformly so that outputs remain bit
-identical to the upstream Tawhiri.
+\paragraph{Coordinate system upgrades.} The cosine factor is no longer a
+deferral: horizontal motion is integrated by great-circle rotation and the
+$1/\cos\varphi$ term is gone from the formulation entirely. What remains
+deferred is the \emph{ellipsoid}: migrating from a spherical Earth to
+WGS84 would make distances metric directly, but GFS itself uses a
+spherical Earth and its wind components are not directly portable to the
+ellipsoid. A clean implementation would provide a coordinate-system
+parameter on the profile request; for now the spherical model is used
+uniformly.
 
 \paragraph{Monte Carlo.} GEFS already provides 21 ensemble members per
 epoch. A Monte Carlo prediction would sample $K$ trajectories per
diff --git a/internal/engine/engine_test.go b/internal/engine/engine_test.go
index 8c3fd57..ccf0f88 100644
--- a/internal/engine/engine_test.go
+++ b/internal/engine/engine_test.go
@@ -119,14 +119,14 @@ func TestPiecewiseRate(t *testing.T) {
 		{Until: math.Inf(1), Rate: 0},
 	})
 
-	if r := m(50, State{}); r.Altitude != 5 {
-		t.Errorf("rate at t=50 = %v, want 5", r.Altitude)
+	if r := m(50, State{}); r.Vertical != 5 {
+		t.Errorf("rate at t=50 = %v, want 5", r.Vertical)
 	}
-	if r := m(150, State{}); r.Altitude != 3 {
-		t.Errorf("rate at t=150 = %v, want 3", r.Altitude)
+	if r := m(150, State{}); r.Vertical != 3 {
+		t.Errorf("rate at t=150 = %v, want 3", r.Vertical)
 	}
-	if r := m(300, State{}); r.Altitude != 0 {
-		t.Errorf("rate at t=300 = %v, want 0", r.Altitude)
+	if r := m(300, State{}); r.Vertical != 0 {
+		t.Errorf("rate at t=300 = %v, want 0", r.Vertical)
 	}
 }
 
@@ -149,11 +149,11 @@ func TestPiecewiseReferenceResolution(t *testing.T) {
 	ctx := StageContext{ProfileStart: 1000, PropagatorStart: 5000}
 	m := built.Build(ctx)
 	// Until=100 from propagator_start=5000 → absolute 5100.
-	if r := m(5050, State{}); r.Altitude != 5 {
-		t.Errorf("rate at t=5050 = %v, want 5", r.Altitude)
+	if r := m(5050, State{}); r.Vertical != 5 {
+		t.Errorf("rate at t=5050 = %v, want 5", r.Vertical)
 	}
-	if r := m(5150, State{}); r.Altitude != 3 {
-		t.Errorf("rate at t=5150 = %v, want 3", r.Altitude)
+	if r := m(5150, State{}); r.Vertical != 3 {
+		t.Errorf("rate at t=5150 = %v, want 3", r.Vertical)
 	}
 }
 
@@ -166,22 +166,20 @@ func (w fixedWind) Wind(_ float64, _, _, _ float64) (weather.Sample, error) {
 func (fixedWind) Epoch() time.Time { return time.Unix(0, 0) }
 func (fixedWind) Source() string   { return "test-fixed" }
 
-func TestWindTransportUnitConversion(t *testing.T) {
-	wind := WindTransport(fixedWind{u: 10, v: 0}, nil)
-	d := wind(0, State{Lat: 0, Lng: 0, Altitude: 0})
-	wantLng := (180.0 / math.Pi) * 10.0 / 6371009.0
-	if math.Abs(d.Lng-wantLng) > 1e-12 {
-		t.Errorf("dlng = %v, want %v", d.Lng, wantLng)
-	}
-	if math.Abs(d.Lat) > 1e-12 {
-		t.Errorf("dlat = %v, want 0 for u=10 v=0", d.Lat)
-	}
+func TestWindTransportPassesWindThroughUnchanged(t *testing.T) {
+	// The wind field already gives a horizontal velocity, so the propagator
+	// receives it verbatim. This replaces a test of the old deg/s conversion,
+	// whose 1/cos(lat) factor is exactly what made the poles unusable.
+	wind := WindTransport(fixedWind{u: 10, v: -4}, nil)
 
-	wind2 := WindTransport(fixedWind{u: 0, v: 5}, nil)
-	d = wind2(0, State{Lat: 60, Lng: 0, Altitude: 0})
-	wantLat := (180.0 / math.Pi) * 5.0 / 6371009.0
-	if math.Abs(d.Lat-wantLat) > 1e-12 {
-		t.Errorf("dlat at lat=60 = %v, want %v", d.Lat, wantLat)
+	for _, lat := range []float64{0, 45, 60, 89, 89.999} {
+		r := wind(0, State{Lat: lat, Lng: 0, Altitude: 0})
+		if r.East != 10 || r.North != -4 {
+			t.Errorf("lat %g: rate = %+v, want East=10 North=-4", lat, r)
+		}
+		if r.Vertical != 0 {
+			t.Errorf("lat %g: wind must not produce vertical motion, got %v", lat, r.Vertical)
+		}
 	}
 }
 
diff --git a/internal/engine/models.go b/internal/engine/models.go
index 9a738d8..7399961 100644
--- a/internal/engine/models.go
+++ b/internal/engine/models.go
@@ -15,10 +15,10 @@ func Sum(models ...Model) Model {
 	if len(models) == 1 {
 		return models[0]
 	}
-	return func(t float64, s State) State {
-		var sum State
+	return func(t float64, s State) numerics.Rate {
+		var sum numerics.Rate
 		for _, m := range models {
-			sum = numerics.AddGeo(sum, m(t, s))
+			sum = numerics.AddRate(sum, m(t, s))
 		}
 		return sum
 	}
@@ -27,7 +27,7 @@ func Sum(models ...Model) Model {
 // ConstantRate returns a model with a constant vertical velocity (m/s).
 // Positive rates are upward.
 func ConstantRate(rate float64) Model {
-	return func(_ float64, _ State) State { return State{Altitude: rate} }
+	return func(_ float64, _ State) numerics.Rate { return numerics.Rate{Vertical: rate} }
 }
 
 // ParachuteDescent returns a model where vertical velocity grows with
@@ -40,8 +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 {
-	return func(_ float64, s State) State {
-		return State{Altitude: numerics.DragTerminalVelocity(seaLevelRate, s.Altitude)}
+	return func(_ float64, s State) numerics.Rate {
+		return numerics.Rate{Vertical: numerics.DragTerminalVelocity(seaLevelRate, s.Altitude)}
 	}
 }
 
@@ -64,33 +64,34 @@ func Piecewise(segments []RateSegment) Model {
 	sort.Slice(sorted, func(i, j int) bool { return sorted[i].Until < sorted[j].Until })
 	finalRate := sorted[len(sorted)-1].Rate
 
-	return func(t float64, _ State) State {
+	return func(t float64, _ State) numerics.Rate {
 		idx := sort.Search(len(sorted), func(i int) bool { return sorted[i].Until > t })
 		if idx == len(sorted) {
-			return State{Altitude: finalRate}
+			return numerics.Rate{Vertical: finalRate}
 		}
-		return State{Altitude: sorted[idx].Rate}
+		return numerics.Rate{Vertical: sorted[idx].Rate}
 	}
 }
 
 // WindTransport returns a model that moves laterally at the wind velocity
 // 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.
+// wind is handed straight to the integrator as a horizontal velocity.
 //
 // 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 {
-	return func(t float64, s State) State {
+	return func(t float64, s State) numerics.Rate {
 		sample, err := field.Wind(t, s.Lat, s.Lng, s.Altitude)
 		if err != nil {
-			return State{}
+			return numerics.Rate{}
 		}
 		if sample.AboveModel && events != nil {
 			events.Emit("above_model", t, s,
 				"altitude exceeded the highest pressure level of the wind dataset; samples extrapolated")
 		}
-		dLat, dLng := numerics.WindToGeoRate(sample.U, sample.V, s.Lat, s.Altitude)
-		return State{Lat: dLat, Lng: dLng}
+		// The wind is already a horizontal velocity; it is handed over as-is.
+		// Converting it to deg/s here is what used to blow up near the poles.
+		return numerics.Rate{East: sample.U, North: sample.V}
 	}
 }
diff --git a/internal/engine/propagator.go b/internal/engine/propagator.go
index 19b080e..3763ec7 100644
--- a/internal/engine/propagator.go
+++ b/internal/engine/propagator.go
@@ -71,7 +71,7 @@ func (p *Propagator) run(ctx StageContext, t0 float64, s0 State, globals []Const
 		constraints = p.BuildConstraints(ctx)
 	}
 
-	field := numerics.Field(model)
+	field := numerics.RateField(model)
 
 	out := Result{Propagator: p.Name, Outcome: OutcomeContinued, Path: numerics.NewPath(estimatedSteps)}
 	out.Path.Append(t0, s0)
diff --git a/internal/engine/types.go b/internal/engine/types.go
index 6050c97..91de9b8 100644
--- a/internal/engine/types.go
+++ b/internal/engine/types.go
@@ -20,11 +20,15 @@ import "predictor-refactored/internal/numerics"
 // the numeric core share one hot-path value type without conversions.
 type State = numerics.GeoVec
 
-// Model returns the time derivative of state at (t, s).
+// Model returns the rate of change of state at (t, s), as a velocity in the
+// local horizontal frame (metres per second east/north/up).
 //
-// The derivative is direction-independent; the integrator applies the
-// sign of dt for reverse propagation.
-type Model func(t float64, s State) State
+// It is deliberately not a lat/lon derivative: that form carries a 1/cos(lat)
+// factor in longitude which diverges at the poles. See numerics.Rate.
+//
+// The rate is direction-independent; the integrator applies the sign of dt for
+// reverse propagation.
+type Model func(t float64, s State) numerics.Rate
 
 // Direction is the time direction of integration.
 type Direction int8
diff --git a/internal/numerics/doc.go b/internal/numerics/doc.go
index 807ba3d..8deaf1b 100644
--- a/internal/numerics/doc.go
+++ b/internal/numerics/doc.go
@@ -1,11 +1,19 @@
 // Package numerics provides the numerical primitives used by the trajectory
-// engine: regular-grid multilinear interpolation, monotone bisection, and
-// a generic explicit Runge-Kutta-4 integrator with binary-search refinement
-// of a termination point.
+// engine: regular-grid multilinear interpolation, monotone bisection, spherical
+// kinematics, and a Runge-Kutta-4 integrator with binary-search refinement of a
+// termination point.
 //
-// The package has no dependencies on any domain type. State and derivative
-// types are generic, and all coordinate-wrap or unit-conversion semantics
-// live in the caller.
+// Positions are carried as GeoVec (degrees and metres) and advanced with
+// GeoStep, which rotates the position along a great circle. Rates are velocities
+// (Rate: metres per second east/north/up), never degrees per second — a
+// longitude derivative carries a 1/cos(lat) factor that diverges at the poles,
+// so the singularity is absent from the formulation rather than guarded against
+// by a threshold latitude. A step that reaches a pole continues down the far
+// side on its own.
+//
+// The package has no dependencies on any domain type. Longitude wrapping into
+// [0, 360) and the sphere geometry live here; unit conversions specific to a
+// data source stay in the caller.
 //
 // All algorithms are documented in docs/numerics.tex.
 package numerics
diff --git a/internal/numerics/grid.go b/internal/numerics/grid.go
index 720618d..b9e019c 100644
--- a/internal/numerics/grid.go
+++ b/internal/numerics/grid.go
@@ -1,6 +1,9 @@
 package numerics
 
-import "fmt"
+import (
+	"fmt"
+	"math"
+)
 
 // Axis describes a regularly-spaced grid axis with N grid points,
 // values left, left+step, left+2*step, ..., left+(N-1)*step.
@@ -28,27 +31,43 @@ func (e *AxisError) Error() string {
 
 // Bracket holds the two surrounding grid indices and the fractional position
 // of a value within an axis. The weight at Lo is (1 - Frac); the weight at Hi
-// is Frac. Frac lies in [0, 1).
+// is Frac. Frac lies in [0, 1].
 type Bracket struct {
 	Lo, Hi int
 	Frac   float64
 }
 
 // Locate returns the bracket containing value within the axis.
-// For a non-wrapping axis, value must lie in [Left, Left + (N-1)*Step);
-// for a wrapping axis, value must lie in [Left, Left + N*Step).
+// The accepted range is closed at both ends: [Left, Left + (N-1)*Step] for a
+// non-wrapping axis, [Left, Left + N*Step] for a wrapping one.
+//
+// The upper end is closed deliberately. On the GFS latitude axis it is the
+// north pole, whose row holds real data — NCEP resolves it per longitude, so it
+// is a usable grid row like any other. Rejecting it used to abort the wind
+// lookup, and because the caller discarded that error the whole prediction
+// silently froze at latitude 90. The same argument applies to the last forecast
+// hour and the topmost pressure level.
 func (a Axis) Locate(value float64) (Bracket, error) {
 	pos := (value - a.Left) / a.Step
-	lo := int(pos) // truncates toward zero; pos is non-negative for valid inputs
 
 	maxLo := a.N - 2
 	if a.Wrap {
 		maxLo = a.N - 1
 	}
-	if lo < 0 || lo > maxLo {
+	// Bound-check in float space. Checking the truncated index instead would
+	// let values just below Left through: int() truncates toward zero, so a pos
+	// of -0.002 became index 0 rather than -1 and extrapolated off the end.
+	if pos < 0 || pos > float64(maxLo+1) {
 		return Bracket{}, &AxisError{Axis: a.Name, Value: value}
 	}
 
+	lo := int(math.Floor(pos))
+	// The exact upper bound belongs to the top cell at Frac 1, rather than
+	// opening a cell that has no neighbour above it.
+	if lo > maxLo {
+		lo = maxLo
+	}
+
 	hi := lo + 1
 	if a.Wrap && hi == a.N {
 		hi = 0
diff --git a/internal/numerics/grid_test.go b/internal/numerics/grid_test.go
index 342d39c..62820f4 100644
--- a/internal/numerics/grid_test.go
+++ b/internal/numerics/grid_test.go
@@ -23,9 +23,11 @@ func TestAxisLocate(t *testing.T) {
 		t.Errorf("Locate(-89.75) = %+v, %v; want frac=0.5", b, err)
 	}
 
-	// 90 is exactly on the upper boundary — there's no Hi above it
-	if _, err := a.Locate(90); err == nil {
-		t.Errorf("Locate(90) should error, got nil")
+	// 90 is exactly on the upper boundary. It is now accepted as the far edge of
+	// the last cell: on the GFS latitude axis that is the north pole, whose row
+	// carries real data. Rejecting it used to freeze predictions there.
+	if b, err := a.Locate(90); err != nil || b.Lo != 359 || b.Hi != 360 || math.Abs(b.Frac-1) > 1e-12 {
+		t.Errorf("Locate(90) = %+v, %v; want {359 360 1}", b, err)
 	}
 
 	if _, err := a.Locate(-91); err == nil {
@@ -47,9 +49,16 @@ func TestAxisLocateWrap(t *testing.T) {
 		t.Errorf("Locate(359.75) = %+v, %v; want {719 0 0.5}", b, err)
 	}
 
-	// 360 is outside the half-open interval
-	if _, err := a.Locate(360); err == nil {
-		t.Errorf("Locate(360) should error, got nil")
+	// 360 is the wrap point and now resolves to it: the far edge of the last
+	// cell, whose Hi is index 0. Weight 1 there means exactly 0 degrees, which
+	// is what 360 means. Callers normalise longitude anyway, so this is a
+	// consistency property rather than a path anyone relies on.
+	if b, err := a.Locate(360); err != nil || b.Lo != 719 || b.Hi != 0 || math.Abs(b.Frac-1) > 1e-12 {
+		t.Errorf("Locate(360) = %+v, %v; want {719 0 1}", b, err)
+	}
+
+	if _, err := a.Locate(360.5); err == nil {
+		t.Errorf("Locate(360.5) should error, got nil")
 	}
 }
 
@@ -92,3 +101,56 @@ func TestLerp(t *testing.T) {
 		t.Errorf("Lerp(10, 20, 0.25) != 12.5")
 	}
 }
+
+// The GFS latitude axis runs -90..90 at 0.5 deg (N=361), so the north pole is
+// the axis's exact upper bound. Bracketing must accept it: the pole row holds
+// real data (NCEP resolves it per longitude via POLFIXV), and refusing it made
+// the whole prediction freeze silently at latitude 90. The same applies to the
+// last forecast hour and the topmost pressure level.
+func TestAxisLocateAcceptsExactUpperBound(t *testing.T) {
+	t.Parallel()
+
+	lat := Axis{Left: -90, Step: 0.5, N: 361, Name: "lat"}
+
+	tests := []struct {
+		name     string
+		value    float64
+		wantLo   int
+		wantHi   int
+		wantFrac float64
+	}{
+		{name: "exact lower bound", value: -90, wantLo: 0, wantHi: 1, wantFrac: 0},
+		{name: "interior point", value: 0.25, wantLo: 180, wantHi: 181, wantFrac: 0.5},
+		{name: "one cell below the top", value: 89.5, wantLo: 359, wantHi: 360, wantFrac: 0},
+		{name: "inside the top cell", value: 89.75, wantLo: 359, wantHi: 360, wantFrac: 0.5},
+		// The case that used to error: the far edge of the last cell.
+		{name: "exact upper bound is the top of the last cell", value: 90, wantLo: 359, wantHi: 360, wantFrac: 1},
+	}
+
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			t.Parallel()
+			b, err := lat.Locate(tt.value)
+			if err != nil {
+				t.Fatalf("Locate(%v) returned error: %v", tt.value, err)
+			}
+			if b.Lo != tt.wantLo || b.Hi != tt.wantHi {
+				t.Errorf("Locate(%v) = lo %d hi %d, want lo %d hi %d", tt.value, b.Lo, b.Hi, tt.wantLo, tt.wantHi)
+			}
+			if math.Abs(b.Frac-tt.wantFrac) > 1e-12 {
+				t.Errorf("Locate(%v) frac = %v, want %v", tt.value, b.Frac, tt.wantFrac)
+			}
+		})
+	}
+}
+
+func TestAxisLocateStillRejectsOutOfRange(t *testing.T) {
+	t.Parallel()
+
+	lat := Axis{Left: -90, Step: 0.5, N: 361, Name: "lat"}
+	for _, v := range []float64{-90.001, 90.001, 91, -100} {
+		if _, err := lat.Locate(v); err == nil {
+			t.Errorf("Locate(%v) accepted a value outside the axis", v)
+		}
+	}
+}
diff --git a/internal/numerics/motion_test.go b/internal/numerics/motion_test.go
deleted file mode 100644
index 3bc8051..0000000
--- a/internal/numerics/motion_test.go
+++ /dev/null
@@ -1,58 +0,0 @@
-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/ode.go b/internal/numerics/ode.go
index 9a199ce..504b67e 100644
--- a/internal/numerics/ode.go
+++ b/internal/numerics/ode.go
@@ -1,31 +1,8 @@
 package numerics
 
-// Field returns the time derivative of a geographic state at (t, y).
-// The derivative is direction-independent; the integrator applies the sign
-// of dt for reverse-time integration.
-type Field func(t float64, y GeoVec) GeoVec
-
 // Crossed reports whether a termination condition holds at (t, y).
 type Crossed func(t float64, y GeoVec) bool
 
-// RK4Step performs one classical Runge-Kutta-4 step from (t, y) with step dt.
-// dt may be negative to integrate backwards in time. Longitude wrapping is
-// applied at every intermediate add via GeoAdd, matching the reference
-// integrator. The function performs no heap allocation.
-func RK4Step(t float64, y GeoVec, dt float64, f Field) GeoVec {
-	half := dt / 2
-	k1 := f(t, y)
-	k2 := f(t+half, GeoAdd(y, half, k1))
-	k3 := f(t+half, GeoAdd(y, half, k2))
-	k4 := f(t+dt, GeoAdd(y, dt, k3))
-
-	y2 := GeoAdd(y, dt/6, k1)
-	y2 = GeoAdd(y2, dt/3, k2)
-	y2 = GeoAdd(y2, dt/3, k3)
-	y2 = GeoAdd(y2, dt/6, k4)
-	return y2
-}
-
 // RefineCrossing locates a crossing between (t1, y1) (not crossed) and
 // (t2, y2) (crossed) by binary search in the linear-interpolation parameter
 // space, stopping when the parameter interval is narrower than tol.
diff --git a/internal/numerics/ode_test.go b/internal/numerics/ode_test.go
index 0a0697b..9aa9933 100644
--- a/internal/numerics/ode_test.go
+++ b/internal/numerics/ode_test.go
@@ -7,7 +7,7 @@ import (
 
 func TestRK4ExponentialDecay(t *testing.T) {
 	// dAlt/dt = -Alt → exact: Alt(t) = Alt0 * exp(-t).
-	f := func(_ float64, y GeoVec) GeoVec { return GeoVec{Altitude: -y.Altitude} }
+	f := func(_ float64, y GeoVec) Rate { return Rate{Vertical: -y.Altitude} }
 
 	y := GeoVec{Altitude: 1}
 	tnow, dt := 0.0, 0.01
@@ -23,7 +23,7 @@ func TestRK4ExponentialDecay(t *testing.T) {
 
 func TestRK4ReverseTime(t *testing.T) {
 	// dAlt/dt = Alt → exact: Alt(t) = Alt0 * exp(t).
-	f := func(_ float64, y GeoVec) GeoVec { return GeoVec{Altitude: y.Altitude} }
+	f := func(_ float64, y GeoVec) Rate { return Rate{Vertical: y.Altitude} }
 
 	y := GeoVec{Altitude: math.E}
 	tnow, dt := 1.0, -0.01
@@ -50,13 +50,6 @@ func TestRefineCrossing(t *testing.T) {
 	}
 }
 
-func TestGeoAddWrapsLongitude(t *testing.T) {
-	y := GeoAdd(GeoVec{Lng: 350}, 1, GeoVec{Lng: 20})
-	if math.Abs(y.Lng-10) > 1e-9 {
-		t.Errorf("GeoAdd wrap: lng = %v, want 10", y.Lng)
-	}
-}
-
 func TestGeoLerpWrap(t *testing.T) {
 	mid := GeoLerp(GeoVec{Lng: 350}, GeoVec{Lng: 10}, 0.5)
 	if math.Abs(mid.Lng) > 1e-9 && math.Abs(mid.Lng-360) > 1e-9 {
diff --git a/internal/numerics/spherical.go b/internal/numerics/spherical.go
new file mode 100644
index 0000000..c7c5dd9
--- /dev/null
+++ b/internal/numerics/spherical.go
@@ -0,0 +1,159 @@
+package numerics
+
+import "math"
+
+// Spherical kinematics for trajectory integration.
+//
+// Positions are carried as GeoVec (degrees) because constraints, path recording
+// and the API all speak latitude and longitude. Motion, however, is expressed
+// as a velocity in metres per second and applied by rotating the position
+// vector along a great circle. Nothing in this file forms a longitude
+// derivative, which is what removes the polar singularity: dLng/dt carries a
+// 1/cos(lat) factor that reached 1.46e12 deg/s at 90 degrees and made any
+// fixed-step integrator meaningless there.
+//
+// Consequences worth knowing:
+//   - A step that reaches a pole continues down the far side, and longitude
+//     picks up 180 degrees on its own. No special case, no threshold latitude.
+//   - Latitude cannot leave [-90, 90] by construction, because it is read back
+//     from a unit vector rather than accumulated.
+
+const (
+	degToRad = math.Pi / 180
+	radToDeg = 180 / math.Pi
+)
+
+// Rate is a velocity in the local horizontal frame at a point: metres per
+// second toward east and north, plus metres per second upward.
+type Rate struct {
+	East     float64
+	North    float64
+	Vertical float64
+}
+
+// AddRate sums two rates componentwise. Rates compose linearly; positions do
+// not, which is why they are advanced with GeoStep instead.
+func AddRate(a, b Rate) Rate {
+	return Rate{East: a.East + b.East, North: a.North + b.North, Vertical: a.Vertical + b.Vertical}
+}
+
+// basis returns the earth-centred unit vectors at (lat, lng): outward radial,
+// east and north.
+//
+// All three are unit length at every latitude, the poles included. What happens
+// at a pole is not a degeneracy but a genuine ambiguity: east and north there
+// depend on which meridian the longitude names. That matches the data — NCEP
+// resolves the GFS pole row per longitude for exactly this reason, so any
+// choice yields the same physical vector.
+func basis(latDeg, lngDeg float64) (radial, east, north [3]float64) {
+	sinLat, cosLat := math.Sincos(latDeg * degToRad)
+	sinLng, cosLng := math.Sincos(lngDeg * degToRad)
+	radial = [3]float64{cosLat * cosLng, cosLat * sinLng, sinLat}
+	east = [3]float64{-sinLng, cosLng, 0}
+	north = [3]float64{-sinLat * cosLng, -sinLat * sinLng, cosLat}
+	return radial, east, north
+}
+
+// toGeo reads a position back off a unit vector.
+func toGeo(v [3]float64, altitude float64) GeoVec {
+	return GeoVec{
+		Lat:      math.Asin(math.Max(-1, math.Min(1, v[2]))) * radToDeg,
+		Lng:      PyMod(math.Atan2(v[1], v[0])*radToDeg, 360),
+		Altitude: altitude,
+	}
+}
+
+// GeoStep advances a position by rate over dt seconds along a great circle.
+// dt may be negative, which travels the same arc in the opposite direction.
+//
+// It is exact for a constant rate: the position vector is rotated in the plane
+// it spans with the direction of travel. Since that direction is orthogonal to
+// the radial by construction, the result stays on the unit sphere without
+// renormalisation.
+func GeoStep(y GeoVec, rate Rate, dt float64) GeoVec {
+	altitude := y.Altitude + rate.Vertical*dt
+	speed := math.Hypot(rate.East, rate.North)
+	if speed == 0 {
+		return GeoVec{Lat: y.Lat, Lng: y.Lng, Altitude: altitude}
+	}
+
+	radial, east, north := basis(y.Lat, y.Lng)
+	var tangent [3]float64
+	for i := range tangent {
+		tangent[i] = (rate.East*east[i] + rate.North*north[i]) / speed
+	}
+
+	// Angle subtended at the earth's centre by the arc travelled.
+	angle := speed * dt / (EarthRadius + y.Altitude)
+	sinA, cosA := math.Sincos(angle)
+	var out [3]float64
+	for i := range out {
+		out[i] = radial[i]*cosA + tangent[i]*sinA
+	}
+	return toGeo(out, altitude)
+}
+
+// GreatCircleMetres is the surface distance between two positions, ignoring
+// altitude.
+//
+// Uses the chord rather than acos(dot): for nearby points acos loses most of
+// its significant digits, and these distances are checked to sub-millimetre
+// tolerances in tests.
+func GreatCircleMetres(a, b GeoVec) float64 {
+	ra, _, _ := basis(a.Lat, a.Lng)
+	rb, _, _ := basis(b.Lat, b.Lng)
+	var chordSq float64
+	for i := range ra {
+		d := ra[i] - rb[i]
+		chordSq += d * d
+	}
+	return 2 * EarthRadius * math.Asin(math.Min(1, math.Sqrt(chordSq)/2))
+}
+
+// RateField returns the rate of change of state at (t, y). The rate is
+// direction-independent; the integrator applies the sign of dt for reverse-time
+// integration.
+type RateField func(t float64, y GeoVec) Rate
+
+// RK4Step performs one classical Runge-Kutta-4 step along the sphere.
+//
+// The four stage rates are combined as earth-centred vectors rather than as
+// local east/north pairs. That distinction matters: the local frame rotates
+// between stages, and near a pole it rotates fast enough that averaging
+// components directly would reintroduce the very error this formulation exists
+// to remove. The combined velocity is then read back in the frame at the
+// starting point — which also discards the small radial component averaging
+// introduces — and applied as a single great-circle step.
+func RK4Step(t float64, y GeoVec, dt float64, f RateField) GeoVec {
+	half := dt / 2
+
+	k1 := f(t, y)
+	y2 := GeoStep(y, k1, half)
+	k2 := f(t+half, y2)
+	y3 := GeoStep(y, k2, half)
+	k3 := f(t+half, y3)
+	y4 := GeoStep(y, k3, dt)
+	k4 := f(t+dt, y4)
+
+	points := [4]GeoVec{y, y2, y3, y4}
+	rates := [4]Rate{k1, k2, k3, k4}
+	weights := [4]float64{1.0 / 6, 1.0 / 3, 1.0 / 3, 1.0 / 6}
+
+	var vx, vy, vz, vertical float64
+	for i := range points {
+		_, east, north := basis(points[i].Lat, points[i].Lng)
+		w := weights[i]
+		vx += w * (rates[i].East*east[0] + rates[i].North*north[0])
+		vy += w * (rates[i].East*east[1] + rates[i].North*north[1])
+		vz += w * (rates[i].East*east[2] + rates[i].North*north[2])
+		vertical += w * rates[i].Vertical
+	}
+
+	_, east0, north0 := basis(y.Lat, y.Lng)
+	mean := Rate{
+		East:     vx*east0[0] + vy*east0[1] + vz*east0[2],
+		North:    vx*north0[0] + vy*north0[1] + vz*north0[2],
+		Vertical: vertical,
+	}
+	return GeoStep(y, mean, dt)
+}
diff --git a/internal/numerics/spherical_test.go b/internal/numerics/spherical_test.go
new file mode 100644
index 0000000..5dc0c61
--- /dev/null
+++ b/internal/numerics/spherical_test.go
@@ -0,0 +1,230 @@
+package numerics
+
+import (
+	"math"
+	"testing"
+)
+
+// Closed-form checks for great-circle stepping. Each expectation is an exact
+// analytic result, not a golden value copied from a previous run.
+
+const tolDeg = 1e-9
+
+func TestGeoStep(t *testing.T) {
+	t.Parallel()
+
+	// Angular distance covered by `speed` for `dt` at sea level, in degrees.
+	arcDeg := func(speed, dt float64) float64 {
+		return speed * dt / EarthRadius * 180 / math.Pi
+	}
+
+	tests := []struct {
+		name             string
+		start            GeoVec
+		rate             Rate
+		dt               float64
+		wantLat, wantLng float64
+		wantAlt          float64
+	}{
+		{
+			name:    "eastward at the equator advances longitude only",
+			start:   GeoVec{Lat: 0, Lng: 0},
+			rate:    Rate{East: 10},
+			dt:      100,
+			wantLat: 0,
+			wantLng: arcDeg(10, 100),
+		},
+		{
+			name:    "northward at the equator advances latitude only",
+			start:   GeoVec{Lat: 0, Lng: 0},
+			rate:    Rate{North: 10},
+			dt:      100,
+			wantLat: arcDeg(10, 100),
+			wantLng: 0,
+		},
+		{
+			name:    "zero horizontal rate leaves the position alone",
+			start:   GeoVec{Lat: 51.5, Lng: 359.9, Altitude: 1000},
+			rate:    Rate{Vertical: 5},
+			dt:      10,
+			wantLat: 51.5,
+			wantLng: 359.9,
+			wantAlt: 1050,
+		},
+		{
+			// The whole point of the change: 111 m from the pole, a due-north
+			// step must pass over the pole and come down the far meridian.
+			// Start 0.001 deg from the pole, travel 600 m (0.005396 deg), so it
+			// overshoots by 0.004396 deg on longitude 30+180.
+			name:    "due north over the pole flips longitude by 180",
+			start:   GeoVec{Lat: 89.999, Lng: 30},
+			rate:    Rate{North: 10},
+			dt:      60,
+			wantLat: 90 - (arcDeg(10, 60) - 0.001),
+			wantLng: 210,
+		},
+		{
+			name:    "due south over the south pole flips longitude by 180",
+			start:   GeoVec{Lat: -89.999, Lng: 200},
+			rate:    Rate{North: -10},
+			dt:      60,
+			wantLat: -90 + (arcDeg(10, 60) - 0.001),
+			wantLng: 20,
+		},
+		{
+			name:    "longitude stays wrapped into [0,360)",
+			start:   GeoVec{Lat: 0, Lng: 359.999},
+			rate:    Rate{East: 100},
+			dt:      100,
+			wantLat: 0,
+			wantLng: math.Mod(359.999+arcDeg(100, 100), 360),
+		},
+	}
+
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			t.Parallel()
+			got := GeoStep(tt.start, tt.rate, tt.dt)
+			if math.Abs(got.Lat-tt.wantLat) > tolDeg {
+				t.Errorf("Lat = %.12f, want %.12f", got.Lat, tt.wantLat)
+			}
+			if math.Abs(got.Lng-tt.wantLng) > tolDeg {
+				t.Errorf("Lng = %.12f, want %.12f", got.Lng, tt.wantLng)
+			}
+			if math.Abs(got.Altitude-tt.wantAlt) > 1e-9 {
+				t.Errorf("Altitude = %v, want %v", got.Altitude, tt.wantAlt)
+			}
+		})
+	}
+}
+
+// The defect this replaces: dLng went as 1/cos(lat), reaching 1.46e12 deg/s at
+// the pole. Ground displacement must instead stay equal to speed*dt at every
+// latitude, because that is what physically happens.
+func TestGeoStepGroundDistanceIndependentOfLatitude(t *testing.T) {
+	t.Parallel()
+
+	const speed, dt = 10.0, 60.0
+	want := speed * dt // 600 m
+
+	for _, lat := range []float64{0, 45, 60, 85.051129, 89, 89.9, 89.99, 89.999, 90} {
+		start := GeoVec{Lat: lat, Lng: 17}
+		got := GeoStep(start, Rate{East: speed}, dt)
+		d := GreatCircleMetres(start, got)
+		// 0.1 mm. The point is to catch a divergence — the old code was wrong by
+		// a factor of 1e12 here — not to police the last bit: near the pole
+		// cos(lat) is ~1e-5, so double precision limits this to a few microns.
+		if math.Abs(d-want) > 1e-4 {
+			t.Errorf("lat %g: ground distance = %.9f m, want %.9f m", lat, d, want)
+		}
+	}
+}
+
+func TestGeoStepReverseGoesBackAlongTheSameCircle(t *testing.T) {
+	t.Parallel()
+
+	// A negative dt must travel the same arc in the opposite direction. Note it
+	// is NOT true that GeoStep(GeoStep(y, r, dt), r, -dt) == y: the local frame
+	// rotates during the step, so the same east/north pair means a different
+	// physical direction at the arrival point. The invariant that does hold is
+	// that the two endpoints straddle the start on one great circle.
+	const speed, dt = 12.0, 60.0
+	for _, lat := range []float64{0, 60, 89.99} {
+		y := GeoVec{Lat: lat, Lng: 100, Altitude: 5000}
+		r := Rate{East: speed, Vertical: 3}
+		fwd := GeoStep(y, r, dt)
+		back := GeoStep(y, r, -dt)
+
+		// The arc is flown at altitude, so its projection onto the surface —
+		// which is what GreatCircleMetres reports — is shorter by R/(R+alt).
+		wantGround := speed * dt * EarthRadius / (EarthRadius + y.Altitude)
+		if d := GreatCircleMetres(y, back); math.Abs(d-wantGround) > 1e-4 {
+			t.Errorf("lat %g: reverse arc = %.9f m, want %.9f m", lat, d, wantGround)
+		}
+		if d := GreatCircleMetres(fwd, back); math.Abs(d-2*wantGround) > 1e-4 {
+			t.Errorf("lat %g: forward/reverse separation = %.9f m, want %.9f m",
+				lat, d, 2*wantGround)
+		}
+		if math.Abs(back.Altitude-(y.Altitude-3*dt)) > 1e-9 {
+			t.Errorf("lat %g: reverse altitude = %v, want %v", lat, back.Altitude, y.Altitude-3*dt)
+		}
+	}
+}
+
+// Great-circle motion is rotation about a fixed axis, so a field returning the
+// local east/north components of `omega x r` has an exact solution: a single
+// rotation. RK4 must reproduce it, including next to the pole where the local
+// frame spins fastest between stages.
+//
+// A field with *constant* east/north would be a rhumb line, not a great circle,
+// so it cannot be used for this comparison.
+func TestRK4StepMatchesGreatCircle(t *testing.T) {
+	t.Parallel()
+
+	for _, lat := range []float64{0, 60, 89.9, 89.999} {
+		start := GeoVec{Lat: lat, Lng: 40, Altitude: 20000}
+		const speed, dt, n = 35.0, 60.0, 10
+
+		// Rate at the start point defines the great circle; GeoStep over the
+		// whole interval is then the exact answer.
+		initial := Rate{East: speed * 0.8, North: speed * 0.6}
+		exact := GeoStep(start, initial, dt*n)
+
+		// Pointwise field for that same great circle: rotation about the axis
+		// r0 x t0 at constant angular rate.
+		field := greatCircleField(start, initial)
+
+		stepped := start
+		for range n {
+			stepped = RK4Step(0, stepped, dt, field)
+		}
+
+		if d := GreatCircleMetres(stepped, exact); d > 0.01 {
+			t.Errorf("lat %g: RK4 drifted %.6f m from the exact great circle", lat, d)
+		}
+	}
+}
+
+func TestAddRate(t *testing.T) {
+	t.Parallel()
+
+	got := AddRate(Rate{East: 1, North: 2, Vertical: 3}, Rate{East: 10, North: 20, Vertical: 30})
+	want := Rate{East: 11, North: 22, Vertical: 33}
+	if got != want {
+		t.Errorf("AddRate = %+v, want %+v", got, want)
+	}
+}
+
+// greatCircleField builds a rate field whose exact solution is the great circle
+// through `start` with initial rate `initial`: rotation about the fixed axis
+// r0 x t0. At any point it returns the local east/north components of that
+// rotation's velocity, so the field is defined pointwise without needing to
+// know how far along the arc we are.
+func greatCircleField(start GeoVec, initial Rate) RateField {
+	speed := math.Hypot(initial.East, initial.North)
+	r0, e0, n0 := basis(start.Lat, start.Lng)
+	var t0 [3]float64
+	for i := range t0 {
+		t0[i] = (initial.East*e0[i] + initial.North*n0[i]) / speed
+	}
+	axis := cross(r0, t0)
+
+	return func(_ float64, y GeoVec) Rate {
+		r, e, n := basis(y.Lat, y.Lng)
+		v := cross(axis, r)
+		return Rate{
+			East:  speed * dot(v, e),
+			North: speed * dot(v, n),
+		}
+	}
+}
+
+func cross(a, b [3]float64) [3]float64 {
+	return [3]float64{
+		a[1]*b[2] - a[2]*b[1],
+		a[2]*b[0] - a[0]*b[2],
+		a[0]*b[1] - a[1]*b[0],
+	}
+}
+
+func dot(a, b [3]float64) float64 { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] }
diff --git a/internal/numerics/vec.go b/internal/numerics/vec.go
index 08d22a5..8247282 100644
--- a/internal/numerics/vec.go
+++ b/internal/numerics/vec.go
@@ -26,16 +26,6 @@ func PyMod(a, b float64) float64 {
 	return r
 }
 
-// GeoAdd returns y + k*dy with longitude wrapped to [0, 360). Latitude and
-// altitude accumulate linearly. This is the integrator's state-update step.
-func GeoAdd(y GeoVec, k float64, dy GeoVec) GeoVec {
-	return GeoVec{
-		Lat:      y.Lat + k*dy.Lat,
-		Lng:      PyMod(y.Lng+k*dy.Lng, 360),
-		Altitude: y.Altitude + k*dy.Altitude,
-	}
-}
-
 // GeoLerp linearly interpolates two geographic states by parameter l in
 // [0, 1]. Longitude takes the shorter great-circle arc.
 func GeoLerp(a, b GeoVec, l float64) GeoVec {
@@ -65,25 +55,6 @@ 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
-}