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

nyx_space/dynamics/
solarpressure.rs
