Added the "mass-change" setter needed to modify a spacecraft mid-flight

tests/area_change.rs

@@ -0,0 +1,185 @@
+//! A spacecraft whose cross-section changes part-way through an arc.
+//!
+//! The case this exists for is a deployment or a jettison: the vehicle, not its state,
+//! is different after some instant, and drag and radiation pressure see a different
+//! area from then on. Nothing in the crate does this for you — it is a
+//! [`Callback`] that swaps the vehicle through
+//! [`Assembly::context_mut`](astrokit::propagation::Assembly::context_mut).
+//!
+//! Written as an integration test rather than a unit one, and importing from `astrokit`
+//! alone, so it compiles against exactly what a consumer can reach. A unit test would
+//! have the private surface available and would prove nothing about reachability.
+
+use std::path::Path;
+
+use astrokit::kernels;
+use astrokit::prelude::*;
+use astrokit::propagation::Assembly;
+use nalgebra::Vector3;
+
+/// Swaps the vehicle at a fixed time after the epoch.
+struct Deploy {
+    at_seconds: f64,
+    deployed: Spacecraft,
+}
+
+impl Callback<f64, 6, Assembly<f64, 6, Orbit>> for Deploy {
+    /// Zero at the moment of the change, so the integrator lands on it rather than
+    /// straddling it with whatever step it happened to be taking.
+    fn event(&self, t: f64, _y: &nalgebra::SVector<f64, 6>, _p: &Assembly<f64, 6, Orbit>) -> f64 {
+        t - self.at_seconds
+    }
+
+    fn effect(&mut self, ode: &mut ODE<'_, f64, 6, Assembly<f64, 6, Orbit>>) {
+        ode.params
+            .context_mut()
+            .set_spacecraft(self.deployed.clone());
+    }
+}
+
+const ARC: f64 = 86400.0;
+const MASS: f64 = 500.0;
+const SMALL: f64 = 5.0;
+const LARGE: f64 = 20.0;
+
+fn initial() -> State<f64, 6, Orbit> {
+    Orbit
+        .with_position(Vector3::new(6.778e6, 0.0, 0.0))
+        .with_velocity(Vector3::new(0.0, 7.668e3, 0.0))
+}
+
+/// A day of drag at a fixed area, returning where it ends up.
+fn fixed_area(area: f64, almanac: &std::sync::Arc<anise::almanac::Almanac>) -> Vector3<f64> {
+    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
+    let forces = [
+        TwoBody { body: Body::EARTH },
+        Drag {
+            density: DensityModel::harris_priester(),
+        },
+    ];
+    let mut propagation = Propagation::new(epoch, initial(), &forces)
+        .frame(Frame::EARTH_J2000)
+        .spacecraft(Spacecraft::cannonball(MASS, area, 2.2, 1.3))
+        .almanac(almanac.clone())
+        .duration(Duration::from_seconds(ARC))
+        .accuracy(Accuracy::High);
+    Orbit
+        .state(*propagation.run().states.last().unwrap())
+        .position()
+}
+
+/// The same day, small until midway and large after.
+fn deploying(almanac: &std::sync::Arc<anise::almanac::Almanac>) -> Vector3<f64> {
+    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
+    let forces = [
+        TwoBody { body: Body::EARTH },
+        Drag {
+            density: DensityModel::harris_priester(),
+        },
+    ];
+    let mut deploy = Deploy {
+        at_seconds: ARC / 2.0,
+        deployed: Spacecraft::cannonball(MASS, LARGE, 2.2, 1.3),
+    };
+    let mut propagation = Propagation::new(epoch, initial(), &forces)
+        .frame(Frame::EARTH_J2000)
+        .spacecraft(Spacecraft::cannonball(MASS, SMALL, 2.2, 1.3))
+        .almanac(almanac.clone())
+        .duration(Duration::from_seconds(ARC))
+        .accuracy(Accuracy::High)
+        .callback(&mut deploy);
+    Orbit
+        .state(*propagation.run().states.last().unwrap())
+        .position()
+}
+
+/// Deploying halfway must land strictly between never deploying and always deployed.
+///
+/// That is the assertion worth making: it can only hold if the change happened, and
+/// happened in the middle rather than at either end. Comparing against a remembered
+/// displacement would pass just as well with the callback firing at the wrong time, or
+/// with the area never changing at all.
+#[test]
+fn test_a_mid_arc_area_change_lands_between_the_two_fixed_areas() {
+    let directory = kernels::directory()
+        .unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
+    let Ok(almanac) = kernels::load(&directory).map(std::sync::Arc::new) else {
+        eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+        return;
+    };
+
+    let drag_free = {
+        let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
+        let mut plain = Propagation::new(epoch, initial(), &[TwoBody { body: Body::EARTH }])
+            .frame(Frame::EARTH_J2000)
+            .duration(Duration::from_seconds(ARC))
+            .accuracy(Accuracy::High);
+        Orbit.state(*plain.run().states.last().unwrap()).position()
+    };
+
+    let moved = |end: Vector3<f64>| (end - drag_free).norm();
+    let small = moved(fixed_area(SMALL, &almanac));
+    let large = moved(fixed_area(LARGE, &almanac));
+    let changed = moved(deploying(&almanac));
+
+    // Four times the area is four times the drag, so the bracket is wide and the
+    // ordering cannot be an accident of the integrator.
+    assert!(
+        large > small * 3.0,
+        "the two fixed areas should bracket widely: {small:.1} m and {large:.1} m"
+    );
+    assert!(
+        changed > small && changed < large,
+        "deploying halfway moved {changed:.1} m, outside the bracket \
+         [{small:.1}, {large:.1}] m set by never and always"
+    );
+
+    // Where between the two it lands says *when* the change happened, and the expected
+    // value is derivable rather than remembered. Drag removes energy, which shows up as
+    // along-track position error growing with the time left in the arc, so a bit of
+    // drag at time t contributes about (T - t). The extra area acts only over the
+    // second half, so its share is
+    //
+    //     int_{T/2}^{T} (T - t) dt   /   int_{0}^{T} (T - t) dt   =   1/4
+    //
+    // Measured 0.255, against the quarter that predicts. A change at the midpoint is
+    // the only thing that puts it there: firing at the start would give 1 and near the
+    // end would give 0, and both are far outside this band.
+    let fraction = (changed - small) / (large - small);
+    assert!(
+        (0.15..0.40).contains(&fraction),
+        "deploying halfway put it {fraction:.3} of the way from never to always, where \
+         the time-remaining weighting predicts about 0.25"
+    );
+}
+
+/// The vehicle the models see is the one the callback left, not the one it started with.
+///
+/// Separate from the trajectory assertion above because it checks the mechanism rather
+/// than its consequence: a swap that silently failed would still produce a trajectory,
+/// just the wrong one, and this says which.
+#[test]
+fn test_the_swap_reaches_the_models() {
+    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
+    let mut context = astrokit::dynamics::Context::new(epoch, Frame::EARTH_J2000);
+    context.set_spacecraft(Spacecraft::cannonball(MASS, SMALL, 2.2, 1.3));
+    let area = |surface: &SurfaceModel| match surface {
+        SurfaceModel::Cannonball { area, .. } => *area,
+        _ => panic!("a cannonball was attached"),
+    };
+    let before = area(&context.spacecraft().expect("attached").surface);
+
+    context.set_spacecraft(Spacecraft::cannonball(MASS, LARGE, 2.2, 1.3));
+    let after = area(&context.spacecraft().expect("attached").surface);
+
+    assert_eq!(before, SMALL);
+    assert_eq!(after, LARGE);
+
+    // The derived inverse inertia has to follow the vehicle, which is the whole reason
+    // this is a setter rather than a mutable reference to the field.
+    let heavy = Spacecraft::rectangular_box(MASS, Vector3::new(2.0, 1.0, 3.0), 2.2, 1.3);
+    let expected = heavy.inertia.try_inverse().expect("invertible");
+    context.set_spacecraft(heavy);
+    let (_, inverse) = context.inertia().expect("attached and invertible");
+    assert_eq!(*inverse, expected);
+}