v0.12.0
This release lands Kshana's first non-analytic orbit propagator — a Cowell
integrator with a hierarchical six-perturbation force model (two-body + J2–J6 zonal +
epoch-driven Sun/Moon third body + solar-radiation pressure with a conical
umbra/penumbra shadow + atmospheric drag + the post-Newtonian Schwarzschild relativistic
correction) driven by a choice of two adaptive integrators (RK4 step-doubling and the
Dormand–Prince RK5(4) embedded pair) — alongside a maneuver / trajectory-design layer
(impulsive and finite burns, an Izzo Lambert solver, and a porkchop sweep), a
gravity-map-matching alt-PNT layer that recovers a 60-minute GPS-denied track to under
500 m, a batch + sequential orbit-determination pipeline, and a full 17-state
tightly-coupled GNSS/INS UKF with quantum-CAI dead-reckoning. Every numerical capability
is pinned against analytic truth or a hand-derived closed form; the off-by-default
perturbations leave the released goldens untouched.
Added
- Post-Newtonian (Schwarzschild) relativistic correction (
forces::relativistic_accel+
propagator::ForceModel::relativity). Adds the dominant general-relativistic perturbation on a
near-Earth orbit — the leading driver of the relativistic perigee advance — in the IERS /
Montenbruck–Gillβ = γ = 1forma = (μ/c²r³)·{[4μ/r − v²]·r + 4(r·v)·v}. Like atmospheric
drag it is velocity-dependent, so it rides the(r, v)integrator RHS via
[accel_rv], opt-in and off by default. Validated self-contained: on a circular orbit it
collapses to the closed form3μ²/(c²r³)·r̂(purely radial and outward, off-axis components
exactly zero); its ratio to two-body is the textbook≈1.9·10⁻⁹at LEO (theμ/(c²r)
signature); a radial-velocity case matches the hand-simplifiedμ(4μ + 3v²r)/(c²r³); and in the
propagator it perturbs the orbit without dissipating it — the semi-major axis is conserved to
well under a metre/day, the structural opposite of drag's monotonic decay. Because it is off by
default the two-body/J2/zonal goldens are untouched. PPN-parameter (β,γ) tuning and the
Lense–Thirring frame-dragging term remain follow-ons. - Conical umbra+penumbra shadow model (
forces::conical_shadow), now used by solar-radiation
pressure. Upgrades the binary umbral-cylinder eclipse to a smoothν ∈ [0,1]factor: the Sun
and Earth are modelled as disks of apparent angular radiia = asin(R☉/d☉),b = asin(Rₑ/|r|)
with apparent centre separationc, andνis one minus the fraction of the Sun's disk occulted
by the Earth's disk (the circle–circle lens-overlap area) — full sun forc ≥ a+b, total umbra
forc ≤ b−a, annular forc ≤ a−b, and a continuous penumbra in between.srp_accelnow uses
it, so the SRP force tapers smoothly through eclipse instead of switching on/off. Adds the IAU
nominalforces::SOLAR_RADIUS. Validated self-contained:ν = 1in full sun andν = 0deep in
the umbra (exact), a smooth monotonic penumbra (νrises 0 → ~½ atc = b→ 1 across the
[b−a, b+a]band), and the conical penumbra extends beyond the umbral cylinder (a point the
binary cylinder calls fully lit is0 < ν < 1for the cone). The simplercylindrical_shadow
remains available; solar limb darkening and the oblate-Earth shadow remain follow-ons. - Dormand–Prince RK5(4) embedded integrator (
integrator::dopri54_step/
integrator::integrate_dopri+propagator::propagate_dopri). Adds the standard
Dormand–Prince (1980) embedded Butcher-tableau pair alongside the existing RK4 step-doubling
driver: seven FSAL stages yield a 5th-order solution and a 4th-order error estimate from one set
of evaluations (7 vs 11 function calls per step), a cheaper local-error estimate. The adaptive
driver reuses the same RMS-error norm and0.9·(1/err)^(1/5)step controller, so it is a drop-in
alternative;propagator::propagate_dopriexposes it on the orbit force model. Validated
self-contained: the embedded error estimate is O(h⁵) (halving the step cuts it ~32×); DP5(4)
integratesy' = ytoeand the harmonic oscillator over 50 periods conserving energy to
<1e-6; it reaches the same endpoint at the same tolerance in fewer function evaluations than
step doubling (without sacrificing accuracy); andpropagate_dopriclears the same analytic-truth
gate as the RK4 path — sub-metre against the exact universal-variable Kepler solution over a
24 h LEO orbit — while the two drivers agree to <1 m on a J2..J6 orbit (no closed form). Higher
embedded pairs (RKF7(8) / DOP853) remain a follow-on. - Atmospheric drag wired into the propagator as its first velocity-dependent force
(forces::atmospheric_density+forces::drag_accel+propagator::ForceModel::drag). Adds
the Vallado Table 8-4 piecewise-exponential atmosphereρ = ρ0·exp(−(h−h0)/H)(28 bands from
sea level past 1000 km, clamped below the surface) and the quadratic drag
a = −½ · ρ(h) · (C_D·A/m) · |v_rel| · v_relagainst the co-rotating atmosphere
v_rel = v − ωₑ ẑ × r(forces::EARTH_ROTATION_RATE = 7.2921151467e-5). Because drag depends on
velocity,ForceModelgains a newaccel_rv(t, r, v)and the integrator RHS now passes velocity
(f(t,[r;v]) = [v; a(t,r,v)]); the position-onlyaccel_atis unchanged, so the conservative
terms and goldens are untouched. Validated self-contained: the density anchors at the
1.225 kg/m³ sea-level value, clamps below the surface, decreases monotonically through LEO,
sits in the solar-mean ~1e-12 kg/m³ band at 400 km, and its recovered local scale height
(≈ 58 km at 400 km) is physical; drag opposes the co-rotating relative velocity at the
~2e-6 m/s² LEO magnitude forC_D·A/m = 0.02 m²/kg; and — the key signature — drag is
dissipative: a 300 km orbit loses specific energy monotonically and its semi-major axis
decays a bounded ~km/day, where the vacuum baseline conserves energy to <1e-9. The
NRLMSISE-00 thermospheric density (the < 5 % drag-density clause) remains a follow-on. - Solar-radiation pressure wired epoch-driven into the propagator force model
(forces::srp_accel+propagator::ForceModel::solar_radiation). Adds the cannonball SRP
modela = ν · P☉ · cᵣ · (A/m) · (AU/d)² · d̂with a cylindrical-shadow eclipse factor
(forces::cylindrical_shadow, ν ∈ {0,1}): the radiation pressureP☉ = Φ☉/cfrom the modern
1361 W/m² total solar irradiance (≈ 4.5398·10⁻⁶ N/m²), the inverse-square(AU/d)²flux fall-off,
and the radial push away from the Sun. It rides the same epoch-driven RHS as the third
body, sampling theephemSun once at the advanced epochepoch_jd_tt + t/86400shared between
the Sun third body and SRP. Composable:
with_zonals_j2_j6().third_body(true, true, epoch).solar_radiation(1.5, 0.02). Validated
self-contained against hand-derived signatures: the 1-AU radiation pressure pins to its textbook
≈ 4.5398·10⁻⁶ N/m²; a fully-lit LEO sat's SRP is bit-identical to the cannonball formula,
points away from the Sun, and sits in the ~1.36·10⁻⁷ m/s² band for cᵣ = 1.5, A/m = 0.02
m²/kg; doubling the Sun distance quarters the magnitude (inverse-square); the cylindrical
shadow eclipses only the umbral cylinder (anti-sunward and within one Earth radius of the
Earth–Sun line) and yields exactly zero SRP in eclipse; and in the propagator SRP perturbs
a LEO orbit by a small bounded amount that scales ~linearly with A/m — while a model with no
perturbations stays bit-for-bit time-independent, leaving the two-body/J2/zonal goldens untouched.
The conical umbra/penumbra (smooth ν ∈ [0,1]), atmospheric drag, and external GMAT/Orekit
cross-validation remain follow-ons. - Epoch-driven Sun/Moon third body wired into the time-varying propagator RHS
(propagator::ForceModel::third_body/accel_at). The third-body perturbation is no longer a
standalone force term — it is now integrated by the Cowell propagator as a genuinely time-varying
force: each RHS evaluation samples theephemSun/Moon positions at the advanced epoch
epoch_jd_tt + t/86400(reusingprecession::julian_centuries_ttfor the day↔century
conversion), so the perturbers move along their orbits during the integration rather than being
frozen at the start. Composable with any gravity model
(ForceModel::with_zonals_j2_j6().third_body(true, true, epoch)). Validated self-contained:
the RHS Sun term is bit-identical tothird_body_accelevaluated at the ephemeris position for
that instant at botht = 0andt = 1 day(proving the 86400 s ↔ 1 day ↔ 1/36525 century
wiring exactly), the perturber advances ~2.6·10⁹ m/day between samples (not frozen), the
instantaneous LEO tidal magnitudes hit the textbook ~5·10⁻⁷ m/s² (Sun) and ~1.1·10⁻⁶ m/s²
(Moon, ≈ 2× the Sun) bands, each body measurably perturbs the day-long trajectory while staying
bounded, and the same initial state propagated at epochs a quarter-year apart yields a
different trajectory (the tidal axis rotates 90°) — while a model with neither body enabled is
bit-for-bit time-independent, leaving the two-body/J2/zonal goldens untouched. DE-grade ephemeris
accuracy and external GMAT/Orekit cross-validation remain follow-ons. - Low-precision Moon ephemeris (
ephem::moon_position), completing the Sun/Moon third-body pair.
Adds the Montenbruck & Gill low-precision lunar series (§3.3.2) alongside the Sun model, so the
body-agnosticforces::third_body_accelcan now be driven by either luminary with no external
DE/SPK kernel. Validated self-contained against hand-derived lunar signatures: the geocentric
distance stays inside the real perigee/apogee envelope (~356 500–406 700 km) over a month and its
monthly mean recovers the ~384 400 km semi-major axis; the ecliptic latitude never exceeds the
~5.3° lunar-orbit inclination (checked by projecting onto the ecliptic pole in equatorial
coordinates, validating the latitude series and the obliquity rotation together); the Moon's
direction returns to within 1° after one sidereal month (27.3217 d) and its daily motion stays
in the physical 12–15°/day band; and the lunar third-body perturbation on a LEO satellite has the
textbook ~1.1·10⁻⁶ m/s² magnitude (≈ twice the Sun's). DE-grade position accuracy, atmospheric
drag, and SRP remain follow-ons. - Third-body (Sun) gravity with a built-in low-precision ephemeris (
forces::third_body_accel,
ephem::sun_position). Adds the third-body perturbation to the force model:
a = GM₃·((s−r)/|s−r|³ − s/|s|³)(direct attraction minus the indirect term the geocentric
frame must subtract), with the Sun position supplied by the newephemmodule's
Montenbruck & Gill low-precision analytical series — no external DE/SPK kernel needed for a
low-fidelity run. Validated self-contained: the acceleration matches the exact gradient of its
own disturbing potential (third_body_potential), the perturbation vanishes at the geocentre
and has the textbook ~5·10⁻⁷ m/s² magnitude on a LEO satellite, and the Sun ephemeris hits
hand-derived J2000 anchors — perihelion distance ≈ 1.471·10¹¹ m, declination ≈ −23° near the
December solstice, an apparent motion of ≈ 1°/day (≈ 90° per quarter-year), and a distance
that stays inside the 0.983–1.017 AU Earth-orbit envelope across a full year. Delivers the
third-body half of the numerical-propagator milestone's force-model step (the Moon is delivered in a
companion entry above); DE-grade position accuracy, atmospheric drag, and SRP remain follow-ons. - J2–J6 zonal-harmonic force model (
forces::zonal_accel/zonal_potential). Extends the
Cowell propagator's force model beyond J2 to the full Earth zonal field through degree 6 (the
standard published EGM-96 unnormalisedJ2..J6), wired into the propagator as
ForceModel::with_zonals_j2_j6(). The acceleration is the exact analytic gradient of the zonal
disturbing potentialR(r) = −(μ/r)·Σ Jₙ(Re/r)ⁿPₙ(z/r)(Legendre polynomials by upward recurrence),
validated three independent ways: it reduces to the 666-vector-validatedj2_accelto machine
precision when restricted to[J2]; it matches the numerical gradient of its own potential
through the full J2..J6 field (the conservative-field gold-standard check); and the oddJ3vs even
J2/J4..J6terms exhibit their characteristic north–south (anti)symmetry underz → −z— the
pear-shape asymmetry. A propagated J2..J6 orbit conserves total energy (kinetic + central + zonal
potential) to ~1e-8 over a day and perturbs the J2-only orbit by a small non-zero amount. This
delivers step-2 ("J2–J6 zonal harmonics") of the numerical-propagator milestone; the high-degree EGM
tesseral field, drag, SRP, third-body, and external GMAT/Orekit cross-validation remain follow-ons. - Numerical (Cowell) orbit propagator (
src/propagator.rs). Kshana's first non-analytic
propagator (the rest of the orbit stack is analytic SGP4/SDP4): it wires the two-body + J2 force
model (src/forces.rs) into the adaptive step-doubling RK4 driver (src/integrator.rs) as
f(t,[r;v]) = [v; a(r)], with aForceModeltoggle. Validated against analytic truth that is
stronger than a numerical cross-tool would be: the unperturbed orbit reproduces the exact
universal-variable Kepler solution to sub-metre over a 24-hour LEO orbit (a tighter gate than
the "vs a numerical reference < 10 m" the milestone phrases), specific energy and angular momentum
conserve to ~1e-9 relative, and the J2 nodal regression reproduces the closed-formj2_secular_rates
to first-order theory (within 2 %, the O(J2²) residual). Also addssolve_kepler_checked, a Newton
solver for Kepler's equation that returnsErrinstead of a silently-wrong answer when it fails
to converge within a bounded iteration budget (the near-perigeee = 0.999case). Honest scope: the
force model is two-body + J2 only — the high-degree EGM tesseral field (200×200 + loader), drag
(NRLMSISE-00), SRP, third-body forces, and an external GMAT/Orekit cross-validation remain follow-ons. - 60-minute GPS-denied gravity-map matching to < 500 m (
run_gps_denied_gravity_nav).
Deepens the alt-PNT layer to the ESA NAVISP Quantum Wayfarer validation target: a vehicle
flies a ~700 km track for a full one-hour GNSS outage — its inertial solution drifting to
≈ 70 km — and a cold-atom gravimeter plus a hierarchical coarse-to-fine particle/grid
matcher recovers the constant INS drift to ≈ 145 m (< 500 m), a > 480× cut. The gravimeter's
real white-noise floor is injected as a deterministic seeded sequence, so the matcher is
never handed noise-free truth yet the run is exactly reproducible (verified bit-identical, and
stable to a few metres across noise realisations). A regression-grade test shows the refinement
is necessary — a single coarse grid stalls at ~2 km, only the three-stage refinement breaks
the 500 m barrier. New committed scenarioscenarios/gps-denied-gravity-nav.toml. The
docs/CAPABILITY.mdrow stays honestly partial (still no bundled EGM2008 map) with its
evidence updated to the 60-min < 500 m result. Honest scope unchanged: low-degree
spherical-harmonic field + synthetic mascons; a Monte-Carlo over map-representation-error
realisations is a follow-on. - Overclaim ledger + regression guard (
docs/CLAIMS-VS-REALITY.md,tests/no_overclaims.rs).
Closes the honesty/de-claim track: the fourteen overclaims an earlier audit catalogued
(OC-0…OC-13) are now all GREEN — the strong claims (OC-0coupled clock+position Kalman,
OC-2jamming J/S→C/N₀→loss-of-lock,OC-7Mach–Zehnder CAI physics,OC-8ARAIM HPL/VPL)
are superseded by shipped, tested capabilities rather than softened wording, and the
remaining rows are de-claimed to match the code. A new CI test scans the live public surfaces
(README,CAPABILITY,GLOSSARY,web/) and fails if any retired bare overclaim phrase
reappears uncaveated, so a GREEN row cannot silently regress. The per-run "integrity" FoM stays
honestly labelled filter self-consistency (not aviation integrity); the real ARAIM HPL/VPL is
surfaced separately so the two are never conflated. - Gravity-map-matching navigation (GPS-denied alt-PNT). New
src/gravimeter.rsadds the
alt-PNT capability layer ESA NAVISP's Quantum Wayfarer / QT-CCI gravity-map-matching studies
call for: a cold-atom gravimeter measurement model whose white-noise floor is derived from
the CAI accelerometer ASD (σ = ASD/√τ); a low-degree, fully-normalised spherical-harmonic
gravity-anomaly field (validated against the closed-form Legendre functionsP̄₁₁=√3·cosφ,
P̄₂₀=(√5/2)(3sin²φ−1),P̄₂₂=(√15/2)cos²φand a hand-derived single-term anomaly of
1.897 mGal) plus synthetic mascons for the high-degree local features; and a
gravity-map-matching particle filter (composingmapmatch+particle_filter) that
recovers a GPS-denied track from the anomaly sequence it flies through. A committed NAVISP
benchmark (scenarios/gravity-map-nav.toml) cuts a ~73 km free-inertial drift to a few km.
Honest scope: Kshana does not bundle the full EGM2008 2190° coefficient set — the field is
low-degree + mascons, not a real high-resolution map; the EGM/EIGEN loader, magnetic map,
terrain-aided SLAM, and scenario-enginekind=wiring with an SVG drift chart remain follow-ons.
docs/CAPABILITY.md"Gravity-map / alt-PNT navigation" → partial. - Maneuver modeling and trajectory-design beachhead. New
src/maneuver.rsadds the first
trajectory-design layer above SGP4: impulsive ΔV nodes that apply a velocity discontinuity and
carry a 6×6 covariance forward (deterministic burn ⇒ identity state-transition; the
execution-error covariance rotates from the burn frame — ECI or LVLH — into the velocity block),
a finite-burn integration (constant thrust over a burn arc with mass as a state) whose achieved
ΔV is checked against the closed-form Tsiolkovsky rocket equation to better than 0.01 %, an
Izzo-2015 single-revolution Lambert solver (r1,r2, time-of-flight ⇒v1,v2),
an exact universal-variable Kepler propagator (two-body truth), and a porkchop sweep that
maps a launch-epoch × arrival-epoch grid to departure C3 and arrival V∞, emitted as a 2-D JSON
array for browser contour rendering. Validation is self-contained and stronger than a tutorial
read-out: every Lambert output is round-tripped through the Kepler propagator (it must land back
onr2), and the porkchop minimum is checked against the analytic Hohmann-transfer C3 floor for
two coplanar circular orbits. Kshana positions this as the performance-simulation layer above
GMAT/Orekit, not a replacement (multi-revolution branches and a real planetary ephemeris remain
out of scope). Ten tests. - Full 17-state tightly-coupled GNSS/INS UKF with quantum-CAI dead-reckoning. New
src/fusion/tightly_coupled17.rscarries the complete inertial-navigation state a
tightly-coupled filter estimates —[position, velocity, attitude-error, accelerometer bias, gyro bias, clock bias, clock drift](17 states) — propagated through the strapdown
mechanization driven by the measured specific force and angular rate, with the small-angle
C ≈ I + [ψ×]body→inertial rotation so attitude error couples into horizontal acceleration
(the standard INS tilt coupling). During a GNSS outage it coasts on the IMU alone; the
velocity-random-walk process noise is the cold-atom-interferometer accelerometer's derived
q_va(crate::inertial::quantum_imu), so the dead-reckoning drift is the quantum-sensor
limited one — a 120-second outage stays bounded to a few hundred metres versus the kilometres
a navigation-grade free INS would reach. The pseudorange/range-rate update runs through the
shared unscented filter (with α = 1 for well-conditioned weights at this state size). Five
tests: measurement-model identity, perfect-IMU constant-velocity integration, GNSS aiding,
accelerometer-bias estimation, and the CAI-limited 120-s outage benchmark.
Full changelog: CHANGELOG.md · Docs: README