sidereon_core/astro/forces/
relativity.rs

1//! Geocentric Schwarzschild relativistic correction.
2//!
3//! The acceleration is the single central-body term
4//! `GM/(c^2 r^3) * [ (4 GM/r - v^2) r + 4 (r . v) v ]`, as given by
5//! Montenbruck & Gill, "Satellite Orbits", equation 3.146. Position is in km,
6//! velocity is in km/s, and the returned acceleration is km/s^2.
7
8use crate::astro::constants::physics::SPEED_OF_LIGHT_KM_S;
9use crate::astro::constants::MU_EARTH;
10use crate::astro::error::PropagationError;
11use crate::astro::forces::r#trait::ForceModel;
12use crate::astro::propagator::api::PropagationContext;
13use crate::astro::state::CartesianState;
14use nalgebra::Vector3;
15
16/// Schwarzschild correction for a single central body.
17#[derive(Debug, Clone, Copy, PartialEq)]
18pub struct SchwarzschildRelativity {
19    /// Central-body gravitational parameter, km^3/s^2.
20    pub mu_km3_s2: f64,
21    /// Speed of light, km/s.
22    pub c_km_s: f64,
23}
24
25impl Default for SchwarzschildRelativity {
26    fn default() -> Self {
27        Self {
28            mu_km3_s2: MU_EARTH,
29            c_km_s: SPEED_OF_LIGHT_KM_S,
30        }
31    }
32}
33
34impl SchwarzschildRelativity {
35    /// Build a Schwarzschild correction with explicit constants.
36    pub fn new(mu_km3_s2: f64, c_km_s: f64) -> Self {
37        Self { mu_km3_s2, c_km_s }
38    }
39}
40
41impl ForceModel for SchwarzschildRelativity {
42    fn acceleration(
43        &self,
44        state: &CartesianState,
45        _ctx: &PropagationContext,
46    ) -> Result<Vector3<f64>, PropagationError> {
47        if self.c_km_s.is_infinite() && self.c_km_s.is_sign_positive() {
48            return Ok(Vector3::zeros());
49        }
50        if !self.c_km_s.is_finite() || self.c_km_s <= 0.0 {
51            return Err(PropagationError::InvalidInput(
52                "c_km_s must be finite and positive".to_string(),
53            ));
54        }
55
56        let r2 = state.position_km.norm_squared();
57        if r2 == 0.0 {
58            return Err(PropagationError::NumericalFailure(
59                "Zero position magnitude".to_string(),
60            ));
61        }
62        let r = r2.sqrt();
63        let r3 = r2 * r;
64        let v2 = state.velocity_km_s.norm_squared();
65        let r_dot_v = state.position_km.dot(&state.velocity_km_s);
66        let factor = self.mu_km3_s2 / (self.c_km_s * self.c_km_s * r3);
67        let bracket = state.position_km * (4.0 * self.mu_km3_s2 / r - v2)
68            + state.velocity_km_s * (4.0 * r_dot_v);
69
70        Ok(bracket * factor)
71    }
72}
73
74#[cfg(test)]
75mod tests {
76    //! Clean-room validation anchors:
77    //!
78    //! For a circular GNSS-altitude orbit, `r . v = 0` and `v^2 = GM/r`, so
79    //! Montenbruck & Gill equation 3.146 reduces to magnitude
80    //! `3 * GM^2 / (c^2 * r^3)`. At `r = 26560 km`, this is
81    //! `2.830552329290399e-10 m/s^2`.
82
83    use super::*;
84
85    #[test]
86    fn gnss_circular_magnitude_matches_closed_form() {
87        let r = 26_560.0;
88        let v = (MU_EARTH / r).sqrt();
89        let state = CartesianState::new(0.0, [r, 0.0, 0.0], [0.0, v, 0.0]);
90        let accel = SchwarzschildRelativity::default()
91            .acceleration(&state, &PropagationContext::default())
92            .expect("relativistic acceleration");
93        let expected = 2.830_552_329_290_399e-13;
94
95        assert!((accel.x - expected).abs() < 1.0e-27);
96        assert_eq!(accel.y, 0.0);
97        assert_eq!(accel.z, 0.0);
98        assert!((1.0e-10..1.0e-9).contains(&(accel.norm() * 1000.0)));
99    }
100
101    #[test]
102    fn acceleration_vanishes_in_newtonian_limit() {
103        let state = CartesianState::new(0.0, [26_560.0, 0.0, 0.0], [0.0, 3.874, 0.0]);
104        let accel = SchwarzschildRelativity::new(MU_EARTH, f64::INFINITY)
105            .acceleration(&state, &PropagationContext::default())
106            .expect("relativistic acceleration");
107        assert_eq!(accel, Vector3::zeros());
108    }
109}
sidereon_core/astro/forces/relativity.rs

sidereon_core/astro/forces/
relativity.rs
