sidereon_core/astro/forces/
third_body.rs

1//! Sun and Moon third-body point-mass perturbations.
2//!
3//! The acceleration uses the standard direct-minus-indirect expression
4//! `mu_b * ((r_b - r) / |r_b - r|^3 - r_b / |r_b|^3)`, as given by
5//! Montenbruck & Gill, "Satellite Orbits", section 3.2, and Vallado,
6//! "Fundamentals of Astrodynamics and Applications", section 8.6. Body
7//! positions come from the crate's analytic Sun/Moon series at the state epoch.
8
9use crate::astro::bodies::sun_moon::sun_moon_eci;
10use crate::astro::constants::astro::{GM_MOON_KM3_S2, GM_SUN_KM3_S2};
11use crate::astro::constants::time::{DAYS_PER_JULIAN_CENTURY, SECONDS_PER_DAY};
12use crate::astro::constants::units::M_PER_KM;
13use crate::astro::error::PropagationError;
14use crate::astro::forces::r#trait::ForceModel;
15use crate::astro::propagator::api::PropagationContext;
16use crate::astro::state::CartesianState;
17use nalgebra::Vector3;
18
19/// Enabled third bodies for [`ThirdBodyGravity`].
20#[derive(Debug, Clone, Copy, PartialEq, Eq)]
21pub struct ThirdBodyBodies {
22    /// Include the Sun perturbation.
23    pub sun: bool,
24    /// Include the Moon perturbation.
25    pub moon: bool,
26}
27
28impl ThirdBodyBodies {
29    /// No third bodies.
30    pub const NONE: Self = Self {
31        sun: false,
32        moon: false,
33    };
34
35    /// Sun only.
36    pub const SUN: Self = Self {
37        sun: true,
38        moon: false,
39    };
40
41    /// Moon only.
42    pub const MOON: Self = Self {
43        sun: false,
44        moon: true,
45    };
46
47    /// Sun and Moon.
48    pub const SUN_AND_MOON: Self = Self {
49        sun: true,
50        moon: true,
51    };
52}
53
54/// Third-body gravity from analytic Sun and Moon positions.
55#[derive(Debug, Clone, Copy, PartialEq)]
56pub struct ThirdBodyGravity {
57    /// Enabled bodies.
58    pub bodies: ThirdBodyBodies,
59    /// Solar gravitational parameter, km^3/s^2.
60    pub gm_sun_km3_s2: f64,
61    /// Lunar gravitational parameter, km^3/s^2.
62    pub gm_moon_km3_s2: f64,
63}
64
65impl Default for ThirdBodyGravity {
66    fn default() -> Self {
67        Self {
68            bodies: ThirdBodyBodies::SUN_AND_MOON,
69            gm_sun_km3_s2: GM_SUN_KM3_S2,
70            gm_moon_km3_s2: GM_MOON_KM3_S2,
71        }
72    }
73}
74
75impl ThirdBodyGravity {
76    /// Build a third-body model with explicit bodies and gravitational parameters.
77    pub fn new(bodies: ThirdBodyBodies, gm_sun_km3_s2: f64, gm_moon_km3_s2: f64) -> Self {
78        Self {
79            bodies,
80            gm_sun_km3_s2,
81            gm_moon_km3_s2,
82        }
83    }
84
85    /// Sun-only third-body gravity.
86    pub fn sun() -> Self {
87        Self {
88            bodies: ThirdBodyBodies::SUN,
89            ..Self::default()
90        }
91    }
92
93    /// Moon-only third-body gravity.
94    pub fn moon() -> Self {
95        Self {
96            bodies: ThirdBodyBodies::MOON,
97            ..Self::default()
98        }
99    }
100}
101
102impl ForceModel for ThirdBodyGravity {
103    fn acceleration(
104        &self,
105        state: &CartesianState,
106        _ctx: &PropagationContext,
107    ) -> Result<Vector3<f64>, PropagationError> {
108        let bodies = sun_moon_positions_km(state.epoch_tdb_seconds)?;
109        let mut accel = Vector3::zeros();
110        if self.bodies.sun {
111            accel += third_body_acceleration(state.position_km, bodies.sun_km, self.gm_sun_km3_s2)?;
112        }
113        if self.bodies.moon {
114            accel +=
115                third_body_acceleration(state.position_km, bodies.moon_km, self.gm_moon_km3_s2)?;
116        }
117        Ok(accel)
118    }
119}
120
121struct SunMoonKm {
122    sun_km: Vector3<f64>,
123    moon_km: Vector3<f64>,
124}
125
126fn sun_moon_positions_km(epoch_tdb_seconds: f64) -> Result<SunMoonKm, PropagationError> {
127    let t_tt_centuries = epoch_tdb_seconds / (DAYS_PER_JULIAN_CENTURY * SECONDS_PER_DAY);
128    let bodies = sun_moon_eci(t_tt_centuries)
129        .map_err(|error| PropagationError::ForceModelFailure(format!("Sun/Moon: {error}")))?;
130    Ok(SunMoonKm {
131        sun_km: meters_to_km(bodies.sun),
132        moon_km: meters_to_km(bodies.moon),
133    })
134}
135
136fn meters_to_km(position_m: [f64; 3]) -> Vector3<f64> {
137    Vector3::new(
138        position_m[0] / M_PER_KM,
139        position_m[1] / M_PER_KM,
140        position_m[2] / M_PER_KM,
141    )
142}
143
144fn third_body_acceleration(
145    sat_pos_km: Vector3<f64>,
146    body_pos_km: Vector3<f64>,
147    gm_body_km3_s2: f64,
148) -> Result<Vector3<f64>, PropagationError> {
149    let body_r2 = body_pos_km.norm_squared();
150    if body_r2 == 0.0 {
151        return Err(PropagationError::NumericalFailure(
152            "Zero third-body position magnitude".to_string(),
153        ));
154    }
155    let delta = body_pos_km - sat_pos_km;
156    let delta_r2 = delta.norm_squared();
157    if delta_r2 == 0.0 {
158        return Err(PropagationError::NumericalFailure(
159            "Satellite coincident with third body".to_string(),
160        ));
161    }
162
163    let body_r = body_r2.sqrt();
164    let delta_r = delta_r2.sqrt();
165    Ok(gm_body_km3_s2 * (delta / (delta_r2 * delta_r) - body_pos_km / (body_r2 * body_r)))
166}
167
168#[cfg(test)]
169mod tests {
170    //! Clean-room validation anchors:
171    //!
172    //! The exact vector checks use fixed Sun and Moon stand-ins and numeric
173    //! values evaluated independently from the published direct-minus-indirect
174    //! third-body equation in Montenbruck & Gill, section 3.2. The magnitude
175    //! checks are the textbook LEO order of `1e-7..1e-6 m/s^2` for solar and
176    //! lunar third-body perturbations.
177
178    use super::*;
179    use crate::astro::constants::astro::AU_KM;
180
181    fn assert_close_vector(actual: Vector3<f64>, expected: [f64; 3], tolerance: f64) {
182        for axis in 0..3 {
183            assert!(
184                (actual[axis] - expected[axis]).abs() <= tolerance,
185                "axis {axis}: actual={} expected={}",
186                actual[axis],
187                expected[axis]
188            );
189        }
190    }
191
192    #[test]
193    fn fixed_sun_and_moon_vectors_match_independent_closed_form() {
194        let sat = Vector3::new(7000.0, -1210.0, 1300.0);
195        let sun = third_body_acceleration(sat, Vector3::new(AU_KM, 0.0, 0.0), GM_SUN_KM3_S2)
196            .expect("sun acceleration");
197        let moon = third_body_acceleration(sat, Vector3::new(384_400.0, 0.0, 0.0), GM_MOON_KM3_S2)
198            .expect("moon acceleration");
199
200        assert_close_vector(
201            sun,
202            [
203                5.549_999_393_729_882e-10,
204                4.797_132_723_126_184e-11,
205                -5.153_944_247_986_81e-11,
206            ],
207            1.0e-20,
208        );
209        assert_close_vector(
210            moon,
211            [
212                1.241_117_020_105_779_8e-9,
213                1.103_594_278_427_440_5e-10,
214                -1.185_679_803_269_151e-10,
215            ],
216            1.0e-20,
217        );
218    }
219
220    #[test]
221    fn indirect_term_cancels_at_geocenter() {
222        let state = CartesianState::new(0.0, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]);
223        let accel = ThirdBodyGravity::default()
224            .acceleration(&state, &PropagationContext::default())
225            .expect("third-body acceleration");
226        assert_eq!(accel, Vector3::zeros());
227    }
228
229    #[test]
230    fn leo_magnitudes_match_textbook_order() {
231        let state = CartesianState::new(0.0, [7000.0, 0.0, 0.0], [0.0, 7.5, 0.0]);
232        let sun = ThirdBodyGravity::sun()
233            .acceleration(&state, &PropagationContext::default())
234            .expect("sun acceleration")
235            .norm()
236            * M_PER_KM;
237        let moon = ThirdBodyGravity::moon()
238            .acceleration(&state, &PropagationContext::default())
239            .expect("moon acceleration")
240            .norm()
241            * M_PER_KM;
242
243        assert!(
244            (1.0e-7..1.0e-6).contains(&sun),
245            "sun perturbation magnitude {sun} m/s^2"
246        );
247        assert!(
248            (1.0e-7..2.0e-6).contains(&moon),
249            "moon perturbation magnitude {moon} m/s^2"
250        );
251    }
252}
sidereon_core/astro/forces/third_body.rs

sidereon_core/astro/forces/
third_body.rs
