Orekit fixture generation

tests/common/mod.rs

@@ -0,0 +1,163 @@
+//! The reference-trajectory format, and reading it.
+//!
+//! A fixture is one propagated arc from another tool, with enough of that tool's setup
+//! recorded alongside it to tell whether the comparison is fair. The format is this
+//! crate's to define — a generator conforms to it, not the other way round — so it is
+//! specified here and the loader is the specification's executable half.
+//!
+//! # Matching is the whole difficulty
+//!
+//! A disagreement between two propagators is only evidence of anything if they were asked
+//! the same question. Every one of these has to agree before a number means anything, and
+//! each has bitten someone:
+//!
+//! - **The gravity field, by name and truncation.** Not "EGM96 to degree 70" but the same
+//!   coefficient file: this crate found a 180x180 table that was not EGM96 and differed in
+//!   the fourth digit, which is 1.6 m/day of along-track at a low altitude.
+//! - **The tide convention.** Tide-free against zero-tide moves `C(2,0)` alone, and that
+//!   term dominates everything.
+//! - **The reference radius and `mu` the coefficients are defined against**, which are not
+//!   the body's WGS84 figures. This crate uses EGM96's own 6378136.3 m.
+//! - **The third-body ephemeris**, by source and not only by body. DE440 against DE430 is
+//!   a real difference at the meter level over days.
+//! - **Earth orientation.** Whether EOP corrections are applied at all, and from which
+//!   series, changes where the body-fixed frame points, which changes the gravity
+//!   gradient and the drag wind.
+//! - **The atmosphere**, including which solar and geomagnetic drivers were read, since a
+//!   density model given different weather is a different model.
+//! - **The drag and radiation coefficients and areas**, which tools variously take as
+//!   `Cd` and area separately or as the product.
+//! - **The integrator and its tolerances**, since a comparison at 1 m has to be made
+//!   between two arcs each converged well inside that.
+//!
+//! The configuration block below is deliberately free-form text per entry rather than a
+//! typed mirror of this crate's `Dynamics`: it records what the *other* tool was told, in
+//! that tool's own vocabulary, so that a mismatch is visible rather than silently
+//! translated away.
+
+use serde::{Deserialize, Serialize};
+
+/// One state along a reference arc.
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub struct ReferenceState {
+    /// Seconds of uniform time after the arc's epoch.
+    pub elapsed: f64,
+    /// Position in the arc's frame, m.
+    pub position: [f64; 3],
+    /// Velocity in the arc's frame, m/s.
+    pub velocity: [f64; 3],
+}
+
+/// What produced a reference arc, in the producing tool's own terms.
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct ReferenceSetup {
+    /// The tool, for example `orekit`.
+    pub tool: String,
+    /// Its version, exactly as the tool reports it.
+    pub version: String,
+    /// The date the fixture was generated, `YYYY-MM-DD`.
+    pub generated: String,
+    /// The force model as the tool was configured, one entry per force, in that tool's
+    /// vocabulary rather than translated into this crate's.
+    pub forces: Vec<String>,
+    /// The integrator and its tolerances, likewise verbatim.
+    pub integrator: String,
+    /// Anything else that had to be matched: Earth orientation, tide convention,
+    /// ephemeris source, and so on.
+    pub notes: Vec<String>,
+}
+
+/// One reference arc.
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct ReferenceArc {
+    /// A short name, which is also the fixture's file stem.
+    pub name: String,
+    /// What produced it.
+    pub setup: ReferenceSetup,
+    /// Epoch of the first state, ISO 8601 with its time scale.
+    pub epoch: String,
+    /// The frame every state is expressed in, for example `EME2000`.
+    pub frame: String,
+    /// The states, earliest first, the first being the initial condition.
+    pub states: Vec<ReferenceState>,
+}
+
+/// Why a fixture could not be used.
+#[derive(Debug)]
+pub enum FixtureError {
+    /// The file could not be read.
+    Unreadable(std::io::Error),
+    /// The file is not the format.
+    Malformed(serde_json::Error),
+    /// The file parses but says something impossible.
+    Invalid(String),
+}
+
+impl std::fmt::Display for FixtureError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Self::Unreadable(error) => write!(f, "could not read the fixture: {error}"),
+            Self::Malformed(error) => write!(f, "not a reference arc: {error}"),
+            Self::Invalid(why) => write!(f, "unusable reference arc: {why}"),
+        }
+    }
+}
+
+impl ReferenceArc {
+    /// Reads an arc, checking it says something usable.
+    ///
+    /// # Errors
+    ///
+    /// If the file is missing, is not the format, or is internally inconsistent.
+    pub fn load(path: &std::path::Path) -> Result<Self, FixtureError> {
+        let text = std::fs::read_to_string(path).map_err(FixtureError::Unreadable)?;
+        let arc: Self = serde_json::from_str(&text).map_err(FixtureError::Malformed)?;
+        arc.check()?;
+        Ok(arc)
+    }
+
+    /// The checks a fixture has to pass before it is worth comparing against.
+    ///
+    /// Not style: each of these would otherwise show up as a physics disagreement much
+    /// later. An arc whose times run backwards, or which records no setup, produces a
+    /// number that looks like a result.
+    fn check(&self) -> Result<(), FixtureError> {
+        let invalid = |why: String| Err(FixtureError::Invalid(why));
+        if self.states.len() < 2 {
+            return invalid(format!("{}: {} states", self.name, self.states.len()));
+        }
+        if self.states[0].elapsed != 0.0 {
+            return invalid(format!(
+                "{}: starts at {} s rather than at its own epoch",
+                self.name, self.states[0].elapsed
+            ));
+        }
+        if self.setup.forces.is_empty() {
+            return invalid(format!("{}: records no force model", self.name));
+        }
+        for pair in self.states.windows(2) {
+            if pair[1].elapsed <= pair[0].elapsed {
+                return invalid(format!(
+                    "{}: times are not increasing, {} then {}",
+                    self.name, pair[0].elapsed, pair[1].elapsed
+                ));
+            }
+        }
+        for state in &self.states {
+            let finite = state
+                .position
+                .iter()
+                .chain(&state.velocity)
+                .all(|x| x.is_finite());
+            if !finite {
+                return invalid(format!("{}: a state is not finite", self.name));
+            }
+        }
+        Ok(())
+    }
+
+    /// How long the arc covers, in seconds.
+    pub fn duration(&self) -> f64 {
+        self.states.last().map(|s| s.elapsed).unwrap_or(0.0)
+    }
+}

tests/orekit_comparison.rs

@@ -0,0 +1,71 @@
+//! This crate's propagator against Orekit's, on the fixtures in `tests/fixtures/orekit`.
+//!
+//! Generated by `scripts/generate-orekit-fixtures.sh`. Each fixture records the setup
+//! 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.
+
+mod common;
+
+use astrokit::prelude::*;
+use common::ReferenceArc;
+use nalgebra::Vector3;
+use std::path::Path;
+
+/// 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) {
+    let start = &arc.states[0];
+    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
+    let mut propagation = Propagation::new(
+        epoch,
+        Orbit
+            .with_position(Vector3::from_column_slice(&start.position))
+            .with_velocity(Vector3::from_column_slice(&start.velocity)),
+        forces,
+    )
+    .frame(Frame::EARTH_J2000)
+    .duration(Duration::from_seconds(arc.duration()))
+    .accuracy(Accuracy::High);
+
+    let solution = propagation.run();
+    let (mut worst_position, mut worst_velocity) = (0.0_f64, 0.0_f64);
+    for state in &arc.states {
+        let ours = Orbit.state(solution.interpolate(state.elapsed));
+        worst_position = worst_position
+            .max((ours.position() - Vector3::from_column_slice(&state.position)).norm());
+        worst_velocity = worst_velocity
+            .max((ours.velocity() - Vector3::from_column_slice(&state.velocity)).norm());
+    }
+    (worst_position, worst_velocity)
+}
+
+#[test]
+fn test_a_point_mass_arc_agrees_with_orekit() {
+    let path = Path::new("tests/fixtures/orekit/two-body-low-earth.json");
+    let arc = ReferenceArc::load(path).expect("the fixture is committed");
+
+    assert_eq!(arc.setup.tool, "orekit");
+    assert_eq!(arc.duration(), 86400.0);
+    assert_eq!(arc.states.len(), 1441);
+
+    let (position, velocity) = compare(&arc, &[TwoBody { body: Body::EARTH }]);
+
+    // 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:
+    // a later fixture that disagrees by more than this is disagreeing about physics.
+    assert!(
+        position < 1.0e-3,
+        "worst position disagreement {position:.3e} m over the day"
+    );
+    assert!(
+        velocity < 1.0e-6,
+        "worst velocity disagreement {velocity:.3e} m/s over the day"
+    );
+    println!("point mass: {position:.3e} m, {velocity:.3e} m/s");
+}

tests/reference_format.rs

@@ -0,0 +1,157 @@
+//! The reference-trajectory format: that it round trips, and that it refuses a fixture
+//! which would otherwise produce a number looking like a result.
+//!
+//! No fixture files are committed yet — generating them needs a JVM, which this machine
+//! does not have. What is tested here is the format itself, which a generator has to
+//! conform to.
+
+mod common;
+
+use common::{FixtureError, ReferenceArc, ReferenceSetup, ReferenceState};
+
+fn setup() -> ReferenceSetup {
+    ReferenceSetup {
+        tool: "orekit".to_string(),
+        version: "13.1.8".to_string(),
+        generated: "2026-10-02".to_string(),
+        forces: vec!["HolmesFeatherstoneAttractionModel(EGM96, 70, 70)".to_string()],
+        integrator: "DormandPrince853, abs 1e-12, rel 1e-12".to_string(),
+        notes: vec!["EOP: none applied".to_string()],
+    }
+}
+
+fn arc(states: Vec<ReferenceState>) -> ReferenceArc {
+    ReferenceArc {
+        name: "low-earth-one-day".to_string(),
+        setup: setup(),
+        epoch: "2026-01-01T00:00:00 UTC".to_string(),
+        frame: "EME2000".to_string(),
+        states,
+    }
+}
+
+fn state(elapsed: f64) -> ReferenceState {
+    ReferenceState {
+        elapsed,
+        position: [6.978e6, 0.0, 0.0],
+        velocity: [0.0, 7.56e3, 0.0],
+    }
+}
+
+fn write(arc: &ReferenceArc) -> tempfile_path::Written {
+    tempfile_path::write(&serde_json::to_string_pretty(arc).unwrap())
+}
+
+/// A minimal stand-in for a temporary file, so the test needs no extra dependency.
+mod tempfile_path {
+    use std::path::{Path, PathBuf};
+    use std::sync::atomic::{AtomicUsize, Ordering};
+
+    pub struct Written(PathBuf);
+
+    impl Written {
+        pub fn path(&self) -> &Path {
+            &self.0
+        }
+    }
+
+    impl Drop for Written {
+        fn drop(&mut self) {
+            let _ = std::fs::remove_file(&self.0);
+        }
+    }
+
+    pub fn write(text: &str) -> Written {
+        static COUNT: AtomicUsize = AtomicUsize::new(0);
+        let path = std::env::temp_dir().join(format!(
+            "astrokit-fixture-{}-{}.json",
+            std::process::id(),
+            COUNT.fetch_add(1, Ordering::Relaxed)
+        ));
+        std::fs::write(&path, text).expect("the temporary directory is writable");
+        Written(path)
+    }
+}
+
+#[test]
+fn test_an_arc_round_trips_through_the_format() {
+    let original = arc(vec![state(0.0), state(60.0), state(120.0)]);
+    let written = write(&original);
+    assert_eq!(ReferenceArc::load(written.path()).unwrap(), original);
+}
+
+#[test]
+fn test_the_setup_is_carried_verbatim() {
+    // The point of the block: it records what the *other* tool was told, in that tool's
+    // own words, so a mismatch is visible rather than translated away.
+    let original = arc(vec![state(0.0), state(60.0)]);
+    let written = write(&original);
+    let loaded = ReferenceArc::load(written.path()).unwrap();
+
+    assert_eq!(loaded.setup.tool, "orekit");
+    assert_eq!(loaded.setup.version, "13.1.8");
+    assert!(loaded.setup.forces[0].contains("EGM96"));
+    assert!(loaded.setup.integrator.contains("1e-12"));
+    assert_eq!(loaded.duration(), 60.0);
+}
+
+#[test]
+fn test_a_fixture_that_would_mislead_is_refused() {
+    let unusable = |arc: ReferenceArc, expected: &str| {
+        let written = write(&arc);
+        match ReferenceArc::load(written.path()) {
+            Err(FixtureError::Invalid(why)) => {
+                assert!(why.contains(expected), "expected {expected:?} in {why:?}")
+            }
+            other => panic!("should have been refused: {other:?}"),
+        }
+    };
+
+    // One state is not an arc.
+    unusable(arc(vec![state(0.0)]), "1 states");
+
+    // Not starting at its own epoch means every comparison is offset by an unknown.
+    unusable(
+        arc(vec![state(30.0), state(60.0)]),
+        "rather than at its own epoch",
+    );
+
+    // Times that do not advance mean the states cannot be matched to anything.
+    unusable(
+        arc(vec![state(0.0), state(60.0), state(60.0)]),
+        "not increasing",
+    );
+
+    // An arc with no recorded force model cannot be matched, so it cannot be compared.
+    let mut no_forces = arc(vec![state(0.0), state(60.0)]);
+    no_forces.setup.forces.clear();
+    unusable(no_forces, "records no force model");
+}
+
+/// A non-finite state cannot survive the format at all: JSON has no NaN, so it is written
+/// as `null` and refused on the way back in.
+#[test]
+fn test_a_non_finite_state_cannot_round_trip() {
+    let mut broken = arc(vec![state(0.0), state(60.0)]);
+    broken.states[1].position[0] = f64::NAN;
+    let written = write(&broken);
+    assert!(matches!(
+        ReferenceArc::load(written.path()),
+        Err(FixtureError::Malformed(_))
+    ));
+}
+
+#[test]
+fn test_a_malformed_file_is_told_apart_from_an_unusable_one() {
+    let written = tempfile_path::write("{\"name\": \"half a fixture\"}");
+    assert!(matches!(
+        ReferenceArc::load(written.path()),
+        Err(FixtureError::Malformed(_))
+    ));
+
+    let missing = std::path::Path::new("tests/fixtures/orekit/not-generated-yet.json");
+    assert!(matches!(
+        ReferenceArc::load(missing),
+        Err(FixtureError::Unreadable(_))
+    ));
+}