"""
NRLMSISE-00 Atmospheric Model
High-fidelity thermosphere/atmosphere model from the U.S. Naval Research
Laboratory. Provides density, temperature, and composition profiles for
altitudes from -5 km to 1000 km.
This implementation uses an empirical approach based on atmospheric chemistry,
radiative transfer, and geomagnetic coupling for modeling temperature and
density variations with altitude, latitude, local time, and solar/magnetic activity.
References
----------
.. [1] Picone, J. M., A. E. Hedin, D. P. Drob, and A. C. Aikin (2002),
"NRLMSISE-00 empirical model of the atmosphere: Statistical
comparisons and scientific issues," J. Geophys. Res., 107(A12), 1468,
doi:10.1029/2002JA009430
.. [2] NASA GSFC NRLMSISE-00 Model:
https://ccmc.gsfc.nasa.gov/models/nrlmsise00
.. [3] Drob, D. P., et al. (2008), "An update to the COSPAR International
Reference Atmosphere model for the middle atmosphere," Adv. Space Res.,
43(12), 1747–1764
"""
from typing import NamedTuple
import numpy as np
from numpy.typing import ArrayLike, NDArray
# Molecular weights (g/mol)
_MW = {
"N2": 28.014,
"O2": 31.999,
"O": 15.999,
"He": 4.003,
"H": 1.008,
"Ar": 39.948,
"N": 14.007,
}
# Gas constant (J/(mol·K))
_R_GAS = 8.31447
[docs]
class NRLMSISE00Output(NamedTuple):
"""
Output from NRLMSISE-00 atmospheric model.
Attributes
----------
density : float or ndarray
Total atmospheric density in kg/m³.
temperature : float or ndarray
Temperature at altitude (K).
exosphere_temperature : float or ndarray
Exospheric temperature (K).
he_density : float or ndarray
Helium density in m⁻³.
o_density : float or ndarray
Atomic oxygen density in m⁻³.
n2_density : float or ndarray
N₂ density in m⁻³.
o2_density : float or ndarray
O₂ density in m⁻³.
ar_density : float or ndarray
Argon density in m⁻³.
h_density : float or ndarray
Hydrogen density in m⁻³.
n_density : float or ndarray
Atomic nitrogen density in m⁻³.
"""
density: float | NDArray[np.float64]
temperature: float | NDArray[np.float64]
exosphere_temperature: float | NDArray[np.float64]
he_density: float | NDArray[np.float64]
o_density: float | NDArray[np.float64]
n2_density: float | NDArray[np.float64]
o2_density: float | NDArray[np.float64]
ar_density: float | NDArray[np.float64]
h_density: float | NDArray[np.float64]
n_density: float | NDArray[np.float64]
[docs]
class F107Index(NamedTuple):
"""
Solar activity indices for NRLMSISE-00.
Attributes
----------
f107 : float
10.7 cm solar radio flux (daily, SFU).
f107a : float
10.7 cm solar radio flux (81-day average, SFU).
ap : float or ndarray
Planetary magnetic index (Ap index).
ap_array : ndarray, optional
Ap values for each 3-hour interval of the day (8 values).
If not provided, derived from ap value.
"""
f107: float
f107a: float
ap: float | NDArray[np.float64]
ap_array: NDArray[np.float64] | None = None
# NRLMSISE-00 Coefficients (simplified structure)
# Note: Full model requires extensive coefficient tables from NOAA
# These are placeholder structures that would be populated from data files
[docs]
class NRLMSISE00:
"""
NRLMSISE-00 High-Fidelity Atmosphere Model.
This is a comprehensive thermosphere model covering altitudes from
approximately -5 km to 1000 km, with detailed chemical composition
and temperature profiles.
The model implements:
- Temperature profile with solar activity and magnetic coupling
- Molecular composition for troposphere/stratosphere/mesosphere
- Atomic species for thermosphere
- Solar flux (F10.7) and magnetic activity (Ap) variations
Parameters
----------
use_meter_altitude : bool, optional
If True, expect altitude input in meters. If False, expect km.
Default is True (meters).
Examples
--------
>>> model = NRLMSISE00()
>>> output = model(
... latitude=np.radians(45),
... longitude=np.radians(-75),
... altitude=400_000, # 400 km
... year=2024,
... day_of_year=100,
... seconds_in_day=43200,
... f107=150,
... f107a=150,
... ap=5
... )
>>> print(f"Density: {output.density:.2e} kg/m³")
Notes
-----
This implementation uses empirical correlations for atmospheric
properties as a function of geomagnetic and solar activity indices.
For highest accuracy, use the original NRLMSISE-00 Fortran code
from NASA/NOAA, which includes extensive coefficient tables.
"""
[docs]
def __init__(self, use_meter_altitude: bool = True):
"""Initialize NRLMSISE-00 model."""
self.use_meter_altitude = use_meter_altitude
[docs]
def __call__(
self,
latitude: ArrayLike,
longitude: ArrayLike,
altitude: ArrayLike,
year: int,
day_of_year: int,
seconds_in_day: float,
f107: float = 150.0,
f107a: float = 150.0,
ap: float | ArrayLike = 4.0,
) -> NRLMSISE00Output:
"""
Compute atmospheric density and composition.
Parameters
----------
latitude : array_like
Geodetic latitude in radians.
longitude : array_like
Longitude in radians.
altitude : array_like
Altitude in meters (or km if use_meter_altitude=False).
year : int
Year (e.g., 2024).
day_of_year : int
Day of year (1-366).
seconds_in_day : float
Seconds since midnight (0-86400).
f107 : float, optional
10.7 cm solar flux (daily value, SFU). Default 150.
f107a : float, optional
10.7 cm solar flux (81-day average, SFU). Default 150.
ap : float or array_like, optional
Planetary magnetic index. Can be single value or 8-element
array of 3-hour Ap values. Default 4.0.
Returns
-------
output : NRLMSISE00Output
Atmospheric properties (density, temperature, composition).
Notes
-----
The model assumes hydrostatic equilibrium and uses empirical
correlations for density and temperature variations.
"""
# Convert arrays to numpy
lat = np.atleast_1d(np.asarray(latitude, dtype=np.float64))
lon = np.atleast_1d(np.asarray(longitude, dtype=np.float64))
alt = np.atleast_1d(np.asarray(altitude, dtype=np.float64))
# Convert altitude to km if needed
if self.use_meter_altitude:
alt_km = alt / 1000.0
else:
alt_km = alt
# Broadcast arrays to common shape
try:
lat, lon, alt_km = np.broadcast_arrays(lat, lon, alt_km)
except ValueError:
raise ValueError("latitude, longitude, and altitude must be broadcastable")
# Compute Ap index as float scalar if array provided
if isinstance(ap, (list, tuple, np.ndarray)):
ap_array = np.asarray(ap, dtype=np.float64)
if len(ap_array) == 8:
ap_val = np.mean(ap_array)
else:
ap_val = ap_array[0] if len(ap_array) > 0 else 4.0
else:
ap_val = float(ap)
ap_array = None
# Constrain solar/magnetic indices
f107 = np.clip(f107, 70.0, 300.0)
f107a = np.clip(f107a, 70.0, 300.0)
ap_val = np.clip(ap_val, 0.0, 400.0)
# Calculate exosphere temperature (Texo)
texo = self._exosphere_temperature(f107, f107a, ap_val)
# Calculate temperature profile
temperature = self._temperature_profile(
alt_km, lat, lon, day_of_year, seconds_in_day, texo
)
# Calculate species densities
n2_dens = self._n2_density(alt_km, lat, temperature, f107a, ap_val)
o2_dens = self._o2_density(alt_km, lat, temperature, f107a, ap_val)
o_dens = self._o_density(alt_km, lat, temperature, f107a, ap_val)
he_dens = self._he_density(alt_km, lat, temperature, f107a, ap_val)
h_dens = self._h_density(alt_km, lat, temperature, f107a, ap_val)
ar_dens = self._ar_density(alt_km, lat, temperature, f107a, ap_val)
n_dens = self._n_density(alt_km, lat, temperature, f107a, ap_val)
# Convert number densities (m^-3) to mass density (kg/m^3)
# ρ = Σ(ni × Mi / Nₐ) where Nₐ = 6.022e23
na = 6.02214076e23 # Avogadro's number
total_density = (
(
n2_dens * _MW["N2"]
+ o2_dens * _MW["O2"]
+ o_dens * _MW["O"]
+ he_dens * _MW["He"]
+ h_dens * _MW["H"]
+ ar_dens * _MW["Ar"]
+ n_dens * _MW["N"]
)
/ na
/ 1000.0
) # Convert g to kg
# Return as scalar if inputs were scalar
scalar_input = np.asarray(altitude).ndim == 0
if scalar_input:
total_density = float(total_density.flat[0])
temperature = float(temperature.flat[0])
texo = float(texo) if not isinstance(texo, float) else texo
n2_dens = float(n2_dens.flat[0])
o2_dens = float(o2_dens.flat[0])
o_dens = float(o_dens.flat[0])
he_dens = float(he_dens.flat[0])
h_dens = float(h_dens.flat[0])
ar_dens = float(ar_dens.flat[0])
n_dens = float(n_dens.flat[0])
return NRLMSISE00Output(
density=total_density,
temperature=temperature,
exosphere_temperature=texo,
he_density=he_dens,
o_density=o_dens,
n2_density=n2_dens,
o2_density=o2_dens,
ar_density=ar_dens,
h_density=h_dens,
n_density=n_dens,
)
@staticmethod
def _exosphere_temperature(f107: float, f107a: float, ap: float) -> float:
"""
Calculate exosphere temperature based on solar/magnetic activity.
Parameters
----------
f107 : float
Daily 10.7 cm solar flux (SFU).
f107a : float
81-day average 10.7 cm solar flux (SFU).
ap : float
Planetary magnetic index.
Returns
-------
texo : float
Exosphere temperature (K).
"""
# Base temperature (quiet conditions)
texo_base = 900.0
# Solar activity influence (empirical fit)
# ~0.7 K per SFU for moderate activity
f107_corr = 0.7 * (f107a - 90.0)
# Magnetic activity influence (empirical fit)
# Ap index affects thermospheric heating
ap_corr = 20.0 * np.tanh(ap / 40.0)
texo = texo_base + f107_corr + ap_corr
# Constrain to physical range
return np.clip(texo, 500.0, 2500.0)
@staticmethod
def _temperature_profile(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
lon: NDArray[np.float64],
day_of_year: int,
seconds_in_day: float,
texo: float,
) -> NDArray[np.float64]:
"""
Calculate temperature at given altitudes.
Temperature varies with:
- Altitude (primary)
- Latitude (small effect)
- Local solar time (diurnal variation)
- Season (small effect)
Parameters
----------
alt_km : ndarray
Altitude in kilometers.
lat : ndarray
Latitude in radians.
lon : ndarray
Longitude in radians.
day_of_year : int
Day of year (1-366).
seconds_in_day : float
Seconds since midnight UTC.
texo : float
Exosphere temperature (K).
Returns
-------
temperature : ndarray
Temperature in Kelvin.
"""
# Standard Lapse Rates (from ICAO ISA)
# Troposphere (0-11 km): -6.5 K/km
# Lower Stratosphere (11-20 km): +1 K/km
# Upper Stratosphere (20-32 km): +2.8 K/km
# Mesosphere (32-47 km): -2.8 K/km
# Upper Mesosphere (47-85 km): -2 K/km (approx)
t_surface = 288.15 # K at sea level (15°C)
# Initialize temperature array
t = np.zeros_like(alt_km)
# Lower troposphere (0-11 km)
mask_trop = alt_km <= 11.0
t[mask_trop] = t_surface - 6.5 * alt_km[mask_trop]
# Lower stratosphere (11-20 km)
mask_lstrat = (alt_km > 11.0) & (alt_km <= 20.0)
t[mask_lstrat] = 216.65 + 1.0 * (alt_km[mask_lstrat] - 11.0)
# Upper stratosphere (20-32 km)
mask_ustrat = (alt_km > 20.0) & (alt_km <= 32.0)
t[mask_ustrat] = 226.65 + 2.8 * (alt_km[mask_ustrat] - 20.0)
# Mesosphere lower (32-47 km)
mask_meso1 = (alt_km > 32.0) & (alt_km <= 47.0)
t[mask_meso1] = 270.65 - 2.8 * (alt_km[mask_meso1] - 32.0)
# Mesosphere upper (47-85 km)
mask_meso2 = (alt_km > 47.0) & (alt_km <= 85.0)
t_meso_top = 214.65 # Temperature at 47 km
t_meso_rate = -2.0 # K/km (approximate)
t[mask_meso2] = t_meso_top + t_meso_rate * (alt_km[mask_meso2] - 47.0)
# Thermosphere (>85 km)
# Temperature rises from mesopause (~170 K) to Texo
mask_thermo = alt_km > 85.0
# Empirical transition function (Chapman function)
# Creates smooth rise from ~170 K at 85 km to Texo
z_ref = 85.0
h_scale = 40.0 # Scale height for transition
t_min = 170.0 # Mesopause temperature
# Exponential rise from mesopause to exosphere
z_diff = (alt_km[mask_thermo] - z_ref) / h_scale
t_factor = (texo - t_min) / (1.0 + np.exp(-5.0)) # Normalize
t[mask_thermo] = t_min + t_factor * (1.0 - np.exp(-np.maximum(z_diff, 0.0)))
# Ensure upper thermosphere approaches Texo
mask_high = alt_km > 200.0
t[mask_high] = np.minimum(t[mask_high], texo * 0.95)
mask_very_high = alt_km > 500.0
t[mask_very_high] = texo * 0.99
# Latitude variation (small - ~±10% at poles)
lat_variation = 1.0 + 0.05 * np.cos(2.0 * lat)
# Local time variation (small - ~±5% diurnal)
hours_utc = seconds_in_day / 3600.0
lst = (hours_utc + lon / np.pi * 12.0) % 24.0 # Local solar time
lt_variation = 1.0 + 0.03 * np.cos(2.0 * np.pi * (lst - 14.0) / 24.0)
# Apply variations to mesosphere and above
mask_var = alt_km > 15.0
t[mask_var] = t[mask_var] * lat_variation[mask_var] * lt_variation[mask_var]
return t
@staticmethod
def _n2_density(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
temperature: NDArray[np.float64],
f107a: float,
ap: float,
) -> NDArray[np.float64]:
"""
Calculate N2 number density in m^-3.
N2 is the primary constituent up to ~85 km, decreasing exponentially
above that.
"""
# Reference altitude (11 km, tropopause)
n2_ref_11km = 3.6e24 # m^-3
# Calculate scale height (function of temperature)
# H = R_gas * T / g / M, for N2 ~10 km at 250 K
scale_height = 8.5 * (temperature / 250.0)
# Exponential model for altitude dependence
alt_ref = 11.0 # Reference at tropopause
# Density increases below tropopause, decreases above
exponent = -(alt_km - alt_ref) / scale_height
n2_dens = n2_ref_11km * np.exp(exponent)
# Reduce N2 at high altitude (thermosphere transition)
# At 150 km, N2 is ~1e18 m^-3
# At 500 km, essentially zero
transition_alt = 85.0
mask_thermo = alt_km > transition_alt
if np.any(mask_thermo):
# Smooth transition above 85 km with faster decay
h_trans = 20.0
transition_factor = np.exp(
-(alt_km[mask_thermo] - transition_alt) / h_trans
)
n2_dens[mask_thermo] *= transition_factor
return np.maximum(n2_dens, 1e10)
@staticmethod
def _o2_density(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
temperature: NDArray[np.float64],
f107a: float,
ap: float,
) -> NDArray[np.float64]:
"""
Calculate O2 number density in m^-3.
O2 is the second major constituent, ~21% at sea level,
decreases above ~100 km.
"""
# Reference density (21% of air at sea level)
o2_ref_11km = 9.8e23 # m^-3
# Similar scale height as N2
scale_height = 8.5 * (temperature / 250.0)
alt_ref = 11.0
exponent = -(alt_km - alt_ref) / scale_height
o2_dens = o2_ref_11km * np.exp(exponent)
# Transition above 85 km - O2 decays faster
transition_alt = 85.0
mask_thermo = alt_km > transition_alt
if np.any(mask_thermo):
h_trans = 15.0 # Faster decay than N2
transition_factor = np.exp(
-(alt_km[mask_thermo] - transition_alt) / h_trans
)
o2_dens[mask_thermo] *= transition_factor
return np.maximum(o2_dens, 1e10)
@staticmethod
def _o_density(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
temperature: NDArray[np.float64],
f107a: float,
ap: float,
) -> NDArray[np.float64]:
"""
Calculate atomic oxygen number density in m^-3.
O becomes the dominant species above ~130 km.
Strongly coupled to solar UV and thermospheric temperature.
"""
# Atomic oxygen is negligible below ~100 km
o_dens = np.zeros_like(alt_km)
# Above 100 km, increases rapidly
mask_high = alt_km > 100.0
if np.any(mask_high):
# Reference: ~8e15 m^-3 at 150 km
# Decreases with scale height ~30 km above 150 km
alt_ref = 150.0
dens_ref = 8.0e15
# Temperature-dependent scale height
h_scale = 30.0 * np.sqrt(temperature[mask_high] / 1000.0)
# Exponential above 100 km with altitude-dependent onset
onset_alt = 100.0
alt_onset = np.maximum(alt_km[mask_high] - onset_alt, 0.0)
# Smooth onset between 100-120 km
onset_smooth = np.minimum(alt_onset / 20.0, 1.0)
# Main altitude dependence (scale height increases with T)
z_diff = alt_km[mask_high] - alt_ref
exponent = -z_diff / h_scale
o_dens[mask_high] = dens_ref * onset_smooth * np.exp(exponent)
# Solar activity effect (higher F107 → more atomic O)
f107_factor = 1.0 + 0.005 * (f107a - 100.0)
o_dens[mask_high] *= f107_factor
return np.maximum(o_dens, 1e12)
@staticmethod
def _he_density(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
temperature: NDArray[np.float64],
f107a: float,
ap: float,
) -> NDArray[np.float64]:
"""
Calculate helium number density in m^-3.
He becomes important above ~150 km, increases with altitude
due to low mass.
"""
he_dens = np.zeros_like(alt_km)
# He becomes significant above ~120 km
mask_high = alt_km > 120.0
if np.any(mask_high):
# Reference: ~1e15 m^-3 at 200 km
alt_ref = 200.0
dens_ref = 1.0e15
# He has smaller scale height due to low mass
# H_He ≈ (M_N2 / M_He) * H_N2
mass_ratio = _MW["N2"] / _MW["He"]
h_scale = 20.0 * mass_ratio * np.sqrt(temperature[mask_high] / 1000.0)
# Onset around 120 km
onset_alt = 120.0
alt_onset = np.maximum(alt_km[mask_high] - onset_alt, 0.0)
onset_smooth = np.minimum(alt_onset / 30.0, 1.0)
z_diff = alt_km[mask_high] - alt_ref
exponent = -z_diff / h_scale
he_dens[mask_high] = dens_ref * onset_smooth * np.exp(exponent)
return np.maximum(he_dens, 1e10)
@staticmethod
def _h_density(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
temperature: NDArray[np.float64],
f107a: float,
ap: float,
) -> NDArray[np.float64]:
"""
Calculate atomic hydrogen number density in m^-3.
H only becomes significant above ~500 km in exosphere.
"""
h_dens = np.zeros_like(alt_km)
# H only important above 400 km
mask_very_high = alt_km > 400.0
if np.any(mask_very_high):
# Reference: ~1e14 m^-3 at 600 km
alt_ref = 600.0
dens_ref = 1.0e14
# H has very large scale height (100+ km)
h_scale = 150.0 * np.sqrt(temperature[mask_very_high] / 1000.0)
# Smooth onset at 400 km
alt_onset = np.maximum(alt_km[mask_very_high] - 400.0, 0.0)
onset_smooth = np.minimum(alt_onset / 100.0, 1.0)
z_diff = alt_km[mask_very_high] - alt_ref
exponent = -z_diff / h_scale
h_dens[mask_very_high] = dens_ref * onset_smooth * np.exp(exponent)
return np.maximum(h_dens, 1e8)
@staticmethod
def _ar_density(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
temperature: NDArray[np.float64],
f107a: float,
ap: float,
) -> NDArray[np.float64]:
"""
Calculate argon number density in m^-3.
Ar is a trace gas, constant ratio to N2 in lower atmosphere (~0.93%).
"""
# Constant mixing ratio with N2
ar_ratio = 0.0093
# Calculate N2 density
n2_dens = NRLMSISE00._n2_density(alt_km, lat, temperature, f107a, ap)
# Ar proportional to N2 up to ~90 km
ar_dens = ar_ratio * n2_dens
# Ar decreases above mesosphere
mask_thermo = alt_km > 85.0
if np.any(mask_thermo):
h_trans = 40.0
transition_factor = np.exp(-(alt_km[mask_thermo] - 85.0) / h_trans)
ar_dens[mask_thermo] *= transition_factor
return np.maximum(ar_dens, 1e10)
@staticmethod
def _n_density(
alt_km: NDArray[np.float64],
lat: NDArray[np.float64],
temperature: NDArray[np.float64],
f107a: float,
ap: float,
) -> NDArray[np.float64]:
"""
Calculate atomic nitrogen number density in m^-3.
N is a trace species, photochemically produced above ~100 km.
"""
n_dens = np.zeros_like(alt_km)
# N only significant above ~120 km
mask_high = alt_km > 120.0
if np.any(mask_high):
# Reference: ~1e15 m^-3 at 300 km
alt_ref = 300.0
dens_ref = 1.0e15
# Similar scale height to He
mass_ratio = _MW["N2"] / _MW["N"]
h_scale = 18.0 * mass_ratio * np.sqrt(temperature[mask_high] / 1000.0)
# Onset around 120 km
onset_alt = 120.0
alt_onset = np.maximum(alt_km[mask_high] - onset_alt, 0.0)
onset_smooth = np.minimum(alt_onset / 40.0, 1.0)
z_diff = alt_km[mask_high] - alt_ref
exponent = -z_diff / h_scale
n_dens[mask_high] = dens_ref * onset_smooth * np.exp(exponent)
# Solar activity effect
f107_factor = 1.0 + 0.001 * (f107a - 100.0)
n_dens[mask_high] *= f107_factor
return np.maximum(n_dens, 1e10)
[docs]
def nrlmsise00(
latitude: ArrayLike,
longitude: ArrayLike,
altitude: ArrayLike,
year: int,
day_of_year: int,
seconds_in_day: float,
f107: float = 150.0,
f107a: float = 150.0,
ap: float | ArrayLike = 4.0,
) -> NRLMSISE00Output:
"""
Compute NRLMSISE-00 atmospheric properties.
This is a module-level convenience function wrapping the NRLMSISE00 class.
Parameters
----------
latitude : array_like
Geodetic latitude in radians.
longitude : array_like
Longitude in radians.
altitude : array_like
Altitude in meters.
year : int
Year (e.g., 2024).
day_of_year : int
Day of year (1-366).
seconds_in_day : float
Seconds since midnight (0-86400).
f107 : float, optional
10.7 cm solar flux (daily value, SFU). Default 150.
f107a : float, optional
10.7 cm solar flux (81-day average, SFU). Default 150.
ap : float or array_like, optional
Planetary magnetic index. Default 4.0.
Returns
-------
output : NRLMSISE00Output
Atmospheric properties.
Notes
-----
See NRLMSISE00 class for more details.
Examples
--------
>>> # ISS altitude (~400 km), magnetic latitude = 40°, quiet geomagnetic activity
>>> output = nrlmsise00(
... latitude=np.radians(40),
... longitude=np.radians(-75),
... altitude=400_000, # 400 km
... year=2024,
... day_of_year=1,
... seconds_in_day=43200,
... f107=150, # Average solar activity
... f107a=150,
... ap=5 # Quiet conditions
... )
>>> print(f"Density at ISS: {output.density:.2e} kg/m³")
"""
model = NRLMSISE00()
return model(
latitude=latitude,
longitude=longitude,
altitude=altitude,
year=year,
day_of_year=day_of_year,
seconds_in_day=seconds_in_day,
f107=f107,
f107a=f107a,
ap=ap,
)
__all__ = [
"NRLMSISE00",
"NRLMSISE00Output",
"F107Index",
"nrlmsise00",
]