nyx_space/dynamics/drag/nrlmsise00/
mod.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
19//! NRLMSISE-00 empirical atmosphere model.
20//!
21//! Clean-room implementation based on the following references:
22//! - Picone, J.M. et al. (2002), "NRLMSISE-00 empirical model of the atmosphere:
23//!   Statistical comparisons and scientific issues", J. Geophys. Res., 107(A12), 1468,
24//!   doi:10.1029/2002JA009430
25//! - Hedin, A.E. (1991), "Extension of the MSIS thermosphere model into the middle
26//!   and lower atmosphere", J. Geophys. Res., 96(A2), 1159-1172.
27//! - Hedin, A.E. (1987), "MSIS-86 thermospheric model",
28//!   J. Geophys. Res., 92(A5), 4649-4662.
29//!
30//! Coefficient values are from the official NRL distribution, treated as published data.
31//! NRLMSISE-00 is believed to be in the public domain as a U.S. Government work
32//! (17 U.S.C. § 105), though no explicit license was provided by NRL.
33//!
34//! Note: MSIS is a registered trademark. This module uses the name "NRLMSISE-00"
35//! for nominative fair use (identifying compatibility with the NRL model).
36//!
37//! Validated against `pymsis` (official NRL Fortran wrapper, `version=0`).
38
39use crate::dynamics::DynamicsError;
40pub use crate::io::space_weather::Msise00DailyWeather;
41use hifitime::Epoch;
42use serde::{Deserialize, Serialize};
43use serde_dhall::StaticType;
44
45#[cfg(feature = "python")]
46use pyo3::prelude::*;
47
48pub mod coefficients;
49mod model;
50
51/// Full output of the NRLMSISE-00 model.
52///
53/// Includes temperatures and all species number densities.
54#[derive(Debug, Clone)]
55pub struct Nrlmsise00Output {
56    /// Exospheric temperature [K].
57    pub temp_exo_k: f64,
58    /// Temperature at altitude [K].
59    pub temp_alt_k: f64,
60    /// He number density [cm⁻³].
61    pub density_he_per_cm3: f64,
62    /// O number density [cm⁻³].
63    pub density_o_per_cm3: f64,
64    /// N₂ number density [cm⁻³].
65    pub density_n2_per_cm3: f64,
66    /// O₂ number density [cm⁻³].
67    pub density_o2_per_cm3: f64,
68    /// Ar number density [cm⁻³].
69    pub density_ar_per_cm3: f64,
70    /// H number density [cm⁻³].
71    pub density_h_per_cm3: f64,
72    /// N number density [cm⁻³].
73    pub density_n_per_cm3: f64,
74    /// Anomalous oxygen number density [cm⁻³].
75    pub density_anomalous_o_per_cm3: f64,
76    /// Total mass density [kg/m³].
77    pub total_mass_density_kg_m3: f64,
78}
79
80/// Input parameters for a single NRLMSISE-00 evaluation.
81#[derive(Debug, Clone)]
82pub struct Nrlmsise00Input {
83    /// Day of year [1-366].
84    pub day_of_year: u32,
85    /// Universal time [seconds since midnight].
86    pub ut_seconds: f64,
87    /// Geodetic altitude [km].
88    pub altitude_km: f64,
89    /// Geodetic latitude [degrees, -90 to 90].
90    pub latitude_deg: f64,
91    /// Geodetic longitude [degrees, 0 to 360 or -180 to 180].
92    pub longitude_deg: f64,
93    /// Local apparent solar time [hours, 0-24].
94    pub local_solar_time_hours: f64,
95    /// Previous day's F10.7 [SFU].
96    pub f107_daily: f64,
97    /// 81-day centered average F10.7 [SFU].
98    pub f107_avg: f64,
99    /// Daily Ap index.
100    pub ap_daily: f64,
101    /// 7-element Ap array for magnetic activity variations.
102    pub ap_array: [f64; 7],
103}
104
105#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, StaticType)]
106#[cfg_attr(feature = "python", pyclass(from_py_object, eq, eq_int))]
107pub enum GeomagneticMode {
108    Off,
109    StandardDailyAp,
110    ExtendedHistory57h, // Sets sw[9] = -1.0
111}
112
113#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, StaticType)]
114#[cfg_attr(feature = "python", pyclass(from_py_object, get_all, set_all))]
115pub struct Nrlmsise00Flags {
116    pub geomagnetic: GeomagneticMode,
117    pub f107_solar_flux: bool,
118    pub time_independent: bool,
119    pub annual_harmonics: bool,
120    pub semiannual_harmonics: bool,
121    pub diurnal_tides: bool,
122    pub semidiurnal_tides: bool,
123    pub terdiurnal_tides: bool,
124    pub ut_and_longitude: bool,
125    pub exospheric_temp_variations: bool,
126    pub lower_boundary_temp_variations: bool,
127    pub gradient_variations: bool,
128    pub departures_from_diffusive_equilibrium: bool,
129    pub lower_thermosphere_temp_variations: bool,
130    pub upper_stratosphere_temp_variations: bool,
131    pub boundary_density_variations: bool,
132    pub lower_mesosphere_temp_variations: bool,
133    pub turbopause_scale_height_variations: bool,
134}
135
136impl Default for Nrlmsise00Flags {
137    fn default() -> Self {
138        Self {
139            geomagnetic: GeomagneticMode::StandardDailyAp,
140            f107_solar_flux: true,
141            time_independent: true,
142            annual_harmonics: true,
143            semiannual_harmonics: true,
144            diurnal_tides: true,
145            semidiurnal_tides: true,
146            terdiurnal_tides: true,
147            ut_and_longitude: true,
148            exospheric_temp_variations: true,
149            lower_boundary_temp_variations: true,
150            gradient_variations: true,
151            departures_from_diffusive_equilibrium: true,
152            lower_thermosphere_temp_variations: true,
153            upper_stratosphere_temp_variations: true,
154            boundary_density_variations: true,
155            lower_mesosphere_temp_variations: true,
156            turbopause_scale_height_variations: true,
157        }
158    }
159}
160
161impl Nrlmsise00Flags {
162    /// Compiles the high-level flags into the raw 24-element float array consumed by the kernel.
163    pub(crate) fn to_switches(self) -> [f64; 24] {
164        let mut sw = [1.0f64; 24];
165
166        // NOTE Unit selection is ALWAYS set to 1.0 because the calculation code does not even check it.
167
168        // Geomagnetic storm mode
169        sw[9] = match self.geomagnetic {
170            GeomagneticMode::Off => 0.0,
171            GeomagneticMode::StandardDailyAp => 1.0,
172            GeomagneticMode::ExtendedHistory57h => -1.0,
173        };
174
175        // Specific feature toggles
176        if !self.f107_solar_flux {
177            sw[1] = 0.0;
178        }
179        if !self.time_independent {
180            sw[2] = 0.0;
181        }
182        if !self.annual_harmonics {
183            sw[3] = 0.0;
184            sw[5] = 0.0;
185        }
186        if !self.semiannual_harmonics {
187            sw[4] = 0.0;
188            sw[6] = 0.0;
189        }
190        if !self.diurnal_tides {
191            sw[7] = 0.0;
192        }
193        if !self.semidiurnal_tides {
194            sw[8] = 0.0;
195        }
196        if !self.terdiurnal_tides {
197            sw[14] = 0.0;
198        }
199        if !self.ut_and_longitude {
200            sw[10] = 0.0;
201            sw[11] = 0.0;
202            sw[12] = 0.0;
203            sw[13] = 0.0;
204        }
205        if !self.exospheric_temp_variations {
206            sw[16] = 0.0;
207        }
208        if !self.lower_boundary_temp_variations {
209            sw[17] = 0.0;
210        }
211        if !self.gradient_variations {
212            sw[19] = 0.0;
213        }
214        if !self.departures_from_diffusive_equilibrium {
215            sw[15] = 0.0;
216        }
217        if !self.lower_thermosphere_temp_variations {
218            sw[18] = 0.0;
219        }
220        if !self.upper_stratosphere_temp_variations {
221            sw[20] = 0.0;
222        }
223        if !self.boundary_density_variations {
224            sw[21] = 0.0;
225        }
226        if !self.lower_mesosphere_temp_variations {
227            sw[22] = 0.0;
228        }
229        if !self.turbopause_scale_height_variations {
230            sw[23] = 0.0;
231        }
232
233        sw
234    }
235}
236
237#[cfg(feature = "python")]
238#[pymethods]
239impl Nrlmsise00Flags {
240    #[new]
241    #[pyo3(signature = (
242        geomagnetic = None,
243        f107_solar_flux = true,
244        time_independent = true,
245        annual_harmonics = true,
246        semiannual_harmonics = true,
247        diurnal_tides = true,
248        semidiurnal_tides = true,
249        terdiurnal_tides = true,
250        ut_and_longitude = true,
251        exospheric_temp_variations = true,
252        lower_boundary_temp_variations = true,
253        gradient_variations = true,
254        departures_from_diffusive_equilibrium = true,
255        lower_thermosphere_temp_variations = true,
256        upper_stratosphere_temp_variations = true,
257        boundary_density_variations = true,
258        lower_mesosphere_temp_variations = true,
259        turbopause_scale_height_variations = true,
260    ))]
261    #[allow(clippy::too_many_arguments)]
262    fn py_new(
263        geomagnetic: Option<GeomagneticMode>,
264        f107_solar_flux: bool,
265        time_independent: bool,
266        annual_harmonics: bool,
267        semiannual_harmonics: bool,
268        diurnal_tides: bool,
269        semidiurnal_tides: bool,
270        terdiurnal_tides: bool,
271        ut_and_longitude: bool,
272        exospheric_temp_variations: bool,
273        lower_boundary_temp_variations: bool,
274        gradient_variations: bool,
275        departures_from_diffusive_equilibrium: bool,
276        lower_thermosphere_temp_variations: bool,
277        upper_stratosphere_temp_variations: bool,
278        boundary_density_variations: bool,
279        lower_mesosphere_temp_variations: bool,
280        turbopause_scale_height_variations: bool,
281    ) -> Self {
282        Self {
283            geomagnetic: geomagnetic.unwrap_or(GeomagneticMode::StandardDailyAp),
284            f107_solar_flux,
285            time_independent,
286            annual_harmonics,
287            semiannual_harmonics,
288            diurnal_tides,
289            semidiurnal_tides,
290            terdiurnal_tides,
291            ut_and_longitude,
292            exospheric_temp_variations,
293            lower_boundary_temp_variations,
294            gradient_variations,
295            departures_from_diffusive_equilibrium,
296            lower_thermosphere_temp_variations,
297            upper_stratosphere_temp_variations,
298            boundary_density_variations,
299            lower_mesosphere_temp_variations,
300            turbopause_scale_height_variations,
301        }
302    }
303
304    fn __repr__(&self) -> String {
305        format!("{:?}", self)
306    }
307
308    fn __str__(&self) -> String {
309        format!("{:?} @ {self:p}", self)
310    }
311}
312
313/// Compute full NRLMSISE-00 output for the given input parameters.
314///
315/// Returns temperatures and all species number densities.
316fn calculate(input: &Nrlmsise00Input, flags: Nrlmsise00Flags) -> Nrlmsise00Output {
317    let sw = flags.to_switches();
318    let (d, temp_exo, temp_alt) = model::compute(input, &sw);
319    // d[0..8]: He, O, N2, O2, Ar, total_mass(g/cm³), H, N, anomO
320    // Total mass density: d[5] is in g/cm³, convert to kg/m³ (* 1000)
321    Nrlmsise00Output {
322        temp_exo_k: temp_exo,
323        temp_alt_k: temp_alt,
324        density_he_per_cm3: d[0],
325        density_o_per_cm3: d[1],
326        density_n2_per_cm3: d[2],
327        density_o2_per_cm3: d[3],
328        density_ar_per_cm3: d[4],
329        density_h_per_cm3: d[6],
330        density_n_per_cm3: d[7],
331        density_anomalous_o_per_cm3: d[8],
332        // Convert g/cm^3 to kg/m^3: g->kg <=> 1e-3; cm^3 -> m^3 <-> 1e6 => 1e3
333        total_mass_density_kg_m3: d[5] * 1e3,
334    }
335}
336
337/// Compute full atmospheric composition from geodetic coordinates and epoch.
338///
339/// Returns the complete NRLMSISE-00 output including:
340/// - Total mass density \[kg/m³\]
341/// - Number densities \[cm⁻³\] for 9 species: He, O, N₂, O₂, Ar, H, N, anomalous O
342/// - Exospheric and local temperatures \[K\]
343///
344/// This is the high-level API that takes pre-computed geodetic coordinates.
345/// For direct low-level access with explicit NRLMSISE-00 input parameters,
346/// use [`Nrlmsise00::calculate()`].
347pub fn msise00_density(
348    sw: Msise00DailyWeather,
349    lst_h: f64,
350    latitude_deg: f64,
351    longitude_deg: f64,
352    altitude_km: f64,
353    epoch: Epoch,
354    flags: Nrlmsise00Flags,
355) -> Result<Nrlmsise00Output, DynamicsError> {
356    let at_midnight = epoch.with_hms(0, 0, 0);
357    let ut_seconds = (epoch - at_midnight).to_seconds();
358
359    let input = Nrlmsise00Input {
360        day_of_year: at_midnight.day_of_year() as u32,
361        ut_seconds,
362        altitude_km,
363        latitude_deg,
364        longitude_deg,
365        local_solar_time_hours: lst_h,
366        f107_daily: sw.f107_daily_sfu,
367        f107_avg: sw.f107_avg_sfu,
368        ap_daily: sw.ap_daily,
369        ap_array: sw.ap_3hour_history,
370    };
371
372    Ok(calculate(&input, flags))
373}
374
375#[cfg(test)]
376mod tests {
377    use super::*;
378
379    #[test]
380    fn test_nrlmsise00_flags() {
381        let default_flags = Nrlmsise00Flags::default();
382        let default_switches = default_flags.to_switches();
383
384        // Spot-check standard switches mapping
385        assert_eq!(default_switches[0], 1.0);
386        assert_eq!(default_switches[9], 1.0);
387        assert_eq!(default_switches[1], 1.0);
388        assert_eq!(default_switches[2], 1.0);
389
390        // Customize some flags
391        let mut custom_flags = Nrlmsise00Flags {
392            geomagnetic: GeomagneticMode::StandardDailyAp,
393            f107_solar_flux: false,
394            time_independent: false,
395            annual_harmonics: false,
396            semiannual_harmonics: false,
397            diurnal_tides: false,
398            semidiurnal_tides: false,
399            terdiurnal_tides: false,
400            ut_and_longitude: false,
401            exospheric_temp_variations: false,
402            lower_boundary_temp_variations: false,
403            gradient_variations: false,
404            departures_from_diffusive_equilibrium: false,
405            lower_thermosphere_temp_variations: false,
406            upper_stratosphere_temp_variations: false,
407            boundary_density_variations: false,
408            lower_mesosphere_temp_variations: false,
409            turbopause_scale_height_variations: false,
410        };
411
412        let custom_switches = custom_flags.to_switches();
413        assert_eq!(custom_switches[0], 1.0); // Always 1.0
414        assert_eq!(custom_switches[9], 1.0); // GeomagneticMode::StandardDailyAp -> 1.0
415        assert_eq!(custom_switches[1], 0.0);
416        assert_eq!(custom_switches[2], 0.0);
417        assert_eq!(custom_switches[3], 0.0);
418        assert_eq!(custom_switches[5], 0.0);
419        assert_eq!(custom_switches[4], 0.0);
420        assert_eq!(custom_switches[6], 0.0);
421        assert_eq!(custom_switches[7], 0.0);
422        assert_eq!(custom_switches[8], 0.0);
423        assert_eq!(custom_switches[14], 0.0);
424        assert_eq!(custom_switches[10], 0.0);
425        assert_eq!(custom_switches[11], 0.0);
426        assert_eq!(custom_switches[12], 0.0);
427        assert_eq!(custom_switches[13], 0.0);
428        assert_eq!(custom_switches[16], 0.0);
429        assert_eq!(custom_switches[17], 0.0);
430        assert_eq!(custom_switches[19], 0.0);
431        assert_eq!(custom_switches[15], 0.0);
432        assert_eq!(custom_switches[18], 0.0);
433        assert_eq!(custom_switches[20], 0.0);
434        assert_eq!(custom_switches[21], 0.0);
435        assert_eq!(custom_switches[22], 0.0);
436        assert_eq!(custom_switches[23], 0.0);
437
438        // Test GeomagneticMode::Off
439        custom_flags.geomagnetic = GeomagneticMode::Off;
440        let switches_off = custom_flags.to_switches();
441        assert_eq!(switches_off[9], 0.0);
442    }
443}
nyx_space/dynamics/drag/nrlmsise00/mod.rs

nyx_space/dynamics/drag/nrlmsise00/
mod.rs
