Orekit gravity fixtures

tests/orekit_comparison.rs

@@ -4,22 +4,34 @@
 //! Orekit was given, so a disagreement can be attributed to the models rather than left
 //! as a mystery; see `tests/common/mod.rs` for what has to match and why.
 //!
-//! The only fixture so far is a point mass. That is deliberate as a first step: with no
-//! field, no third body, no atmosphere and an inertial frame, there is nothing left to
-//! mismatch, so the figure it produces measures the comparison itself — frames, time
-//! scale, units and `mu` — rather than any model. Anything added after this is measured
-//! against a harness already known to be sound.
+//! The point-mass fixture comes first and is not about gravity at all: with no field, no
+//! third body, no atmosphere and an inertial frame, there is nothing left to mismatch, so
+//! the figure it produces measures the comparison itself — frames, time scale, units and
+//! `mu`. Every later scenario is read against that floor.
 
 mod common;
 
+use astrokit::kernels;
 use astrokit::prelude::*;
 use common::ReferenceArc;
 use nalgebra::Vector3;
 use std::path::Path;
+use std::sync::Arc;
+
+/// The kernels, or `None` if they have not been fetched.
+fn almanac() -> Option<Arc<anise::almanac::Almanac>> {
+    let directory = kernels::directory()
+        .unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
+    kernels::load(&directory).ok().map(Arc::new)
+}
 
 /// Propagates the fixture's initial state with this crate and returns the worst position
 /// and velocity disagreement over the arc.
-fn compare(arc: &ReferenceArc, forces: &[Dynamics]) -> (f64, f64) {
+fn compare(
+    arc: &ReferenceArc,
+    forces: &[Dynamics],
+    almanac: Option<Arc<anise::almanac::Almanac>>,
+) -> (f64, f64) {
     let start = &arc.states[0];
     let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
     let mut propagation = Propagation::new(
@@ -32,6 +44,9 @@ fn compare(arc: &ReferenceArc, forces: &[Dynamics]) -> (f64, f64) {
     .frame(Frame::EARTH_J2000)
     .duration(Duration::from_seconds(arc.duration()))
     .accuracy(Accuracy::High);
+    if let Some(almanac) = almanac {
+        propagation = propagation.almanac(almanac);
+    }
 
     let solution = propagation.run();
     let (mut worst_position, mut worst_velocity) = (0.0_f64, 0.0_f64);
@@ -54,7 +69,7 @@ fn test_a_point_mass_arc_agrees_with_orekit() {
     assert_eq!(arc.duration(), 86400.0);
     assert_eq!(arc.states.len(), 1441);
 
-    let (position, velocity) = compare(&arc, &[TwoBody { body: Body::EARTH }]);
+    let (position, velocity) = compare(&arc, &[TwoBody { body: Body::EARTH }], None);
 
     // Sub-millimeter over a day. Nothing in this case can disagree except the two
     // integrators and the arithmetic, so this is the floor the comparison can measure to:
@@ -69,3 +84,45 @@ fn test_a_point_mass_arc_agrees_with_orekit() {
     );
     println!("point mass: {position:.3e} m, {velocity:.3e} m/s");
 }
+
+/// EGM96, read by both sides from the same coefficient file.
+///
+/// The one knob matched by intent rather than by construction is the body-fixed frame:
+/// Orekit turns the Earth with ITRF under IERS 2010 conventions and the Earth orientation
+/// parameters in `orekit-data`, while this crate uses ITRF93 from a SPICE binary PCK. The
+/// agreement below is what that difference costs, together with everything else.
+fn gravity(degree: usize, allowed_meters: f64) {
+    let Some(almanac) = almanac() else {
+        eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+        return;
+    };
+    let path = format!("tests/fixtures/orekit/gravity-{degree}-low-earth.json");
+    let arc = ReferenceArc::load(Path::new(&path)).expect("the fixture is committed");
+    assert!(arc.setup.forces[0].contains("EGM96"));
+
+    let (position, velocity) = compare(
+        &arc,
+        &[Gravity {
+            field: GravityModel::Egm96,
+            degree,
+            order: degree,
+        }],
+        Some(almanac),
+    );
+    assert!(
+        position < allowed_meters,
+        "degree {degree}: worst position disagreement {position:.3e} m over the day"
+    );
+    println!("gravity {degree}: {position:.3e} m, {velocity:.3e} m/s");
+}
+
+#[test]
+fn test_a_low_degree_field_agrees_with_orekit() {
+    // An order of magnitude inside the one meter a day this is aiming at.
+    gravity(8, 0.1);
+}
+
+#[test]
+fn test_a_full_field_agrees_with_orekit() {
+    gravity(70, 0.2);
+}