Record the measured validation and performance figures in the readme

The acceptance criteria ask for numbers rather than claims, and every row but
the three still-unbuilt ones has been measured. Splits validation into the
Orekit comparison and the analytic invariants, and gives the performance table
step counts alongside times so a ratio of times can be read as per-step cost.

The three invariant tests and the MSIS drag comparison now print what they
achieved, so the readme's figures can be reproduced by running them.

Also corrects the ordinary-diffeq version in the foundations table, and the
claim that the fixture generator produces observations, which it does not yet.

readme.md

@@ -35,17 +35,17 @@ These are depended on, not reimplemented.
 |---|---|---|
 | [`hifitime`](https://crates.io/crates/hifitime) 4.3 | MPL-2.0 | TAI/TT/ET/TDB/UTC/GPS time scales, leap seconds |
 | [`anise`](https://crates.io/crates/anise) 0.10 | MPL-2.0 | Reference frames, planetary ephemerides, SPICE-format kernels |
-| [`ordinary-diffeq`](https://crates.io/crates/ordinary-diffeq) 0.4 | MIT | Adaptive integration, dense output, callbacks, autodiff |
+| [`ordinary-diffeq`](https://crates.io/crates/ordinary-diffeq) 0.8 | MIT | Adaptive integration, dense output, callbacks, autodiff |
 
 ## Validation
 
 Before we can make a highly-performant propagator, we have to make a *correct* propagator.
 
 1. **Orekit reference fixtures.** A committed, rerunnable generation script
-   produces reference trajectories and observations; the data lands under
+   produces reference trajectories; the data lands under `tests/fixtures/` with the
-   `tests/fixtures/` with the Orekit version, force-model configuration and
+   Orekit version, force-model configuration and generation date recorded beside it,
-   generation date recorded beside it; Rust tests assert agreement within
+   and Rust tests assert agreement within declared tolerances. Observations are still
-   declared tolerances.
+   to come.
 2. **Analytic invariants.**: ie, "does a two-body orbit close back in on itself after one-period?" and similar conservation-tests, meant to run fast as unit tests to ensure that nothing has broken.
 3. **ANISE cross-checks.** Frame transformations and ephemeris lookups against
    ANISE's own validated values, confirming that we are calling it correctly.
@@ -54,28 +54,89 @@ Every fixture comparison matches gravity degree and order, density model,
 third-body ephemeris source and integrator tolerances; anything that cannot be
 matched is documented and its contribution bounded.
 
+### Against Orekit
+
+One force at a time over a one-day low Earth arc, worst disagreement across all 1441
+states, against Orekit 13.1.8. Each force is isolated so that a disagreement is
+attributable; the figures are read against the point-mass row, where there is no model
+left to mismatch and what is measured is the comparison itself — frames, time scale,
+units and `mu`.
+
 | Case | Target | Achieved | Test |
 |---|---|---|---|
-| 1-day LEO arc, full matched force model | 1 m position, 1 mm/s velocity | not yet measured | — |
+| Point mass | 1 m, 1 mm/s | 0.77 mm, 8.2e-7 m/s | `point_mass_arc` |
+| Sun and Moon, over a point mass | 1 m, 1 mm/s | 0.84 mm, 8.5e-7 m/s | `sun_and_moon` |
+| Cannonball radiation pressure | 1 m, 1 mm/s | 1.00 mm, 1.0e-6 m/s | `radiation_pressure` |
+| EGM96 8x8 | 1 m, 1 mm/s | 1.09 cm, 9.1e-6 m/s | `low_degree_field` |
+| EGM96 20x20 | 1 m, 1 mm/s | 1.10 cm, 9.0e-6 m/s | `higher_degree_field` |
+| Drag, Harris-Priester | 1% of the drag effect | 7.86 m, 8.9e-3 m/s — 0.025% | `drag` |
+| Drag, NRLMSISE-00 | 1e-3% of the drag effect | 2.7 cm, 3.1e-5 m/s — 0.0001% | `msis_drag` |
 | 6-day LEO arc | to be set from the 1-day result | not yet measured | — |
 | GEO arc, including a time-varying-area change | to be set | not yet measured | — |
-| Two-body energy and angular momentum | integrator tolerance | not yet measured | — |
-| Two-body closure over N periods | integrator tolerance | not yet measured | — |
-| J2 secular nodal regression and apsidal rotation | to be set | not yet measured | — |
 | RA/DEC, range, range-rate, az/el against Orekit | to be set | not yet measured | — |
 
+Tests are in `tests/orekit_comparison.rs`; run them with `--nocapture` to see the
+figures above printed.
+
+The two drag rows are bounded as a fraction of what drag does over the arc rather than
+in meters, because the meters depend on the arc and the fraction does not. They are also
+the one place the two sides knowingly model different things: both read the same
+tabulated Harris-Priester profiles with the same bulge lag and exponent, and the
+densities agree to 0.03%, but Orekit interpolates the table exponentially as published
+while this crate fits a monotone cubic to the logarithm so that the altitude derivative
+exists. NRLMSISE-00, where both sides implement the same reference, agrees 250 times
+closer.
+
+### Analytic invariants
+
+| Case | Target | Achieved | Test |
+|---|---|---|---|
+| Two-body closure over one period | integrator tolerance | 9.8e-7 m, 1.1e-9 m/s | `orbit_closes_over_one_period` |
+| Two-body energy and angular momentum, five periods | integrator tolerance | 1.1e-13, 5.3e-14 relative | `conserves_energy_and_angular_momentum` |
+| J2 secular nodal regression and apsidal rotation | 1% | 3.8e-3, 6.2e-3 relative | `j2_secular_rates_match_the_closed_form` |
+
+The J2 residual is systematic rather than noise: the closed form is first order and
+describes mean elements where the fit is on osculating ones, and the field's symmetry
+axis is the true pole while the inclination is measured in J2000. It is a check against a
+sign, normalization or coefficient error, all of which would miss by far more than a
+percent. The tight comparison is the Orekit one above.
+
 ## Performance
 
 Measured against `nyx-space` and Orekit on the same benchmarks with matched force
 models, on the same hardware. Neither publishes usable figures, so the references
 are measured here rather than cited.
 
+The scenarios are the ones the fixtures above were generated from, so the trajectories
+being timed are known to agree. One day of arc, steps and milliseconds, measured
+2026-10-02. Orekit runs under a JVM warmed ten times before anything is recorded, and
+the minimum of thirty runs is reported rather than the mean.
+
 | Benchmark | astrokit | nyx-space | Orekit |
 |---|---|---|---|
-| Two-body throughput (state-steps/s) | not yet measured | — | — |
+| Two-body | 1752 steps, 0.49 ms | not yet measured | 1932 steps, 4.60 ms |
-| Full-force LEO, 1 day | not yet measured | — | — |
+| EGM96 8x8 | 1968 steps, 31.9 ms | not yet measured | 1938 steps, 34.2 ms |
+| EGM96 20x20 | 2473 steps, 54.5 ms | not yet measured | 1940 steps, 52.9 ms |
+| Sun and Moon | 1794 steps, 35.9 ms | not yet measured | 1932 steps, 36.1 ms |
+| Radiation pressure | 2067 steps, 28.0 ms | not yet measured | 2035 steps, 133.6 ms |
+| Drag, Harris-Priester | 1829 steps, 49.7 ms | not yet measured | 2013 steps, 116.8 ms |
+| Drag, NRLMSISE-00 | 1859 steps, 52.0 ms | not yet measured | 2092 steps, 357.2 ms |
 | Ensemble of N states, scaling vs. cores | not yet measured | — | — |
 
+Step counts are reported alongside the times because the tolerances are not literally
+equal: Orekit is asked for a nanometer of position accuracy and this runs at
+`Accuracy::High`. The two do comparable work — within a quarter of each other on every
+row — so a ratio of times is close to a ratio of per-step cost. A future scenario whose
+step counts diverge needs them brought back together before its ratio means anything.
+
+The radiation row is not like for like however the steps come out. Orekit's radiation
+pressure registers eclipse event detectors, so its integrator locates the shadow
+boundaries by root finding, and this crate integrates straight through them. That is
+work Orekit does and this does not.
+
+Reproduce with `cargo bench --bench reference_scenarios` and
+`scripts/generate-orekit-fixtures.sh --time`.
+
 ## Kernels
 
 Ephemerides and Earth orientation come from NAIF kernels, which are too large to

src/dynamics/models/gravity.rs

@@ -230,6 +230,11 @@ mod tests {
         // approximating.
         assert_relative_eq!(node_rate, expected_node, max_relative = 1e-2);
         assert_relative_eq!(periapsis_rate, expected_periapsis, max_relative = 1e-2);
+        println!(
+            "J2 secular rates: node {:.2e}, periapsis {:.2e} relative",
+            (node_rate - expected_node).abs() / expected_node.abs(),
+            (periapsis_rate - expected_periapsis).abs() / expected_periapsis.abs()
+        );
 
         // And the signs, which first-order theory gets exactly right.
         assert!(node_rate < 0.0, "a prograde node should regress");

src/propagation.rs

@@ -602,6 +602,7 @@ mod tests {
             velocity_error < 1e-8,
             "velocity error {velocity_error:e} m/s"
         );
+        println!("one-period closure: {position_error:.2e} m, {velocity_error:.2e} m/s");
     }
 
     /// Two-body motion conserves specific energy and angular momentum. Both are
@@ -620,11 +621,18 @@ mod tests {
 
         let e0 = energy(&leo());
         let h0 = momentum(&leo());
+        let (mut worst_energy, mut worst_momentum) = (0.0_f64, 0.0_f64);
         for state in &solution.states {
             let y = Orbit.state(*state);
             assert_relative_eq!(energy(&y), e0, max_relative = 1e-12);
             assert_relative_eq!(momentum(&y).norm(), h0.norm(), max_relative = 1e-12);
+            worst_energy = worst_energy.max((energy(&y) - e0).abs() / e0.abs());
+            worst_momentum =
+                worst_momentum.max((momentum(&y).norm() - h0.norm()).abs() / h0.norm());
         }
+        println!(
+            "five-period drift: energy {worst_energy:.2e}, angular momentum {worst_momentum:.2e} relative"
+        );
     }
 
     /// The same spec and the same dynamics, propagated on three layouts of three

tests/orekit_comparison.rs

@@ -335,11 +335,13 @@ fn test_msis_drag_agrees_with_orekit() {
     .accuracy(Accuracy::High);
 
     let solution = propagation.run();
-    let mut position = 0.0_f64;
+    let (mut position, mut velocity) = (0.0_f64, 0.0_f64);
     for state in &arc.states {
         let ours = Orbit.state(solution.interpolate(state.elapsed));
         position =
             position.max((ours.position() - Vector3::from_column_slice(&state.position)).norm());
+        velocity =
+            velocity.max((ours.velocity() - Vector3::from_column_slice(&state.velocity)).norm());
     }
 
     let drag_free = {
@@ -353,7 +355,7 @@ fn test_msis_drag_agrees_with_orekit() {
     let effect = (Vector3::from_column_slice(&last.position) - drag_free).norm();
 
     println!(
-        "msis drag: {position:.3e} m, {:.4}% of a {effect:.0} m effect",
+        "msis drag: {position:.3e} m, {velocity:.3e} m/s, {:.4}% of a {effect:.0} m effect",
         position / effect * 100.0
     );
     assert!(