nyx_space/dynamics/
drag.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 anise::almanac::Almanac;
20use anise::constants::frames::IAU_EARTH_FRAME;
21use snafu::ResultExt;
22
23use super::{
24    DynamicsAlmanacSnafu, DynamicsAstroSnafu, DynamicsError, DynamicsPlanetarySnafu, ForceModel,
25};
26use crate::cosmic::{AstroError, AstroPhysicsSnafu, Frame, Spacecraft};
27use crate::linalg::{Matrix4x3, Vector3};
28use serde::{Deserialize, Serialize};
29use serde_dhall::StaticType;
30use std::fmt;
31use std::sync::Arc;
32
33/// Density in kg/m^3 and altitudes in meters, not kilometers!
34#[derive(Clone, Copy, Debug, Serialize, Deserialize, StaticType)]
35pub enum AtmDensity {
36    Constant(f64),
37    Exponential { rho0: f64, r0: f64, ref_alt_m: f64 },
38    StdAtm { max_alt_m: f64 },
39}
40
41/// `ConstantDrag` implements a constant drag model as defined in Vallado, 4th ed., page 551, with an important caveat.
42///
43/// **WARNING:** This basic model assumes that the velocity of the spacecraft is identical to the velocity of the upper atmosphere.
44/// This is a **bad** assumption and **should not** be used for high fidelity simulations.
45/// This will be resolved after https://github.com/nyx-space/nyx/issues/317 is implemented.
46#[derive(Clone)]
47pub struct ConstantDrag {
48    /// atmospheric density in kg/m^3
49    pub rho: f64,
50    /// Frame causing the drag
51    pub drag_frame: Frame,
52    /// Set to true to estimate the coefficient of drag
53    pub estimate: bool,
54}
55
56impl fmt::Display for ConstantDrag {
57    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
58        write!(
59            f,
60            "\tConstant Drag rho = {} kg/m^3 in frame {}",
61            self.rho, self.drag_frame
62        )
63    }
64}
65
66impl ForceModel for ConstantDrag {
67    fn estimation_index(&self) -> Option<usize> {
68        if self.estimate { Some(7) } else { None }
69    }
70
71    fn eom(&self, ctx: &Spacecraft, almanac: Arc<Almanac>) -> Result<Vector3<f64>, DynamicsError> {
72        let osc = almanac
73            .transform_to(ctx.orbit, self.drag_frame, None)
74            .context(DynamicsAlmanacSnafu {
75                action: "transforming into drag frame",
76            })?;
77
78        let velocity = osc.velocity_km_s;
79        // Note the 1e3 factor to convert drag units from ((kg * km^2 * s^-2) / m^1) to (kg * km * s^-2)
80        Ok(-0.5
81            * 1e3
82            * self.rho
83            * ctx.drag.coeff_drag
84            * ctx.drag.area_m2
85            * velocity.norm()
86            * velocity)
87    }
88
89    fn gradient(
90        &self,
91        _osc_ctx: &Spacecraft,
92        _almanac: Arc<Almanac>,
93    ) -> Result<(Vector3<f64>, Matrix4x3<f64>), DynamicsError> {
94        Err(DynamicsError::DynamicsAstro {
95            source: AstroError::PartialsUndefined,
96        })
97    }
98}
99
100/// `Drag` implements all three drag models.
101#[derive(Copy, Clone, Debug, Serialize, Deserialize, StaticType)]
102pub struct Drag {
103    /// Density computation method
104    pub density: AtmDensity,
105    /// Frame to compute the drag in
106    pub drag_frame: Frame,
107    /// Set to true to estimate the coefficient of drag
108    pub estimate: bool,
109}
110
111impl Drag {
112    /// Common exponential drag model for the Earth
113    pub fn earth_exp(almanac: Arc<Almanac>) -> Result<Arc<Self>, DynamicsError> {
114        Ok(Arc::new(Self {
115            density: AtmDensity::Exponential {
116                rho0: 3.614e-13,
117                r0: 700_000.0,
118                ref_alt_m: 88_667.0,
119            },
120            drag_frame: almanac.frame_info(IAU_EARTH_FRAME).context({
121                DynamicsPlanetarySnafu {
122                    action: "planetary data from third body not loaded",
123                }
124            })?,
125            estimate: false,
126        }))
127    }
128
129    /// Drag model which uses the standard atmosphere 1976 model for atmospheric density
130    pub fn std_atm1976(almanac: Arc<Almanac>) -> Result<Arc<Self>, DynamicsError> {
131        Ok(Arc::new(Self {
132            density: AtmDensity::StdAtm {
133                max_alt_m: 1_000_000.0,
134            },
135            drag_frame: almanac.frame_info(IAU_EARTH_FRAME).context({
136                DynamicsPlanetarySnafu {
137                    action: "planetary data from third body not loaded",
138                }
139            })?,
140            estimate: false,
141        }))
142    }
143}
144
145impl fmt::Display for Drag {
146    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
147        write!(
148            f,
149            "\tDrag density {:?} in frame {}",
150            self.density, self.drag_frame
151        )
152    }
153}
154
155impl ForceModel for Drag {
156    fn estimation_index(&self) -> Option<usize> {
157        if self.estimate { Some(7) } else { None }
158    }
159
160    fn eom(&self, ctx: &Spacecraft, almanac: Arc<Almanac>) -> Result<Vector3<f64>, DynamicsError> {
161        let integration_frame = ctx.orbit.frame;
162
163        let osc_drag_frame = almanac
164            .transform_to(ctx.orbit, self.drag_frame, None)
165            .context(DynamicsAlmanacSnafu {
166                action: "transforming into drag frame",
167            })?;
168
169        match self.density {
170            AtmDensity::Constant(rho) => {
171                let velocity = osc_drag_frame.velocity_km_s;
172                // Note the 1e3 factor to convert drag units from ((kg * km^2 * s^-2) / m^1) to (kg * km * s^-2)
173                Ok(-0.5
174                    * 1e3
175                    * rho
176                    * ctx.drag.coeff_drag
177                    * ctx.drag.area_m2
178                    * velocity.norm()
179                    * velocity)
180            }
181
182            AtmDensity::Exponential {
183                rho0,
184                r0,
185                ref_alt_m,
186            } => {
187                // Compute rho in the drag frame.
188                let rho = rho0
189                    * (-(osc_drag_frame.rmag_km()
190                        - (r0
191                            + self
192                                .drag_frame
193                                .mean_equatorial_radius_km()
194                                .context(AstroPhysicsSnafu)
195                                .context(DynamicsAstroSnafu)?))
196                        / ref_alt_m)
197                        .exp();
198
199                // TODO: Drag modeling will be improved in https://github.com/nyx-space/nyx/issues/317
200                // The frame will be double checked in this PR as well.
201                let velocity_integr_frame = almanac
202                    .transform_to(osc_drag_frame, integration_frame, None)
203                    .context(DynamicsAlmanacSnafu {
204                        action: "rotating into the integration frame",
205                    })?
206                    .velocity_km_s;
207
208                let velocity = velocity_integr_frame - osc_drag_frame.velocity_km_s;
209                // Note the 1e3 factor to convert drag units from ((kg * km^2 * s^-2) / m^1) to (kg * km * s^-2)
210                Ok(-0.5
211                    * 1e3
212                    * rho
213                    * ctx.drag.coeff_drag
214                    * ctx.drag.area_m2
215                    * velocity.norm()
216                    * velocity)
217            }
218
219            AtmDensity::StdAtm { max_alt_m } => {
220                let altitude_km = osc_drag_frame.rmag_km()
221                    - self
222                        .drag_frame
223                        .mean_equatorial_radius_km()
224                        .context(AstroPhysicsSnafu)
225                        .context(DynamicsAstroSnafu)?;
226                let rho = if altitude_km > max_alt_m / 1_000.0 {
227                    // Use a constant density
228                    10.0_f64.powf((-7e-5) * altitude_km - 14.464)
229                } else {
230                    // Code from AVS/Schaub's Basilisk
231                    // Calculating the density based on a scaled 6th order polynomial fit to the log of density
232                    let scale = (altitude_km - 526.8000) / 292.8563;
233                    let logdensity =
234                        0.34047 * scale.powi(6) - 0.5889 * scale.powi(5) - 0.5269 * scale.powi(4)
235                            + 1.0036 * scale.powi(3)
236                            + 0.60713 * scale.powi(2)
237                            - 2.3024 * scale
238                            - 12.575;
239
240                    /* Calculating density by raising 10 to the log of density */
241                    10.0_f64.powf(logdensity)
242                };
243
244                let velocity_integr_frame = almanac
245                    .transform_to(osc_drag_frame, integration_frame, None)
246                    .context(DynamicsAlmanacSnafu {
247                        action: "rotating into the integration frame",
248                    })?
249                    .velocity_km_s;
250
251                let velocity = velocity_integr_frame - osc_drag_frame.velocity_km_s;
252                // Note the 1e3 factor to convert drag units from ((kg * km^2 * s^-2) / m^1) to (kg * km * s^-2)
253                Ok(-0.5
254                    * 1e3
255                    * rho
256                    * ctx.drag.coeff_drag
257                    * ctx.drag.area_m2
258                    * velocity.norm()
259                    * velocity)
260            }
261        }
262    }
263
264    fn gradient(
265        &self,
266        _osc_ctx: &Spacecraft,
267        _almanac: Arc<Almanac>,
268    ) -> Result<(Vector3<f64>, Matrix4x3<f64>), DynamicsError> {
269        Err(DynamicsError::DynamicsAstro {
270            source: AstroError::PartialsUndefined,
271        })
272    }
273}
nyx_space/dynamics/drag.rs

nyx_space/dynamics/
drag.rs
