Skip to main content

sidereon_core/astro/forces/
two_body.rs

1use crate::astro::constants::MU_EARTH;
2use crate::astro::error::PropagationError;
3use crate::astro::forces::r#trait::ForceModel;
4use crate::astro::propagator::api::PropagationContext;
5use crate::astro::state::CartesianState;
6use nalgebra::Vector3;
7
8pub struct TwoBodyGravity {
9    pub mu: f64,
10}
11
12impl Default for TwoBodyGravity {
13    fn default() -> Self {
14        Self { mu: MU_EARTH }
15    }
16}
17
18impl ForceModel for TwoBodyGravity {
19    fn acceleration(
20        &self,
21        state: &CartesianState,
22        _ctx: &PropagationContext,
23    ) -> Result<Vector3<f64>, PropagationError> {
24        let r_mag2 = state.position_km.norm_squared();
25        if r_mag2 == 0.0 {
26            return Err(PropagationError::NumericalFailure(
27                "Zero position magnitude".to_string(),
28            ));
29        }
30        let r_mag = r_mag2.sqrt();
31        Ok(state.position_km * (-self.mu / (r_mag2 * r_mag)))
32    }
33}
34
35#[cfg(test)]
36mod tests {
37    use super::*;
38    use crate::astro::propagator::api::PropagationContext;
39    use crate::astro::state::CartesianState;
40
41    #[test]
42    fn acceleration_matches_orbis_force_wrapper_bits() {
43        let state = CartesianState::new(0.0, [7000.0, -1210.0, 1300.0], [0.0, 0.0, 0.0]);
44        let acceleration = TwoBodyGravity::default()
45            .acceleration(&state, &PropagationContext::default())
46            .unwrap();
47
48        assert_eq!(acceleration.x.to_bits(), 13_798_562_943_973_640_097);
49        assert_eq!(acceleration.y.to_bits(), 4_563_548_234_789_153_053);
50        assert_eq!(acceleration.z.to_bits(), 13_787_359_517_156_423_902);
51    }
52}