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, TimeScale};
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/// Defines all of the available flags in the NRLMSISE00 model.
114/// NOTE By default, Nyx will use the mean local solar time computation. Set mean_lst=false
115/// to use the apparent local solar time.
116#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, StaticType)]
117#[cfg_attr(feature = "python", pyclass(from_py_object, get_all, set_all))]
118pub struct Nrlmsise00Flags {
119    pub geomagnetic: GeomagneticMode,
120    pub f107_solar_flux: bool,
121    pub time_independent: bool,
122    pub annual_harmonics: bool,
123    pub semiannual_harmonics: bool,
124    pub diurnal_tides: bool,
125    pub semidiurnal_tides: bool,
126    pub terdiurnal_tides: bool,
127    pub ut_and_longitude: bool,
128    pub exospheric_temp_variations: bool,
129    pub lower_boundary_temp_variations: bool,
130    pub gradient_variations: bool,
131    pub departures_from_diffusive_equilibrium: bool,
132    pub lower_thermosphere_temp_variations: bool,
133    pub upper_stratosphere_temp_variations: bool,
134    pub boundary_density_variations: bool,
135    pub lower_mesosphere_temp_variations: bool,
136    pub turbopause_scale_height_variations: bool,
137    /// Determines if the mean local solar time should be used instead of the apparent solar time.
138    pub mean_lst: bool,
139}
140
141impl Default for Nrlmsise00Flags {
142    fn default() -> Self {
143        Self {
144            geomagnetic: GeomagneticMode::StandardDailyAp,
145            f107_solar_flux: true,
146            time_independent: true,
147            annual_harmonics: true,
148            semiannual_harmonics: true,
149            diurnal_tides: true,
150            semidiurnal_tides: true,
151            terdiurnal_tides: true,
152            ut_and_longitude: true,
153            exospheric_temp_variations: true,
154            lower_boundary_temp_variations: true,
155            gradient_variations: true,
156            departures_from_diffusive_equilibrium: true,
157            lower_thermosphere_temp_variations: true,
158            upper_stratosphere_temp_variations: true,
159            boundary_density_variations: true,
160            lower_mesosphere_temp_variations: true,
161            turbopause_scale_height_variations: true,
162            mean_lst: true,
163        }
164    }
165}
166
167impl Nrlmsise00Flags {
168    /// Compiles the high-level flags into the raw 24-element float array consumed by the kernel.
169    pub(crate) fn to_switches(self) -> [f64; 24] {
170        let mut sw = [1.0f64; 24];
171
172        // NOTE Unit selection is ALWAYS set to 1.0 because the calculation code does not even check it.
173
174        // Geomagnetic storm mode
175        sw[9] = match self.geomagnetic {
176            GeomagneticMode::Off => 0.0,
177            GeomagneticMode::StandardDailyAp => 1.0,
178            GeomagneticMode::ExtendedHistory57h => -1.0,
179        };
180
181        // Specific feature toggles
182        if !self.f107_solar_flux {
183            sw[1] = 0.0;
184        }
185        if !self.time_independent {
186            sw[2] = 0.0;
187        }
188        if !self.annual_harmonics {
189            sw[3] = 0.0;
190            sw[5] = 0.0;
191        }
192        if !self.semiannual_harmonics {
193            sw[4] = 0.0;
194            sw[6] = 0.0;
195        }
196        if !self.diurnal_tides {
197            sw[7] = 0.0;
198        }
199        if !self.semidiurnal_tides {
200            sw[8] = 0.0;
201        }
202        if !self.terdiurnal_tides {
203            sw[14] = 0.0;
204        }
205        if !self.ut_and_longitude {
206            sw[10] = 0.0;
207            sw[11] = 0.0;
208            sw[12] = 0.0;
209            sw[13] = 0.0;
210        }
211        if !self.exospheric_temp_variations {
212            sw[16] = 0.0;
213        }
214        if !self.lower_boundary_temp_variations {
215            sw[17] = 0.0;
216        }
217        if !self.gradient_variations {
218            sw[19] = 0.0;
219        }
220        if !self.departures_from_diffusive_equilibrium {
221            sw[15] = 0.0;
222        }
223        if !self.lower_thermosphere_temp_variations {
224            sw[18] = 0.0;
225        }
226        if !self.upper_stratosphere_temp_variations {
227            sw[20] = 0.0;
228        }
229        if !self.boundary_density_variations {
230            sw[21] = 0.0;
231        }
232        if !self.lower_mesosphere_temp_variations {
233            sw[22] = 0.0;
234        }
235        if !self.turbopause_scale_height_variations {
236            sw[23] = 0.0;
237        }
238
239        sw
240    }
241}
242
243#[cfg(feature = "python")]
244#[pymethods]
245impl Nrlmsise00Flags {
246    #[new]
247    #[pyo3(signature = (
248        geomagnetic = None,
249        f107_solar_flux = true,
250        time_independent = true,
251        annual_harmonics = true,
252        semiannual_harmonics = true,
253        diurnal_tides = true,
254        semidiurnal_tides = true,
255        terdiurnal_tides = true,
256        ut_and_longitude = true,
257        exospheric_temp_variations = true,
258        lower_boundary_temp_variations = true,
259        gradient_variations = true,
260        departures_from_diffusive_equilibrium = true,
261        lower_thermosphere_temp_variations = true,
262        upper_stratosphere_temp_variations = true,
263        boundary_density_variations = true,
264        lower_mesosphere_temp_variations = true,
265        turbopause_scale_height_variations = true,
266        mean_lst = true
267    ))]
268    #[allow(clippy::too_many_arguments)]
269    fn py_new(
270        geomagnetic: Option<GeomagneticMode>,
271        f107_solar_flux: bool,
272        time_independent: bool,
273        annual_harmonics: bool,
274        semiannual_harmonics: bool,
275        diurnal_tides: bool,
276        semidiurnal_tides: bool,
277        terdiurnal_tides: bool,
278        ut_and_longitude: bool,
279        exospheric_temp_variations: bool,
280        lower_boundary_temp_variations: bool,
281        gradient_variations: bool,
282        departures_from_diffusive_equilibrium: bool,
283        lower_thermosphere_temp_variations: bool,
284        upper_stratosphere_temp_variations: bool,
285        boundary_density_variations: bool,
286        lower_mesosphere_temp_variations: bool,
287        turbopause_scale_height_variations: bool,
288        mean_lst: bool,
289    ) -> Self {
290        Self {
291            geomagnetic: geomagnetic.unwrap_or(GeomagneticMode::StandardDailyAp),
292            f107_solar_flux,
293            time_independent,
294            annual_harmonics,
295            semiannual_harmonics,
296            diurnal_tides,
297            semidiurnal_tides,
298            terdiurnal_tides,
299            ut_and_longitude,
300            exospheric_temp_variations,
301            lower_boundary_temp_variations,
302            gradient_variations,
303            departures_from_diffusive_equilibrium,
304            lower_thermosphere_temp_variations,
305            upper_stratosphere_temp_variations,
306            boundary_density_variations,
307            lower_mesosphere_temp_variations,
308            turbopause_scale_height_variations,
309            mean_lst,
310        }
311    }
312
313    fn __repr__(&self) -> String {
314        format!("{:?}", self)
315    }
316
317    fn __str__(&self) -> String {
318        format!("{:?} @ {self:p}", self)
319    }
320}
321
322/// Compute full NRLMSISE-00 output for the given input parameters.
323///
324/// Returns temperatures and all species number densities.
325fn calculate(input: &Nrlmsise00Input, flags: Nrlmsise00Flags) -> Nrlmsise00Output {
326    let sw = flags.to_switches();
327    let (d, temp_exo, temp_alt) = model::compute(input, &sw);
328    // d[0..8]: He, O, N2, O2, Ar, total_mass(g/cm³), H, N, anomO
329    // Total mass density: d[5] is in g/cm³, convert to kg/m³ (* 1000)
330    Nrlmsise00Output {
331        temp_exo_k: temp_exo,
332        temp_alt_k: temp_alt,
333        density_he_per_cm3: d[0],
334        density_o_per_cm3: d[1],
335        density_n2_per_cm3: d[2],
336        density_o2_per_cm3: d[3],
337        density_ar_per_cm3: d[4],
338        density_h_per_cm3: d[6],
339        density_n_per_cm3: d[7],
340        density_anomalous_o_per_cm3: d[8],
341        // Convert g/cm^3 to kg/m^3: g->kg <=> 1e-3; cm^3 -> m^3 <-> 1e6 => 1e3
342        total_mass_density_kg_m3: d[5] * 1e3,
343    }
344}
345
346/// Compute full atmospheric composition from geodetic coordinates and epoch.
347///
348/// Returns the complete NRLMSISE-00 output including:
349/// - Total mass density \[kg/m³\]
350/// - Number densities \[cm⁻³\] for 9 species: He, O, N₂, O₂, Ar, H, N, anomalous O
351/// - Exospheric and local temperatures \[K\]
352///
353/// This is the high-level API that takes pre-computed geodetic coordinates.
354/// For direct low-level access with explicit NRLMSISE-00 input parameters,
355/// use [`Nrlmsise00::calculate()`].
356pub fn msise00_density(
357    sw: Msise00DailyWeather,
358    lst_h: f64,
359    latitude_deg: f64,
360    longitude_deg: f64,
361    altitude_km: f64,
362    mut epoch: Epoch,
363    flags: Nrlmsise00Flags,
364) -> Result<Nrlmsise00Output, DynamicsError> {
365    // Space weather is provided in UTC.
366    epoch = epoch.to_time_scale(TimeScale::UTC);
367    let at_midnight = epoch.with_hms_strict(0, 0, 0);
368    let ut_seconds = (epoch - at_midnight).to_seconds();
369
370    let input = Nrlmsise00Input {
371        day_of_year: at_midnight.day_of_year() as u32,
372        ut_seconds,
373        altitude_km,
374        latitude_deg,
375        longitude_deg,
376        local_solar_time_hours: lst_h,
377        f107_daily: sw.f107_daily_sfu,
378        f107_avg: sw.f107_avg_sfu,
379        ap_daily: sw.ap_daily,
380        ap_array: sw.ap_3hour_history,
381    };
382
383    Ok(calculate(&input, flags))
384}
385
386#[cfg(test)]
387mod tests {
388    use super::*;
389    use approx::assert_relative_eq;
390    use serde::Deserialize;
391    use std::fs;
392    use std::path::PathBuf;
393
394    #[test]
395    fn test_nrlmsise00_flags() {
396        let default_flags = Nrlmsise00Flags::default();
397        let default_switches = default_flags.to_switches();
398
399        // Spot-check standard switches mapping
400        assert_eq!(default_switches[0], 1.0);
401        assert_eq!(default_switches[9], 1.0);
402        assert_eq!(default_switches[1], 1.0);
403        assert_eq!(default_switches[2], 1.0);
404
405        // Customize some flags
406        let mut custom_flags = Nrlmsise00Flags {
407            geomagnetic: GeomagneticMode::StandardDailyAp,
408            f107_solar_flux: false,
409            time_independent: false,
410            annual_harmonics: false,
411            semiannual_harmonics: false,
412            diurnal_tides: false,
413            semidiurnal_tides: false,
414            terdiurnal_tides: false,
415            ut_and_longitude: false,
416            exospheric_temp_variations: false,
417            lower_boundary_temp_variations: false,
418            gradient_variations: false,
419            departures_from_diffusive_equilibrium: false,
420            lower_thermosphere_temp_variations: false,
421            upper_stratosphere_temp_variations: false,
422            boundary_density_variations: false,
423            lower_mesosphere_temp_variations: false,
424            turbopause_scale_height_variations: false,
425            mean_lst: false,
426        };
427
428        let custom_switches = custom_flags.to_switches();
429        assert_eq!(custom_switches[0], 1.0); // Always 1.0
430        assert_eq!(custom_switches[9], 1.0); // GeomagneticMode::StandardDailyAp -> 1.0
431        assert_eq!(custom_switches[1], 0.0);
432        assert_eq!(custom_switches[2], 0.0);
433        assert_eq!(custom_switches[3], 0.0);
434        assert_eq!(custom_switches[5], 0.0);
435        assert_eq!(custom_switches[4], 0.0);
436        assert_eq!(custom_switches[6], 0.0);
437        assert_eq!(custom_switches[7], 0.0);
438        assert_eq!(custom_switches[8], 0.0);
439        assert_eq!(custom_switches[14], 0.0);
440        assert_eq!(custom_switches[10], 0.0);
441        assert_eq!(custom_switches[11], 0.0);
442        assert_eq!(custom_switches[12], 0.0);
443        assert_eq!(custom_switches[13], 0.0);
444        assert_eq!(custom_switches[16], 0.0);
445        assert_eq!(custom_switches[17], 0.0);
446        assert_eq!(custom_switches[19], 0.0);
447        assert_eq!(custom_switches[15], 0.0);
448        assert_eq!(custom_switches[18], 0.0);
449        assert_eq!(custom_switches[20], 0.0);
450        assert_eq!(custom_switches[21], 0.0);
451        assert_eq!(custom_switches[22], 0.0);
452        assert_eq!(custom_switches[23], 0.0);
453
454        // Test GeomagneticMode::Off
455        custom_flags.geomagnetic = GeomagneticMode::Off;
456        let switches_off = custom_flags.to_switches();
457        assert_eq!(switches_off[9], 0.0);
458    }
459
460    #[derive(Deserialize)]
461    struct MsisTestCase {
462        altitude_km: f64,
463        latitude_deg: f64,
464        longitude_deg: f64,
465        day_of_year: u32,
466        ut_seconds: f64,
467        f107_daily: f64,
468        f107_avg: f64,
469        ap_array: [f64; 7],
470        is_storm: bool,
471        expected_total_density_kg_m3: f64,
472        expected_temperature_k: f64,
473    }
474
475    #[test]
476    fn pymsis_validation() {
477        let test_data: PathBuf = [
478            env!("CARGO_MANIFEST_DIR"),
479            "../data/03_tests/nrlmsise00_validation.json",
480        ]
481        .iter()
482        .collect();
483
484        let data = fs::read_to_string(test_data).expect("Failed to read validation JSON");
485        let test_cases: Vec<MsisTestCase> =
486            serde_json::from_str(&data).expect("Failed to deserialize test cases");
487
488        for (i, tc) in test_cases.iter().enumerate() {
489            let input = Nrlmsise00Input {
490                day_of_year: tc.day_of_year,
491                ut_seconds: tc.ut_seconds,
492                altitude_km: tc.altitude_km,
493                latitude_deg: tc.latitude_deg,
494                longitude_deg: tc.longitude_deg,
495                local_solar_time_hours: (tc.ut_seconds / 3600.0 + tc.longitude_deg / 15.0)
496                    .rem_euclid(24.0),
497                f107_daily: tc.f107_daily,
498                f107_avg: tc.f107_avg,
499                ap_daily: tc.ap_array[0],
500                ap_array: tc.ap_array,
501            };
502
503            // Standard switches: 57-hour history enabled
504            let sw = Nrlmsise00Flags {
505                geomagnetic: if tc.is_storm {
506                    GeomagneticMode::ExtendedHistory57h
507                } else {
508                    GeomagneticMode::StandardDailyAp
509                },
510                ..Default::default()
511            };
512
513            let output = calculate(&input, sw);
514
515            let total_density_kg_m3 = output.total_mass_density_kg_m3;
516            let t_alt = output.temp_alt_k;
517
518            // Verify temperature at altitude with a 0.5% relative tolerance
519            println!(
520                "[storm={}] Temperature #{i}: alt={} km, lat={} deg. Rust: {t_alt}, Pymsis: {}",
521                tc.is_storm, tc.altitude_km, tc.latitude_deg, tc.expected_temperature_k
522            );
523            assert_relative_eq!(
524                t_alt,
525                tc.expected_temperature_k,
526                max_relative = 0.001,
527                epsilon = 1e-5
528            );
529
530            // Verify total mass density with a 1.0% relative tolerance
531            // (Density integration magnifies the floating point differences in the exponential term)
532            println!(
533                "[storm={}] Density #{i}: alt={} km. Rust: {total_density_kg_m3:.6e}, Pymsis: {:.6e}",
534                tc.is_storm, tc.altitude_km, tc.expected_total_density_kg_m3
535            );
536            assert_relative_eq!(
537                total_density_kg_m3,
538                tc.expected_total_density_kg_m3,
539                max_relative = 0.004,
540                epsilon = 1e-18,
541            );
542        }
543    }
544}
nyx_space/dynamics/drag/nrlmsise00/mod.rs

nyx_space/dynamics/drag/nrlmsise00/
mod.rs
