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;
42
43pub mod coefficients;
44mod model;
45
46/// Full output of the NRLMSISE-00 model.
47///
48/// Includes temperatures and all species number densities.
49#[derive(Debug, Clone)]
50pub struct Nrlmsise00Output {
51 /// Exospheric temperature [K].
52 pub temp_exo_k: f64,
53 /// Temperature at altitude [K].
54 pub temp_alt_k: f64,
55 /// He number density [cm⁻³].
56 pub density_he_per_cm3: f64,
57 /// O number density [cm⁻³].
58 pub density_o_per_cm3: f64,
59 /// N₂ number density [cm⁻³].
60 pub density_n2_per_cm3: f64,
61 /// O₂ number density [cm⁻³].
62 pub density_o2_per_cm3: f64,
63 /// Ar number density [cm⁻³].
64 pub density_ar_per_cm3: f64,
65 /// H number density [cm⁻³].
66 pub density_h_per_cm3: f64,
67 /// N number density [cm⁻³].
68 pub density_n_per_cm3: f64,
69 /// Anomalous oxygen number density [cm⁻³].
70 pub density_anomalous_o_per_cm3: f64,
71 /// Total mass density [kg/m³].
72 pub total_mass_density_kg_m3: f64,
73}
74
75/// Input parameters for a single NRLMSISE-00 evaluation.
76#[derive(Debug, Clone)]
77pub struct Nrlmsise00Input {
78 /// Day of year [1-366].
79 pub day_of_year: u32,
80 /// Universal time [seconds since midnight].
81 pub ut_seconds: f64,
82 /// Geodetic altitude [km].
83 pub altitude_km: f64,
84 /// Geodetic latitude [degrees, -90 to 90].
85 pub latitude_deg: f64,
86 /// Geodetic longitude [degrees, 0 to 360 or -180 to 180].
87 pub longitude_deg: f64,
88 /// Local apparent solar time [hours, 0-24].
89 pub local_solar_time_hours: f64,
90 /// Previous day's F10.7 [SFU].
91 pub f107_daily: f64,
92 /// 81-day centered average F10.7 [SFU].
93 pub f107_avg: f64,
94 /// Daily Ap index.
95 pub ap_daily: f64,
96 /// 7-element Ap array for magnetic activity variations.
97 pub ap_array: [f64; 7],
98}
99
100/// Compute full NRLMSISE-00 output for the given input parameters.
101///
102/// Returns temperatures and all species number densities.
103fn calculate(input: &Nrlmsise00Input) -> Nrlmsise00Output {
104 let (d, temp_exo, temp_alt) = model::compute(input);
105 // d[0..8]: He, O, N2, O2, Ar, total_mass(g/cm³), H, N, anomO
106 // Total mass density: d[5] is in g/cm³, convert to kg/m³ (* 1000)
107 Nrlmsise00Output {
108 temp_exo_k: temp_exo,
109 temp_alt_k: temp_alt,
110 density_he_per_cm3: d[0],
111 density_o_per_cm3: d[1],
112 density_n2_per_cm3: d[2],
113 density_o2_per_cm3: d[3],
114 density_ar_per_cm3: d[4],
115 density_h_per_cm3: d[6],
116 density_n_per_cm3: d[7],
117 density_anomalous_o_per_cm3: d[8],
118 // Convert g/cm^3 to kg/m^3: g->kg <=> 1e-3; cm^3 -> m^3 <-> 1e6 => 1e3
119 total_mass_density_kg_m3: d[5] * 1e3,
120 }
121}
122
123/// Compute full atmospheric composition from geodetic coordinates and epoch.
124///
125/// Returns the complete NRLMSISE-00 output including:
126/// - Total mass density \[kg/m³\]
127/// - Number densities \[cm⁻³\] for 9 species: He, O, N₂, O₂, Ar, H, N, anomalous O
128/// - Exospheric and local temperatures \[K\]
129///
130/// This is the high-level API that takes pre-computed geodetic coordinates.
131/// For direct low-level access with explicit NRLMSISE-00 input parameters,
132/// use [`Nrlmsise00::calculate()`].
133pub fn msise00_density(
134 sw: Msise00DailyWeather,
135 lst_h: f64,
136 latitude_deg: f64,
137 longitude_deg: f64,
138 altitude_km: f64,
139 epoch: Epoch,
140) -> Result<Nrlmsise00Output, DynamicsError> {
141 let at_midnight = epoch.with_hms(0, 0, 0);
142 let ut_seconds = (epoch - at_midnight).to_seconds();
143
144 let input = Nrlmsise00Input {
145 day_of_year: at_midnight.day_of_year() as u32,
146 ut_seconds,
147 altitude_km,
148 latitude_deg,
149 longitude_deg,
150 local_solar_time_hours: lst_h,
151 f107_daily: sw.f107_daily_sfu,
152 f107_avg: sw.f107_avg_sfu,
153 ap_daily: sw.ap_daily,
154 ap_array: sw.ap_3hour_history,
155 };
156
157 Ok(calculate(&input))
158}