Ground stations are in

tests/ground_station_passes.rs

@@ -0,0 +1,166 @@
+//! A ground station pass, located by the integrator from a callback a user wrote.
+//!
+//! An integration test rather than a unit one on purpose. It compiles against the crate's
+//! public API only, so it is the same thing a script outside the crate would write, and it
+//! fails if any of that surface stops being reachable.
+//!
+//! What it pins is the callback seam. Nothing inside the crate attaches a station event —
+//! passes are meant to be found by post-processing a finished solution — so without this
+//! test nothing would notice if `Assembly::context`, or the epoch, frame and almanac it
+//! reaches, stopped being public. That combination is what makes a user-supplied event
+//! possible at all.
+//!
+//! It imports from `std` and `astrokit` and nothing else, deliberately. A test target can
+//! reach every one of the crate's own dependencies directly, so importing `nalgebra` or
+//! `ordinary_diffeq` here would compile while leaving a consumer of the published crate
+//! unable to write the same thing.
+
+use std::path::Path;
+use std::sync::Arc;
+
+use astrokit::kernels;
+use astrokit::nalgebra::SVector;
+use astrokit::prelude::*;
+use astrokit::propagation::Assembly;
+
+/// Goldstone, at a latitude a 51.6 degree orbit passes over.
+fn goldstone() -> GroundStation {
+    GroundStation::new(
+        "Goldstone",
+        35.4267_f64.to_radians(),
+        -116.8900_f64.to_radians(),
+        1036.0,
+        Body::EARTH,
+    )
+}
+
+/// What a user script would write: note every acquisition and loss of signal.
+struct Passes {
+    station: GroundStation,
+    /// Time of each crossing, and the elevation above the mask the integrator left it at.
+    crossings: Vec<(f64, f64)>,
+}
+
+impl Passes {
+    /// How far the spacecraft is above the station's mask, the quantity whose sign changes
+    /// at a pass boundary.
+    ///
+    /// Everything it needs comes off the assembly's context: the start epoch to turn the
+    /// integrator's `t` into a real one, the frame the state is in, and the kernels to
+    /// rotate it into the frame the station is at rest in.
+    fn above_mask<T, const D: usize, L>(
+        &self,
+        t: T,
+        y: &SVector<T, D>,
+        assembly: &Assembly<T, D, L>,
+    ) -> T
+    where
+        T: Scalar,
+        L: Layout<D> + HasPosition<D> + HasVelocity<D>,
+    {
+        let context = assembly.context();
+        let when = time::advance(context.epoch.epoch(), t.re());
+        let inertial = Cartesian::new(L::position(y), L::velocity(y), context.frame);
+        let fixed = inertial
+            .to_frame(
+                Frame::EARTH_FIXED,
+                when,
+                context
+                    .almanac()
+                    .expect("the propagation was given kernels"),
+            )
+            .expect("an inertial state converts to the body-fixed frame");
+
+        self.station
+            .elevation(&fixed)
+            .expect("body-fixed by construction")
+            - T::lift(self.station.elevation_mask)
+    }
+}
+
+impl<T, const D: usize, L> Callback<T, D, Assembly<T, D, L>> for Passes
+where
+    T: Scalar,
+    L: Layout<D> + HasPosition<D> + HasVelocity<D>,
+{
+    fn event(&self, t: T, y: &SVector<T, D>, assembly: &Assembly<T, D, L>) -> T {
+        self.above_mask(t, y, assembly)
+    }
+
+    fn effect(&mut self, ode: &mut ODE<'_, T, D, Assembly<T, D, L>>) {
+        // The state is not changed: this callback only wants to be told where the boundary
+        // was, which is what makes it an observer rather than part of the dynamics.
+        let residual = self.above_mask(ode.t, &ode.y, &ode.params).re();
+        self.crossings.push((ode.t.re(), residual));
+    }
+}
+
+#[test]
+fn test_a_user_callback_locates_the_station_passes() {
+    let directory = kernels::directory()
+        .unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
+    let Ok(loaded) = kernels::load(&directory) else {
+        eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+        return;
+    };
+    let almanac = Arc::new(loaded);
+
+    let mut passes = Passes {
+        station: goldstone(),
+        crossings: Vec::new(),
+    };
+
+    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
+    let orbit = Keplerian::new(6.978e6, 1.0e-3, 0.9, 0.4, 0.2, 0.0, Frame::EARTH_J2000);
+    let arc = Duration::from_seconds(86400.0);
+
+    // Scoped, because the callback is mutably borrowed for as long as the propagation it
+    // is attached to, so what it recorded cannot be read until that borrow ends.
+    {
+        let mut propagation = Propagation::new(
+            epoch,
+            Orbit.with_keplerian(orbit),
+            &[TwoBody { body: Body::EARTH }],
+        )
+        .almanac(almanac)
+        .callback(&mut passes)
+        .duration(arc);
+        propagation.run();
+    }
+
+    let crossings = &passes.crossings;
+    assert!(
+        crossings.len() >= 4,
+        "only {} crossings over a day",
+        crossings.len()
+    );
+
+    // The point of locating a root rather than sampling: the integrator lands on the
+    // boundary, so the elevation there is the mask.
+    for (at, residual) in crossings {
+        assert!(
+            residual.abs() < 1.0e-9,
+            "{:.1} s is {:.3e} rad off the mask",
+            at,
+            residual
+        );
+    }
+
+    // Strictly ordered, and inside the arc.
+    for pair in crossings.windows(2) {
+        assert!(pair[1].0 > pair[0].0, "crossings out of order");
+    }
+    assert!(crossings[0].0 > 0.0 && crossings.last().unwrap().0 < arc.to_seconds());
+
+    // Consecutive crossings alternate between a pass and the wait for the next one, so the
+    // gaps should separate into minutes and roughly an orbit.
+    let gaps: Vec<f64> = crossings.windows(2).map(|p| p[1].0 - p[0].0).collect();
+    assert!(
+        gaps.iter().any(|gap| (60.0..1200.0).contains(gap)),
+        "no gap short enough to be a pass: {gaps:?}"
+    );
+    assert!(
+        gaps.iter().any(|gap| *gap > 3000.0),
+        "no gap long enough to be the wait between passes: {gaps:?}"
+    );
+}