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/// 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 { Some(6) } else { None }
135    }
136
137    fn eom(&self, ctx: &Spacecraft, almanac: Arc<Almanac>) -> Result<Vector3<f64>, DynamicsError> {
138        let osc = ctx.orbit;
139        // Compute the position of the Sun as seen from the spacecraft
140        let r_sun = almanac
141            .transform_to(ctx.orbit, self.shadow_model.light_source, None)
142            .context(DynamicsAlmanacSnafu {
143                action: "transforming state to vector seen from Sun",
144            })?
145            .radius_km;
146
147        let r_sun_unit = r_sun / r_sun.norm();
148
149        // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor.
150        let occult = self
151            .shadow_model
152            .compute(osc, almanac)
153            .context(DynamicsAlmanacSnafu {
154                action: "solar radiation pressure computation",
155            })?
156            .factor();
157
158        // Compute the illumination factor.
159        let k: f64 = (occult - 1.0).abs();
160
161        let r_sun_au = r_sun.norm() / AU;
162        // in N/(m^2)
163        let flux_pressure = (k * self.phi / SPEED_OF_LIGHT_M_S) * (1.0 / r_sun_au).powi(2);
164
165        // Note the 1e-3 is to convert the SRP from m/s^2 to km/s^2
166        Ok(1e-3 * ctx.srp.coeff_reflectivity * ctx.srp.area_m2 * flux_pressure * r_sun_unit)
167    }
168
169    fn gradient(
170        &self,
171        ctx: &Spacecraft,
172        almanac: Arc<Almanac>,
173    ) -> Result<(Vector3<f64>, Matrix4x3<f64>), DynamicsError> {
174        let osc = ctx.orbit;
175
176        // Compute the position of the Sun as seen from the spacecraft
177        let r_sun = almanac
178            .transform_to(ctx.orbit, self.shadow_model.light_source, None)
179            .context(DynamicsAlmanacSnafu {
180                action: "transforming state to vector seen from Sun",
181            })?
182            .radius_km;
183
184        let r_sun_d: Vector3<OHyperdual<f64, Const<9>>> = hyperspace_from_vector(&r_sun);
185        let r_sun_unit = r_sun_d / norm(&r_sun_d);
186
187        // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor.
188        let occult = self
189            .shadow_model
190            .compute(osc, almanac.clone())
191            .context(DynamicsAlmanacSnafu {
192                action: "solar radiation pressure computation",
193            })?
194            .factor();
195
196        // Compute the illumination factor.
197        let k: f64 = (occult - 1.0).abs();
198
199        let r_sun_au = norm(&r_sun_d) / AU;
200        let inv_r_sun_au = OHyperdual::<f64, Const<9>>::from_real(1.0) / (r_sun_au);
201        let inv_r_sun_au_p2 = inv_r_sun_au.powi(2);
202        // in N/(m^2)
203        let flux_pressure =
204            OHyperdual::<f64, Const<9>>::from_real(k * self.phi / SPEED_OF_LIGHT_M_S)
205                * inv_r_sun_au_p2;
206
207        // Note the 1e-3 is to convert the SRP from m/s^2 to km/s^2
208        let dual_force_scalar = OHyperdual::<f64, Const<9>>::from_real(
209            1e-3 * ctx.srp.coeff_reflectivity * ctx.srp.area_m2,
210        );
211        let mut dual_force: Vector3<OHyperdual<f64, Const<9>>> = Vector3::zeros();
212        dual_force[0] = dual_force_scalar * flux_pressure * r_sun_unit[0];
213        dual_force[1] = dual_force_scalar * flux_pressure * r_sun_unit[1];
214        dual_force[2] = dual_force_scalar * flux_pressure * r_sun_unit[2];
215
216        // Extract result into Vector6 and Matrix6
217        let mut dx = Vector3::zeros();
218        let mut grad = Matrix4x3::zeros();
219        for i in 0..3 {
220            dx[i] += dual_force[i].real();
221            // NOTE: Although the hyperdual state is of size 7, we're only setting the values up to 3 (Matrix3)
222            for j in 0..3 {
223                grad[(i, j)] += dual_force[i][j + 1];
224            }
225        }
226
227        // Compute the partial wrt to Cr.
228        let wrt_cr = self.eom(ctx, almanac)? / ctx.srp.coeff_reflectivity;
229        for j in 0..3 {
230            grad[(3, j)] = wrt_cr[j];
231        }
232
233        Ok((dx, grad))
234    }
235}
236
237impl fmt::Display for SolarPressure {
238    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
239        write!(
240            f,
241            "SRP with φ = {} W/m^2 and eclipse {}",
242            self.phi, self.shadow_model
243        )
244    }
245}
nyx_space/dynamics/solarpressure.rs

nyx_space/dynamics/
solarpressure.rs
