Orekit MSIS drag fixtures

tests/orekit_comparison.rs

@@ -293,3 +293,72 @@ fn test_drag_agrees_with_orekit() {
         position / effect * 100.0
     );
 }
+
+/// Drag through NRLMSISE-00, which both tools implement from the same reference.
+///
+/// Unlike Harris-Priester there is no intended difference in the model here, so this is
+/// the measurement of how close two independent implementations get when they are given
+/// very nearly the same weather. The drivers are not quite identical: for this epoch the
+/// daily flux and every three-hourly ap agree exactly, but CSSI publishes the 81-day
+/// average rounded to a tenth, 154.3 against the 154.3333 this crate computes from the
+/// daily values.
+#[cfg(feature = "msis")]
+#[test]
+fn test_msis_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/msis-drag-low-earth.json");
+    let arc = ReferenceArc::load(path).expect("the fixture is committed");
+    assert!(arc.setup.forces.iter().any(|f| f.contains("NRLMSISE00")));
+
+    let start = &arc.states[0];
+    let initial = Orbit
+        .with_position(Vector3::from_column_slice(&start.position))
+        .with_velocity(Vector3::from_column_slice(&start.velocity));
+    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
+    let mut propagation = Propagation::new(
+        epoch,
+        initial,
+        &[
+            TwoBody { body: Body::EARTH },
+            Drag {
+                density: DensityModel::Nrlmsise00,
+            },
+        ],
+    )
+    .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 = 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());
+    }
+
+    let drag_free = {
+        let mut plain = Propagation::new(epoch, initial, &[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!(
+        "msis drag: {position:.3e} m, {:.4}% of a {effect:.0} m effect",
+        position / effect * 100.0
+    );
+    assert!(
+        position < 1.0e-3 * effect,
+        "worst disagreement {position:.3e} m is {:.4}% of the {effect:.0} m drag moves it",
+        position / effect * 100.0
+    );
+}