sidereon_core/astro/propagator/
driver.rs

1//! High-level numerical-propagation driver.
2//!
3//! A single callable entry that turns a high-level configuration (force-model
4//! choice, integrator, integrator options, initial state, and an output epoch
5//! grid) into a sampled trajectory. It composes the existing
6//! [`StatePropagator`] / [`ForceModelKind`] / [`IntegratorOptions`] primitives
7//! with the canonical Earth constants from [`crate::astro::constants`]; it adds
8//! no integration math of its own. A language binding reduces to: marshal user
9//! options into a [`PropagationConfig`], call [`propagate_states`], then marshal
10//! the returned states back out, with no force-model or integrator policy left
11//! in the binding.
12//!
13//! The defaults come from [`IntegratorOptions::default`] and
14//! [`crate::astro::constants::MU_EARTH`], so the driver's defaults are exactly
15//! the engine's defaults rather than independently chosen numbers.
16
17use crate::astro::constants::{J2_EARTH, MU_EARTH, RE_EARTH};
18use crate::astro::error::PropagationError;
19use crate::astro::propagator::api::IntegratorOptions;
20use crate::astro::propagator::numerical::{ForceModelKind, IntegratorKind, StatePropagator};
21use crate::astro::state::CartesianState;
22
23/// High-level force-model choice for [`PropagationConfig`].
24///
25/// Mirrors the selector every language binding exposes: a point-mass two-body
26/// central force, or two-body plus the Earth J2 oblateness perturbation. The
27/// concrete [`ForceModelKind`] is composed by the driver from this choice and
28/// the configured gravitational parameter, filling the canonical Earth
29/// reference radius and J2 coefficient from [`crate::astro::constants`].
30#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
31pub enum PropagationForceModel {
32    /// Point-mass two-body gravity.
33    #[default]
34    TwoBody,
35    /// Two-body gravity plus the Earth J2 oblateness perturbation.
36    TwoBodyJ2,
37}
38
39/// High-level configuration for the numerical-propagation driver.
40///
41/// Construct with [`PropagationConfig::new`] and override fields; the defaults
42/// match every binding's defaults: two-body gravity, the adaptive DP54
43/// integrator, the canonical [`MU_EARTH`] gravitational parameter, and the
44/// engine-default [`IntegratorOptions`].
45#[derive(Debug, Clone, Copy, PartialEq)]
46pub struct PropagationConfig {
47    /// Initial ECI Cartesian state; its epoch is the propagation start epoch.
48    pub initial: CartesianState,
49    /// Force-model choice.
50    pub force_model: PropagationForceModel,
51    /// Gravitational-parameter override, km^3/s^2. `None` uses the canonical
52    /// [`MU_EARTH`].
53    pub mu_km3_s2: Option<f64>,
54    /// Integrator choice.
55    pub integrator: IntegratorKind,
56    /// Step-size / tolerance controls forwarded to the integrator.
57    pub options: IntegratorOptions,
58}
59
60impl PropagationConfig {
61    /// Build a config from a raw epoch (TDB seconds), ECI position (km), and ECI
62    /// velocity (km/s) with the binding defaults: two-body gravity, the DP54
63    /// integrator, the canonical [`MU_EARTH`] gravitational parameter, and
64    /// [`IntegratorOptions::default`].
65    pub fn new(epoch_tdb_seconds: f64, position_km: [f64; 3], velocity_km_s: [f64; 3]) -> Self {
66        Self {
67            initial: CartesianState::new(epoch_tdb_seconds, position_km, velocity_km_s),
68            force_model: PropagationForceModel::TwoBody,
69            mu_km3_s2: None,
70            integrator: IntegratorKind::Dp54,
71            options: IntegratorOptions::default(),
72        }
73    }
74
75    /// The gravitational parameter the driver will use: the configured override,
76    /// or the canonical [`MU_EARTH`] when none is set.
77    pub fn gravitational_parameter(&self) -> f64 {
78        self.mu_km3_s2.unwrap_or(MU_EARTH)
79    }
80
81    /// Compose the concrete [`ForceModelKind`] from the high-level choice and the
82    /// effective gravitational parameter, filling [`RE_EARTH`] / [`J2_EARTH`] for
83    /// the J2 variant. This is the force-model composition policy the bindings
84    /// each duplicated.
85    pub fn force_model_kind(&self) -> ForceModelKind {
86        let mu_km3_s2 = self.gravitational_parameter();
87        match self.force_model {
88            PropagationForceModel::TwoBody => ForceModelKind::TwoBody { mu_km3_s2 },
89            PropagationForceModel::TwoBodyJ2 => ForceModelKind::TwoBodyJ2 {
90                mu_km3_s2,
91                re_km: RE_EARTH,
92                j2: J2_EARTH,
93            },
94        }
95    }
96
97    /// Assemble the [`StatePropagator`] this config describes. Equivalent to the
98    /// per-binding hand assembly, so a downstream `ephemeris` call is bit-for-bit
99    /// identical to the binding's own.
100    pub fn to_propagator(&self) -> StatePropagator {
101        StatePropagator {
102            initial: self.initial,
103            force_model: self.force_model_kind(),
104            integrator: self.integrator,
105            options: self.options,
106        }
107    }
108}
109
110/// Run the numerical propagator described by `config`, sampling the trajectory
111/// at `epochs_tdb_seconds` (absolute TDB seconds, monotonic in the propagation
112/// direction).
113///
114/// Composes [`PropagationConfig::to_propagator`] with the existing
115/// [`StatePropagator::ephemeris`] sampler; the returned states are exactly what
116/// the engine produces, so the driver is bit-for-bit identical to assembling the
117/// propagator and calling `ephemeris` by hand.
118pub fn propagate_states(
119    config: &PropagationConfig,
120    epochs_tdb_seconds: &[f64],
121) -> Result<Vec<CartesianState>, PropagationError> {
122    config.to_propagator().ephemeris(epochs_tdb_seconds)
123}
124
125#[cfg(test)]
126mod tests {
127    use super::*;
128
129    fn circular_state() -> (f64, [f64; 3], [f64; 3]) {
130        let r: f64 = 7000.0;
131        let v = (MU_EARTH / r).sqrt();
132        (0.0, [r, 0.0, 0.0], [0.0, v, 0.0])
133    }
134
135    fn assert_states_bit_for_bit(left: &[CartesianState], right: &[CartesianState]) {
136        assert_eq!(left.len(), right.len(), "state-count mismatch");
137        for (a, b) in left.iter().zip(right.iter()) {
138            assert_eq!(a.epoch_tdb_seconds.to_bits(), b.epoch_tdb_seconds.to_bits());
139            for axis in 0..3 {
140                assert_eq!(
141                    a.position_array()[axis].to_bits(),
142                    b.position_array()[axis].to_bits(),
143                    "position axis {axis}"
144                );
145                assert_eq!(
146                    a.velocity_array()[axis].to_bits(),
147                    b.velocity_array()[axis].to_bits(),
148                    "velocity axis {axis}"
149                );
150            }
151        }
152    }
153
154    #[test]
155    fn config_defaults_come_from_core_constants() {
156        let (epoch, pos, vel) = circular_state();
157        let config = PropagationConfig::new(epoch, pos, vel);
158
159        assert_eq!(config.force_model, PropagationForceModel::TwoBody);
160        assert_eq!(config.mu_km3_s2, None);
161        assert_eq!(
162            config.gravitational_parameter().to_bits(),
163            MU_EARTH.to_bits()
164        );
165        assert_eq!(config.integrator, IntegratorKind::Dp54);
166        assert_eq!(config.options, IntegratorOptions::default());
167        assert_eq!(
168            config.force_model_kind(),
169            ForceModelKind::TwoBody {
170                mu_km3_s2: MU_EARTH
171            }
172        );
173    }
174
175    #[test]
176    fn force_model_kind_composes_j2_from_canonical_constants() {
177        let (epoch, pos, vel) = circular_state();
178        let config = PropagationConfig {
179            force_model: PropagationForceModel::TwoBodyJ2,
180            ..PropagationConfig::new(epoch, pos, vel)
181        };
182
183        assert_eq!(
184            config.force_model_kind(),
185            ForceModelKind::TwoBodyJ2 {
186                mu_km3_s2: MU_EARTH,
187                re_km: RE_EARTH,
188                j2: J2_EARTH,
189            }
190        );
191    }
192
193    #[test]
194    fn driver_two_body_default_matches_manual_composition_bit_for_bit() {
195        // The driver path: high-level config -> propagate_states.
196        let (epoch, pos, vel) = circular_state();
197        let config = PropagationConfig::new(epoch, pos, vel);
198        let epochs = [0.0, 600.0, 1800.0, 3600.0];
199        let via_driver = propagate_states(&config, &epochs).expect("driver propagation");
200
201        // The hand-assembled path the Python / WASM / C bindings each spell out:
202        // default mu = MU_EARTH, two-body force, DP54, default options.
203        let via_manual = StatePropagator {
204            initial: CartesianState::new(epoch, pos, vel),
205            force_model: ForceModelKind::TwoBody {
206                mu_km3_s2: MU_EARTH,
207            },
208            integrator: IntegratorKind::Dp54,
209            options: IntegratorOptions {
210                abs_tol: 1.0e-9,
211                rel_tol: 1.0e-12,
212                initial_step: 60.0,
213                min_step: 1.0e-6,
214                max_step: 3600.0,
215                max_steps: 1_000_000,
216                dense_output: false,
217            },
218        }
219        .ephemeris(&epochs)
220        .expect("manual propagation");
221
222        assert_states_bit_for_bit(&via_driver, &via_manual);
223    }
224
225    #[test]
226    fn driver_two_body_j2_custom_mu_rk4_matches_manual_composition_bit_for_bit() {
227        let (epoch, pos, vel) = circular_state();
228        let mu = 398_600.5;
229        let options = IntegratorOptions {
230            abs_tol: 1.0e-11,
231            rel_tol: 1.0e-13,
232            initial_step: 30.0,
233            min_step: 1.0e-5,
234            max_step: 120.0,
235            max_steps: 500_000,
236            dense_output: false,
237        };
238        let config = PropagationConfig {
239            force_model: PropagationForceModel::TwoBodyJ2,
240            mu_km3_s2: Some(mu),
241            integrator: IntegratorKind::Rk4,
242            options,
243            ..PropagationConfig::new(epoch, pos, vel)
244        };
245        let epochs = [0.0, 300.0, 900.0];
246        let via_driver = propagate_states(&config, &epochs).expect("driver propagation");
247
248        let via_manual = StatePropagator {
249            initial: CartesianState::new(epoch, pos, vel),
250            force_model: ForceModelKind::TwoBodyJ2 {
251                mu_km3_s2: mu,
252                re_km: RE_EARTH,
253                j2: J2_EARTH,
254            },
255            integrator: IntegratorKind::Rk4,
256            options,
257        }
258        .ephemeris(&epochs)
259        .expect("manual propagation");
260
261        assert_states_bit_for_bit(&via_driver, &via_manual);
262    }
263
264    #[test]
265    fn driver_surfaces_the_integrator_error_unchanged() {
266        // A non-positive initial step is rejected by the integrator itself; the
267        // driver forwards that error verbatim rather than masking it.
268        let (epoch, pos, vel) = circular_state();
269        let config = PropagationConfig {
270            options: IntegratorOptions {
271                initial_step: 0.0,
272                ..IntegratorOptions::default()
273            },
274            ..PropagationConfig::new(epoch, pos, vel)
275        };
276
277        let err = propagate_states(&config, &[0.0, 60.0]).expect_err("non-positive step rejected");
278        match err {
279            PropagationError::InvalidInput(message) => {
280                assert!(message.contains("initial_step"), "{message}");
281                assert!(message.contains("not positive"), "{message}");
282            }
283            other => panic!("expected invalid-input error, got {other:?}"),
284        }
285    }
286}
sidereon_core/astro/propagator/driver.rs

sidereon_core/astro/propagator/
driver.rs
