"""
Barometric thermosphere model.
**This is not NRLMSISE-00.** The module was named ``nrlmsise00`` and described
itself as a high-fidelity NRL model until gh-79. NRLMSISE-00 requires harmonic
coefficient tables from NOAA which this library does not distribute; what is
implemented here is a set of per-species exponential profiles with
temperature-dependent scale heights, driven by F10.7 and Ap, each clamped to a
floor.
It is usable above roughly 200 km, where it agrees with published NRLMSISE-00
densities to within a factor of about two. Below that it is wrong by up to
50x, and :func:`pytcl.atmosphere.us_standard_atmosphere_1976` should be used
instead. :class:`SimplifiedThermosphere` carries the measured comparison.
The references below describe the model this one approximates, not the one
implemented here. They are retained because the temperature parameterisation
and species set follow their structure.
References
----------
- 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
- NASA GSFC NRLMSISE-00 Model:
https://ccmc.gsfc.nasa.gov/models/nrlmsise00
- 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 ThermosphereState(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
# This module does NOT implement NRLMSISE-00. That model needs extensive
# harmonic coefficient tables from NOAA which are not distributed here. What
# follows is a barometric model: per-species exponential profiles with
# temperature-dependent scale heights, each clamped to a floor. It was named
# nrlmsise00 and documented as "a comprehensive thermosphere model", which is
# how it came to return a sea-level density 44% below the true value with
# nothing to warn a caller (gh-79).
[docs]
class SimplifiedThermosphere:
"""
Barometric thermosphere model with solar and geomagnetic coupling.
**This is not NRLMSISE-00.** It was named ``NRLMSISE00`` and described as
"a comprehensive thermosphere model" until gh-79; the real model requires
harmonic coefficient tables from NOAA that this library does not ship.
What this computes is a set of per-species exponential profiles with
temperature-dependent scale heights, driven by F10.7 and Ap, with each
species density clamped to a floor.
Use it above roughly 200 km. Below that it is wrong, in places by more
than an order of magnitude:
======== ============== ================ ==========
Altitude This model US Standard 1976 Ratio
======== ============== ================ ==========
0 km 0.682 kg/m^3 1.225 kg/m^3 0.56
10 km 0.253 0.414 0.61
30 km 0.0261 0.0184 1.42
50 km 1.044e-3 1.027e-3 1.02
80 km 3.69e-7 1.85e-5 **0.02**
======== ============== ================ ==========
For altitudes below about 86 km use
:func:`pytcl.atmosphere.us_standard_atmosphere_1976`, which is validated
against the standard. In the thermosphere this model is reasonable: at
400 km with F10.7 = 150 it gives 2.9e-12 kg/m^3, against a published
NRLMSISE-00 range of roughly 2-4e-12.
The species floors are visible in the output. ``h_density`` returns its
floor of 1e8 m^-3 from 0 to 400 km, and ``he_density`` its floor of
1e10 m^-3 at 0 and 100 km, because the computed profile falls below the
clamp. Those are not measurements of anything.
The model implements:
- Temperature profile with solar activity and magnetic coupling
- Per-species barometric profiles with floors
- 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 = SimplifiedThermosphere()
>>> 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³")
Density: 2.85e-12 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 the simplified thermosphere 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,
) -> ThermosphereState:
"""
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 : ThermosphereState
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 ThermosphereState(
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 = SimplifiedThermosphere._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 simplified_thermosphere(
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,
) -> ThermosphereState:
"""
Compute NRLMSISE-00 atmospheric properties.
This is a module-level convenience function wrapping the SimplifiedThermosphere 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 : ThermosphereState
Atmospheric properties.
Notes
-----
See SimplifiedThermosphere class for more details.
Examples
--------
>>> # ISS altitude (~400 km), magnetic latitude = 40°, quiet geomagnetic activity
>>> output = simplified_thermosphere(
... 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³")
Density at ISS: 2.89e-12 kg/m³
"""
model = SimplifiedThermosphere()
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__ = [
"SimplifiedThermosphere",
"ThermosphereState",
"F107Index",
"simplified_thermosphere",
]