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

sidereon_core/astro/propagator/
driver.rs
