nyx_space/dynamics/
solarpressure.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use super::{DynamicsAlmanacSnafu, DynamicsError, DynamicsPlanetarySnafu, ForceModel};
20use crate::cosmic::eclipse::ShadowModel;
21use crate::cosmic::{AU, Frame, SPEED_OF_LIGHT_M_S, Spacecraft};
22use crate::linalg::{Const, Matrix4x3, Vector3};
23use anise::almanac::Almanac;
24use anise::constants::frames::{EARTH_J2000, SUN_J2000};
25use hyperdual::{Float, OHyperdual, hyperspace_from_vector, linalg::norm};
26use log::warn;
27use serde::{Deserialize, Serialize};
28use serde_dhall::StaticType;
29use snafu::ResultExt;
30use std::fmt;
31use std::sync::Arc;
32
33// Default solar flux in W/m^2
34#[allow(non_upper_case_globals)]
35pub const SOLAR_FLUX_W_m2: f64 = 1367.0;
36
37#[cfg(feature = "python")]
38use pyo3::prelude::*;
39
40/// Computation of solar radiation pressure is based on STK: <http://help.agi.com/stk/index.htm#gator/eq-solar.htm> .
41///
42/// :type shadow_bodies: list[Frame]
43/// :type almanac: Almanac
44/// :type flux_w_m2: float
45/// :type estimate: bool
46#[derive(Clone, Debug, Serialize, Deserialize, StaticType)]
47#[cfg_attr(feature = "python", pyclass(from_py_object, get_all, set_all))]
48pub struct SolarPressure {
49    /// solar flux at 1 AU, in W/m^2
50    pub phi: f64,
51    pub shadow_model: ShadowModel,
52    /// Set to true to estimate the coefficient of reflectivity
53    pub estimate: bool,
54}
55
56impl Default for SolarPressure {
57    /// Important: the default will FAIL at runtime if the shadow model is not manually defined with loaded frames.
58    fn default() -> Self {
59        Self {
60            phi: SOLAR_FLUX_W_m2,
61            estimate: false,
62            shadow_model: ShadowModel {
63                light_source: SUN_J2000,
64                shadow_bodies: vec![EARTH_J2000],
65            },
66        }
67    }
68}
69
70impl SolarPressure {
71    /// Will set the solar flux at 1 AU to: Phi = 1367.0
72    pub fn default_flux_raw(
73        shadow_bodies: Vec<Frame>,
74        almanac: &Almanac,
75    ) -> Result<Self, DynamicsError> {
76        let shadow_model = ShadowModel {
77            light_source: almanac.frame_info(SUN_J2000).context({
78                DynamicsPlanetarySnafu {
79                    action: "planetary data from third body not loaded",
80                }
81            })?,
82            shadow_bodies: shadow_bodies
83                .iter()
84                .filter_map(|object| match almanac.frame_info(object) {
85                    Ok(loaded_obj) => Some(loaded_obj),
86                    Err(e) => {
87                        warn!("when initializing SRP model for {object}, {e}");
88                        None
89                    }
90                })
91                .collect(),
92        };
93        Ok(Self {
94            phi: SOLAR_FLUX_W_m2,
95            shadow_model,
96            estimate: true,
97        })
98    }
99
100    /// Accounts for the shadowing of only one body and will set the solar flux at 1 AU to: Phi = 1367.0
101    pub fn default_flux(shadow_body: Frame, almanac: &Almanac) -> Result<Arc<Self>, DynamicsError> {
102        Ok(Arc::new(Self::default_flux_raw(
103            vec![shadow_body],
104            almanac,
105        )?))
106    }
107
108    /// Accounts for the shadowing of only one body and will set the solar flux at 1 AU to: Phi = 1367.0
109    pub fn default_no_estimation(
110        shadow_bodies: Vec<Frame>,
111        almanac: &Almanac,
112    ) -> Result<Arc<Self>, DynamicsError> {
113        let mut srp = Self::default_flux_raw(shadow_bodies, almanac)?;
114        srp.estimate = false;
115        Ok(Arc::new(srp))
116    }
117
118    /// Must provide the flux in W/m^2
119    pub fn with_flux(
120        flux_w_m2: f64,
121        shadow_bodies: Vec<Frame>,
122        almanac: &Almanac,
123    ) -> Result<Arc<Self>, DynamicsError> {
124        let mut me = Self::default_flux_raw(shadow_bodies, almanac)?;
125        me.phi = flux_w_m2;
126        Ok(Arc::new(me))
127    }
128
129    /// Solar radiation pressure force model accounting for the provided shadow bodies.
130    pub fn new(shadow_bodies: Vec<Frame>, almanac: &Almanac) -> Result<Arc<Self>, DynamicsError> {
131        Ok(Arc::new(Self::default_flux_raw(shadow_bodies, almanac)?))
132    }
133}
134
135impl ForceModel for SolarPressure {
136    fn estimation_index(&self) -> Option<usize> {
137        if self.estimate { Some(6) } else { None }
138    }
139
140    fn eom(&self, ctx: &Spacecraft, almanac: &Almanac) -> Result<Vector3<f64>, DynamicsError> {
141        let osc = ctx.orbit;
142        // Compute the position of the Sun as seen from the spacecraft
143        let r_sun = almanac
144            .transform_to(ctx.orbit, self.shadow_model.light_source, None)
145            .context(DynamicsAlmanacSnafu {
146                action: "transforming state to vector seen from Sun",
147            })?
148            .radius_km;
149
150        let r_sun_unit = r_sun / r_sun.norm();
151
152        // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor.
153        let occult = self
154            .shadow_model
155            .compute(osc, almanac)
156            .context(DynamicsAlmanacSnafu {
157                action: "solar radiation pressure computation",
158            })?
159            .factor();
160
161        // Compute the illumination factor.
162        let k: f64 = (occult - 1.0).abs();
163
164        let r_sun_au = r_sun.norm() / AU;
165        // in N/(m^2)
166        let flux_pressure = (k * self.phi / SPEED_OF_LIGHT_M_S) * (1.0 / r_sun_au).powi(2);
167
168        // Note the 1e-3 is to convert the SRP from m/s^2 to km/s^2
169        Ok(1e-3 * ctx.srp.coeff_reflectivity * ctx.srp.area_m2 * flux_pressure * r_sun_unit)
170    }
171
172    fn gradient(
173        &self,
174        ctx: &Spacecraft,
175        almanac: &Almanac,
176    ) -> Result<(Vector3<f64>, Matrix4x3<f64>), DynamicsError> {
177        let osc = ctx.orbit;
178
179        // Compute the position of the Sun as seen from the spacecraft
180        let r_sun = almanac
181            .transform_to(ctx.orbit, self.shadow_model.light_source, None)
182            .context(DynamicsAlmanacSnafu {
183                action: "transforming state to vector seen from Sun",
184            })?
185            .radius_km;
186
187        let r_sun_d: Vector3<OHyperdual<f64, Const<9>>> = hyperspace_from_vector(&r_sun);
188        let r_sun_unit = r_sun_d / norm(&r_sun_d);
189
190        // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor.
191        let occult = self
192            .shadow_model
193            .compute(osc, almanac)
194            .context(DynamicsAlmanacSnafu {
195                action: "solar radiation pressure computation",
196            })?
197            .factor();
198
199        // Compute the illumination factor.
200        let k: f64 = (occult - 1.0).abs();
201
202        let r_sun_au = norm(&r_sun_d) / AU;
203        let inv_r_sun_au = OHyperdual::<f64, Const<9>>::from_real(1.0) / (r_sun_au);
204        let inv_r_sun_au_p2 = inv_r_sun_au.powi(2);
205        // in N/(m^2)
206        let flux_pressure =
207            OHyperdual::<f64, Const<9>>::from_real(k * self.phi / SPEED_OF_LIGHT_M_S)
208                * inv_r_sun_au_p2;
209
210        // Note the 1e-3 is to convert the SRP from m/s^2 to km/s^2
211        let dual_force_scalar = OHyperdual::<f64, Const<9>>::from_real(
212            1e-3 * ctx.srp.coeff_reflectivity * ctx.srp.area_m2,
213        );
214        let mut dual_force: Vector3<OHyperdual<f64, Const<9>>> = Vector3::zeros();
215        dual_force[0] = dual_force_scalar * flux_pressure * r_sun_unit[0];
216        dual_force[1] = dual_force_scalar * flux_pressure * r_sun_unit[1];
217        dual_force[2] = dual_force_scalar * flux_pressure * r_sun_unit[2];
218
219        // Extract result into Vector6 and Matrix6
220        let mut dx = Vector3::zeros();
221        let mut grad = Matrix4x3::zeros();
222        for i in 0..3 {
223            dx[i] += dual_force[i].real();
224            // NOTE: Although the hyperdual state is of size 7, we're only setting the values up to 3 (Matrix3)
225            for j in 0..3 {
226                grad[(i, j)] += dual_force[i][j + 1];
227            }
228        }
229
230        // Compute the partial wrt to Cr.
231        let wrt_cr = self.eom(ctx, almanac)? / ctx.srp.coeff_reflectivity;
232        for j in 0..3 {
233            grad[(3, j)] = wrt_cr[j];
234        }
235
236        Ok((dx, grad))
237    }
238}
239
240impl fmt::Display for SolarPressure {
241    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
242        write!(
243            f,
244            "SRP with φ = {} W/m^2 and eclipse {}",
245            self.phi, self.shadow_model
246        )
247    }
248}
249
250#[cfg(feature = "python")]
251#[cfg_attr(feature = "python", pymethods)]
252impl SolarPressure {
253    #[pyo3(signature = (shadow_bodies, almanac, flux_w_m2=SOLAR_FLUX_W_m2, estimate=true))]
254    #[new]
255    fn py_new(
256        shadow_bodies: Vec<Frame>,
257        almanac: &Almanac,
258        flux_w_m2: f64,
259        estimate: bool,
260    ) -> Result<Self, DynamicsError> {
261        let mut me = Self::default_flux_raw(shadow_bodies, almanac)?;
262        me.phi = flux_w_m2;
263        me.estimate = estimate;
264
265        Ok(me)
266    }
267
268    fn __str__(&self) -> String {
269        format!("{self}")
270    }
271
272    fn __repr__(&self) -> String {
273        format!("{self} @ {self:p}")
274    }
275}
nyx_space/dynamics/solarpressure.rs

nyx_space/dynamics/
solarpressure.rs
