sidereon_core/astro/forces/
srp.rs

1//! Cannonball solar radiation pressure.
2//!
3//! The acceleration follows the Montenbruck & Gill cannonball SRP model,
4//! `-P_sr * (AU / |r_sun|)^2 * C_R * (A/m) * nu * s_hat`, where `s_hat`
5//! points from the satellite to the Sun. SI pressure and area-to-mass units are
6//! used internally, then converted from m/s^2 to km/s^2 for the propagation
7//! state.
8
9use crate::astro::bodies::sun_moon::sun_moon_eci;
10use crate::astro::constants::astro::{AU_KM, SOLAR_RADIATION_PRESSURE_N_M2};
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::events::eclipse::shadow_fraction;
15use crate::astro::forces::r#trait::ForceModel;
16use crate::astro::propagator::api::PropagationContext;
17use crate::astro::state::CartesianState;
18use nalgebra::Vector3;
19
20/// Cannonball solar radiation pressure parameters.
21#[derive(Debug, Clone, Copy, PartialEq)]
22pub struct SolarRadiationPressure {
23    cr: f64,
24    area_to_mass_m2_kg: f64,
25    pressure_n_m2: f64,
26    au_km: f64,
27}
28
29impl SolarRadiationPressure {
30    /// Build from reflectivity coefficient `C_R` and area-to-mass `A/m` in m^2/kg.
31    pub fn new(cr: f64, area_to_mass_m2_kg: f64) -> Result<Self, PropagationError> {
32        Self::with_pressure(cr, area_to_mass_m2_kg, SOLAR_RADIATION_PRESSURE_N_M2, AU_KM)
33    }
34
35    /// Build with explicit pressure and AU constants.
36    pub fn with_pressure(
37        cr: f64,
38        area_to_mass_m2_kg: f64,
39        pressure_n_m2: f64,
40        au_km: f64,
41    ) -> Result<Self, PropagationError> {
42        validate_positive("cr", cr)?;
43        validate_positive("area_to_mass_m2_kg", area_to_mass_m2_kg)?;
44        validate_positive("pressure_n_m2", pressure_n_m2)?;
45        validate_positive("au_km", au_km)?;
46        Ok(Self {
47            cr,
48            area_to_mass_m2_kg,
49            pressure_n_m2,
50            au_km,
51        })
52    }
53
54    /// Reflectivity coefficient `C_R`.
55    pub fn cr(&self) -> f64 {
56        self.cr
57    }
58
59    /// Area-to-mass ratio `A/m`, m^2/kg.
60    pub fn area_to_mass_m2_kg(&self) -> f64 {
61        self.area_to_mass_m2_kg
62    }
63
64    /// Solar radiation pressure at 1 AU, N/m^2.
65    pub fn pressure_n_m2(&self) -> f64 {
66        self.pressure_n_m2
67    }
68
69    /// Astronomical unit used in the inverse-square scale, km.
70    pub fn au_km(&self) -> f64 {
71        self.au_km
72    }
73}
74
75impl ForceModel for SolarRadiationPressure {
76    fn acceleration(
77        &self,
78        state: &CartesianState,
79        _ctx: &PropagationContext,
80    ) -> Result<Vector3<f64>, PropagationError> {
81        let sun_pos_km = sun_position_km(state.epoch_tdb_seconds)?;
82        srp_acceleration_with_sun(state.position_km, sun_pos_km, *self)
83    }
84}
85
86fn sun_position_km(epoch_tdb_seconds: f64) -> Result<Vector3<f64>, PropagationError> {
87    let t_tt_centuries = epoch_tdb_seconds / (DAYS_PER_JULIAN_CENTURY * SECONDS_PER_DAY);
88    let bodies = sun_moon_eci(t_tt_centuries)
89        .map_err(|error| PropagationError::ForceModelFailure(format!("Sun/Moon: {error}")))?;
90    Ok(Vector3::new(
91        bodies.sun[0] / M_PER_KM,
92        bodies.sun[1] / M_PER_KM,
93        bodies.sun[2] / M_PER_KM,
94    ))
95}
96
97fn srp_acceleration_with_sun(
98    sat_pos_km: Vector3<f64>,
99    sun_pos_km: Vector3<f64>,
100    model: SolarRadiationPressure,
101) -> Result<Vector3<f64>, PropagationError> {
102    let sat_to_sun = sun_pos_km - sat_pos_km;
103    let sat_to_sun_r2 = sat_to_sun.norm_squared();
104    if sat_to_sun_r2 == 0.0 {
105        return Err(PropagationError::NumericalFailure(
106            "Satellite coincident with Sun".to_string(),
107        ));
108    }
109    let sun_r2 = sun_pos_km.norm_squared();
110    if sun_r2 == 0.0 {
111        return Err(PropagationError::NumericalFailure(
112            "Zero Sun position magnitude".to_string(),
113        ));
114    }
115
116    let shadow = shadow_fraction(
117        [sat_pos_km.x, sat_pos_km.y, sat_pos_km.z],
118        [sun_pos_km.x, sun_pos_km.y, sun_pos_km.z],
119    )
120    .map_err(|error| PropagationError::ForceModelFailure(format!("eclipse: {error}")))?;
121    let illumination = 1.0 - shadow;
122    if illumination == 0.0 {
123        return Ok(Vector3::zeros());
124    }
125
126    let s_hat = sat_to_sun / sat_to_sun_r2.sqrt();
127    let pressure_scale = model.pressure_n_m2 * (model.au_km * model.au_km / sun_r2);
128    let accel_m_s2 = -pressure_scale * model.cr * model.area_to_mass_m2_kg * illumination * s_hat;
129    Ok(accel_m_s2 / M_PER_KM)
130}
131
132fn validate_positive(field: &'static str, value: f64) -> Result<(), PropagationError> {
133    if !value.is_finite() || value <= 0.0 {
134        return Err(PropagationError::InvalidInput(format!(
135            "{field} must be finite and positive"
136        )));
137    }
138    Ok(())
139}
140
141#[cfg(test)]
142mod tests {
143    //! Clean-room validation anchors:
144    //!
145    //! The expected 1 AU acceleration magnitude is the published cannonball SRP
146    //! pressure `4.56e-6 N/m^2` multiplied by `C_R * A/m` and converted from
147    //! m/s^2 to km/s^2. The shadow tests use fixed conical geometry and the
148    //! existing eclipse helper, whose return value is geometric shadow fraction.
149
150    use super::*;
151    use crate::astro::constants::astro::AU_KM;
152
153    #[test]
154    fn one_au_magnitude_matches_pressure_area_mass_scale() {
155        let model = SolarRadiationPressure::new(1.2, 0.02).expect("valid SRP");
156        let sat = Vector3::new(7000.0, 0.0, 0.0);
157        let sun = Vector3::new(AU_KM, 0.0, 0.0);
158        let accel = srp_acceleration_with_sun(sat, sun, model).expect("SRP acceleration");
159
160        let expected_km_s2 = 4.56e-6 * 1.2 * 0.02 / M_PER_KM;
161        assert!((accel.norm() - expected_km_s2).abs() < 1.0e-24);
162    }
163
164    #[test]
165    fn acceleration_is_zero_in_deep_umbra() {
166        let model = SolarRadiationPressure::new(1.2, 0.02).expect("valid SRP");
167        let sat = Vector3::new(-7000.0, 0.0, 0.0);
168        let sun = Vector3::new(AU_KM, 0.0, 0.0);
169        let accel = srp_acceleration_with_sun(sat, sun, model).expect("SRP acceleration");
170
171        assert_eq!(accel, Vector3::zeros());
172    }
173
174    #[test]
175    fn direction_is_away_from_sun() {
176        let model = SolarRadiationPressure::new(1.0, 0.01).expect("valid SRP");
177        let sat = Vector3::new(7000.0, 0.0, 0.0);
178        let sun = Vector3::new(AU_KM, 0.0, 0.0);
179        let accel = srp_acceleration_with_sun(sat, sun, model).expect("SRP acceleration");
180
181        assert!(accel.x < 0.0);
182        assert_eq!(accel.y, 0.0);
183        assert_eq!(accel.z, 0.0);
184    }
185}
sidereon_core/astro/forces/srp.rs

sidereon_core/astro/forces/
srp.rs
