Orekit drag fixtures

tests/orekit_comparison.rs

@@ -211,3 +211,85 @@ fn test_radiation_pressure_agrees_with_orekit() {
         "worst position disagreement {position:.3e} m over the day"
     );
 }
+
+/// Drag on a cannonball through Harris-Priester, over a point mass, at 400 km.
+///
+/// This is the one scenario where the two sides are *known* to be modeling different
+/// things, and the comparison measures how much that costs. Both read the same tabulated
+/// profiles with the same 30 degree bulge lag and the same exponent, but Orekit
+/// interpolates the table exponentially, as published, while this crate fits a monotone
+/// cubic to the logarithm of the density so that the altitude derivative exists. That is
+/// worth about 0.2% of the density in this band.
+#[test]
+fn test_drag_agrees_with_orekit() {
+    let Some(almanac) = almanac() else {
+        eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+        return;
+    };
+    let path = Path::new("tests/fixtures/orekit/drag-low-earth.json");
+    let arc = ReferenceArc::load(path).expect("the fixture is committed");
+    assert!(arc
+        .setup
+        .forces
+        .iter()
+        .any(|f| f.contains("HarrisPriester")));
+
+    let start = &arc.states[0];
+    let mut propagation = Propagation::new(
+        Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
+        Orbit
+            .with_position(Vector3::from_column_slice(&start.position))
+            .with_velocity(Vector3::from_column_slice(&start.velocity)),
+        &[
+            TwoBody { body: Body::EARTH },
+            Drag {
+                density: DensityModel::harris_priester(),
+            },
+        ],
+    )
+    .frame(Frame::EARTH_J2000)
+    .spacecraft(Spacecraft::cannonball(500.0, 5.0, 2.2, 1.3))
+    .almanac(almanac)
+    .duration(Duration::from_seconds(arc.duration()))
+    .accuracy(Accuracy::High);
+
+    let solution = propagation.run();
+    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());
+    }
+
+    // Bounded as a fraction of what drag does, not as a number of meters. The
+    // disagreement is a density difference of a few hundredths of a percent and the
+    // along-track position integrates it, so the meters depend on the arc while the
+    // fraction does not.
+    let drag_free = {
+        let mut plain = Propagation::new(
+            Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
+            Orbit
+                .with_position(Vector3::from_column_slice(&start.position))
+                .with_velocity(Vector3::from_column_slice(&start.velocity)),
+            &[TwoBody { body: Body::EARTH }],
+        )
+        .frame(Frame::EARTH_J2000)
+        .duration(Duration::from_seconds(arc.duration()))
+        .accuracy(Accuracy::High);
+        Orbit.state(*plain.run().states.last().unwrap()).position()
+    };
+    let last = arc.states.last().expect("checked on load");
+    let effect = (Vector3::from_column_slice(&last.position) - drag_free).norm();
+
+    println!(
+        "drag: {position:.3e} m, {velocity:.3e} m/s, {:.4}% of a {effect:.0} m effect",
+        position / effect * 100.0
+    );
+    assert!(
+        position < 5.0e-4 * effect,
+        "worst disagreement {position:.3e} m is {:.4}% of the {effect:.0} m drag moves it",
+        position / effect * 100.0
+    );
+}