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::{Frame, Spacecraft, AU, SPEED_OF_LIGHT_M_S};
22use crate::linalg::{Const, Matrix4x3, Vector3};
23use anise::almanac::Almanac;
24use anise::constants::frames::{EARTH_J2000, SUN_J2000};
25use hyperdual::{hyperspace_from_vector, linalg::norm, Float, OHyperdual};
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/// Computation of solar radiation pressure is based on STK: <http://help.agi.com/stk/index.htm#gator/eq-solar.htm> .
38#[derive(Clone, Debug, Serialize, Deserialize, StaticType)]
39pub struct SolarPressure {
40    /// solar flux at 1 AU, in W/m^2
41    pub phi: f64,
42    pub shadow_model: ShadowModel,
43    /// Set to true to estimate the coefficient of reflectivity
44    pub estimate: bool,
45}
46
47impl Default for SolarPressure {
48    /// Important: the default will FAIL at runtime if the shadow model is not manually defined with loaded frames.
49    fn default() -> Self {
50        Self {
51            phi: SOLAR_FLUX_W_m2,
52            estimate: false,
53            shadow_model: ShadowModel {
54                light_source: SUN_J2000,
55                shadow_bodies: vec![EARTH_J2000],
56            },
57        }
58    }
59}
60
61impl SolarPressure {
62    /// Will set the solar flux at 1 AU to: Phi = 1367.0
63    fn default_flux_raw(
64        shadow_bodies: Vec<Frame>,
65        almanac: Arc<Almanac>,
66    ) -> Result<Self, DynamicsError> {
67        let e_loc = ShadowModel {
68            light_source: almanac.frame_info(SUN_J2000).context({
69                DynamicsPlanetarySnafu {
70                    action: "planetary data from third body not loaded",
71                }
72            })?,
73            shadow_bodies: shadow_bodies
74                .iter()
75                .filter_map(|object| match almanac.frame_info(object) {
76                    Ok(loaded_obj) => Some(loaded_obj),
77                    Err(e) => {
78                        warn!("when initializing SRP model for {object}, {e}");
79                        None
80                    }
81                })
82                .collect(),
83        };
84        Ok(Self {
85            phi: SOLAR_FLUX_W_m2,
86            shadow_model: e_loc,
87            estimate: true,
88        })
89    }
90
91    /// Accounts for the shadowing of only one body and will set the solar flux at 1 AU to: Phi = 1367.0
92    pub fn default_flux(
93        shadow_body: Frame,
94        almanac: Arc<Almanac>,
95    ) -> Result<Arc<Self>, DynamicsError> {
96        Ok(Arc::new(Self::default_flux_raw(
97            vec![shadow_body],
98            almanac,
99        )?))
100    }
101
102    /// Accounts for the shadowing of only one body and will set the solar flux at 1 AU to: Phi = 1367.0
103    pub fn default_no_estimation(
104        shadow_bodies: Vec<Frame>,
105        almanac: Arc<Almanac>,
106    ) -> Result<Arc<Self>, DynamicsError> {
107        let mut srp = Self::default_flux_raw(shadow_bodies, almanac)?;
108        srp.estimate = false;
109        Ok(Arc::new(srp))
110    }
111
112    /// Must provide the flux in W/m^2
113    pub fn with_flux(
114        flux_w_m2: f64,
115        shadow_bodies: Vec<Frame>,
116        almanac: Arc<Almanac>,
117    ) -> Result<Arc<Self>, DynamicsError> {
118        let mut me = Self::default_flux_raw(shadow_bodies, almanac)?;
119        me.phi = flux_w_m2;
120        Ok(Arc::new(me))
121    }
122
123    /// Solar radiation pressure force model accounting for the provided shadow bodies.
124    pub fn new(
125        shadow_bodies: Vec<Frame>,
126        almanac: Arc<Almanac>,
127    ) -> Result<Arc<Self>, DynamicsError> {
128        Ok(Arc::new(Self::default_flux_raw(shadow_bodies, almanac)?))
129    }
130}
131
132impl ForceModel for SolarPressure {
133    fn estimation_index(&self) -> Option<usize> {
134        if self.estimate {
135            Some(6)
136        } else {
137            None
138        }
139    }
140
141    fn eom(&self, ctx: &Spacecraft, almanac: Arc<Almanac>) -> Result<Vector3<f64>, DynamicsError> {
142        let osc = ctx.orbit;
143        // Compute the position of the Sun as seen from the spacecraft
144        let r_sun = almanac
145            .transform_to(ctx.orbit, self.shadow_model.light_source, None)
146            .context(DynamicsAlmanacSnafu {
147                action: "transforming state to vector seen from Sun",
148            })?
149            .radius_km;
150
151        let r_sun_unit = r_sun / r_sun.norm();
152
153        // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor.
154        let occult = self
155            .shadow_model
156            .compute(osc, almanac)
157            .context(DynamicsAlmanacSnafu {
158                action: "solar radiation pressure computation",
159            })?
160            .factor();
161
162        // Compute the illumination factor.
163        let k: f64 = (occult - 1.0).abs();
164
165        let r_sun_au = r_sun.norm() / AU;
166        // in N/(m^2)
167        let flux_pressure = (k * self.phi / SPEED_OF_LIGHT_M_S) * (1.0 / r_sun_au).powi(2);
168
169        // Note the 1e-3 is to convert the SRP from m/s^2 to km/s^2
170        Ok(1e-3 * ctx.srp.coeff_reflectivity * ctx.srp.area_m2 * flux_pressure * r_sun_unit)
171    }
172
173    fn dual_eom(
174        &self,
175        ctx: &Spacecraft,
176        almanac: Arc<Almanac>,
177    ) -> Result<(Vector3<f64>, Matrix4x3<f64>), DynamicsError> {
178        let osc = ctx.orbit;
179
180        // Compute the position of the Sun as seen from the spacecraft
181        let r_sun = almanac
182            .transform_to(ctx.orbit, self.shadow_model.light_source, None)
183            .context(DynamicsAlmanacSnafu {
184                action: "transforming state to vector seen from Sun",
185            })?
186            .radius_km;
187
188        let r_sun_d: Vector3<OHyperdual<f64, Const<9>>> = hyperspace_from_vector(&r_sun);
189        let r_sun_unit = r_sun_d / norm(&r_sun_d);
190
191        // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor.
192        let occult = self
193            .shadow_model
194            .compute(osc, almanac.clone())
195            .context(DynamicsAlmanacSnafu {
196                action: "solar radiation pressure computation",
197            })?
198            .factor();
199
200        // Compute the illumination factor.
201        let k: f64 = (occult - 1.0).abs();
202
203        let r_sun_au = norm(&r_sun_d) / AU;
204        let inv_r_sun_au = OHyperdual::<f64, Const<9>>::from_real(1.0) / (r_sun_au);
205        let inv_r_sun_au_p2 = inv_r_sun_au.powi(2);
206        // in N/(m^2)
207        let flux_pressure =
208            OHyperdual::<f64, Const<9>>::from_real(k * self.phi / SPEED_OF_LIGHT_M_S)
209                * inv_r_sun_au_p2;
210
211        // Note the 1e-3 is to convert the SRP from m/s^2 to km/s^2
212        let dual_force_scalar = OHyperdual::<f64, Const<9>>::from_real(
213            1e-3 * ctx.srp.coeff_reflectivity * ctx.srp.area_m2,
214        );
215        let mut dual_force: Vector3<OHyperdual<f64, Const<9>>> = Vector3::zeros();
216        dual_force[0] = dual_force_scalar * flux_pressure * r_sun_unit[0];
217        dual_force[1] = dual_force_scalar * flux_pressure * r_sun_unit[1];
218        dual_force[2] = dual_force_scalar * flux_pressure * r_sun_unit[2];
219
220        // Extract result into Vector6 and Matrix6
221        let mut dx = Vector3::zeros();
222        let mut grad = Matrix4x3::zeros();
223        for i in 0..3 {
224            dx[i] += dual_force[i].real();
225            // NOTE: Although the hyperdual state is of size 7, we're only setting the values up to 3 (Matrix3)
226            for j in 0..3 {
227                grad[(i, j)] += dual_force[i][j + 1];
228            }
229        }
230
231        // Compute the partial wrt to Cr.
232        let wrt_cr = self.eom(ctx, almanac)? / ctx.srp.coeff_reflectivity;
233        for j in 0..3 {
234            grad[(3, j)] = wrt_cr[j];
235        }
236
237        Ok((dx, grad))
238    }
239}
240
241impl fmt::Display for SolarPressure {
242    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
243        write!(
244            f,
245            "SRP with φ = {} W/m^2 and eclipse {}",
246            self.phi, self.shadow_model
247        )
248    }
249}
nyx_space/dynamics/solarpressure.rs

nyx_space/dynamics/
solarpressure.rs
