Skip to main content

sidereon_core/astro/
mod.rs

1//! Numerical astrodynamics engine for orbit propagation, force models, and
2//! future flight-dynamics primitives.
3//!
4//! Current scope:
5//!
6//! - inertial Cartesian state representation
7//! - two-body gravity and J2 perturbation
8//! - fixed-step RK4
9//! - adaptive Dormand-Prince 5(4) (`DP54`)
10//! - propagation results with step statistics
11//!
12//! Planned future work includes dense output, event handling, richer propagation
13//! contexts, additional force models, covariance propagation, estimation, and
14//! maneuver support.
15
16pub mod angles;
17pub mod atmosphere;
18pub mod bodies;
19pub mod cdm;
20pub mod conjunction;
21pub mod constants;
22pub mod covariance;
23pub mod data;
24pub mod error;
25pub mod events;
26pub mod forces;
27pub mod frames;
28pub mod integrators;
29pub mod iod;
30pub mod lambert;
31pub mod math;
32pub mod omm;
33pub mod passes;
34pub mod propagator;
35pub mod rf;
36pub mod sgp4;
37pub mod spk;
38pub mod state;
39pub mod tca;
40pub mod time;
41pub mod tle;
42pub mod tolerances;
43pub mod xml;
44
45pub use spk::{
46    DafByteOrder, DafFileRecord, DafSpk, Spk, SpkError, SpkSegmentDescriptor, SpkState,
47    SpkStateVector,
48};
49
50#[cfg(all(feature = "sgp4-debug-oracle", sgp4_oracle_built))]
51#[doc(hidden)]
52pub mod sgp4_cpp_oracle {
53    //! Test-only oracle bridge to the Vallado C++ implementation.
54    //! Compiled in only when the `sgp4-debug-oracle` feature is on and the
55    //! development-only C++ oracle sources were found by the build script.
56    //! Not part of the public API.
57
58    use std::os::raw::{c_char, c_double, c_int};
59
60    pub const CPP_DUMP_DOUBLE_COUNT: usize = 112;
61    pub const CPP_DUMP_INT_COUNT: usize = 5;
62
63    extern "C" {
64        pub fn cpp_sgp4init_dump(
65            satnum: *const c_char,
66            epoch_sgp4: c_double,
67            bstar: c_double,
68            ndot: c_double,
69            nddot: c_double,
70            ecco: c_double,
71            argpo: c_double,
72            inclo: c_double,
73            mo: c_double,
74            no_kozai: c_double,
75            nodeo: c_double,
76            epochyr: c_int,
77            epochdays: c_double,
78            jdsatepoch: c_double,
79            jdsatepoch_frac: c_double,
80            double_out: *mut c_double,
81            int_out: *mut c_int,
82        ) -> c_int;
83
84        pub fn cpp_sgp4_step(
85            satnum: *const c_char,
86            epoch_sgp4: c_double,
87            bstar: c_double,
88            ndot: c_double,
89            nddot: c_double,
90            ecco: c_double,
91            argpo: c_double,
92            inclo: c_double,
93            mo: c_double,
94            no_kozai: c_double,
95            nodeo: c_double,
96            epochyr: c_int,
97            epochdays: c_double,
98            jdsatepoch: c_double,
99            jdsatepoch_frac: c_double,
100            tsince: c_double,
101            r_out: *mut c_double,
102            v_out: *mut c_double,
103        ) -> c_int;
104    }
105
106    /// Force-reference the C symbols so the linker pulls in the static lib.
107    /// Without this, the rlib has no use of the symbols and the linker
108    /// strips the entire archive when compiling integration tests.
109    #[doc(hidden)]
110    pub fn force_link_oracle() -> usize {
111        let init_dump = cpp_sgp4init_dump as *const ();
112        let step = cpp_sgp4_step as *const ();
113
114        init_dump as usize ^ step as usize
115    }
116}
117
118#[cfg(all(feature = "sgp4-debug-oracle", sgp4_oracle_built))]
119pub use sgp4_cpp_oracle::cpp_sgp4_step;
120
121pub use error::PropagationError;
122pub use state::CartesianState;
123pub use time::Time;
124
125#[cfg(test)]
126mod tests {
127    use super::*;
128    use crate::astro::forces::TwoBodyGravity;
129    use crate::astro::integrators::{Integrator, DP54};
130    use crate::astro::propagator::{api::IntegratorOptions, OrbitalDynamics, PropagationContext};
131    use nalgebra::Vector3;
132
133    #[test]
134    fn test_two_body_dp54_precision() {
135        let r_mag: f64 = 7000.0;
136        let mu: f64 = 398600.4418;
137        let v_mag: f64 = (mu / r_mag).sqrt();
138        let initial_state = CartesianState {
139            epoch_tdb_seconds: 0.0,
140            position_km: Vector3::new(r_mag, 0.0, 0.0),
141            velocity_km_s: Vector3::new(0.0, v_mag, 0.0),
142        };
143
144        let force = TwoBodyGravity::default();
145        let dynamics = OrbitalDynamics {
146            force_model: &force,
147        };
148        let integrator = DP54;
149        let ctx = PropagationContext::default();
150        let opts = IntegratorOptions {
151            abs_tol: 1e-12,
152            rel_tol: 1e-12,
153            initial_step: 1.0,
154            min_step: 1e-15,
155            ..IntegratorOptions::default()
156        };
157
158        let period = 2.0 * std::f64::consts::PI * (r_mag.powi(3) / mu).sqrt();
159        let result = integrator
160            .propagate(initial_state, period, &dynamics, &ctx, &opts)
161            .unwrap();
162
163        let final_pos = result.final_state.position_km;
164        let final_vel = result.final_state.velocity_km_s;
165
166        // Oracle 1: Return to start precision (Sub-millimeter)
167        assert!(
168            (final_pos.x - r_mag).abs() < 1e-7,
169            "Position X error too large: {}",
170            (final_pos.x - r_mag).abs()
171        );
172        assert!(
173            final_pos.y.abs() < 1e-7,
174            "Position Y error too large: {}",
175            final_pos.y.abs()
176        );
177
178        // Oracle 2: Energy conservation (Specific mechanical energy)
179        let initial_energy = v_mag.powi(2) / 2.0 - mu / r_mag;
180        let final_v_mag = final_vel.norm();
181        let final_r_mag = final_pos.norm();
182        let final_energy = final_v_mag.powi(2) / 2.0 - mu / final_r_mag;
183        assert!(
184            (final_energy - initial_energy).abs() < 1e-10,
185            "Energy conservation failure: {}",
186            (final_energy - initial_energy).abs()
187        );
188    }
189}