Atmospheric Density with NRLMSISE-00¶
The NRLMSISE-00 (Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar Exosphere) model is the standard empirical model of Earth's upper atmosphere. It provides neutral temperature and density profiles from the ground to the exosphere (~1000 km).
Atmospheric density is the dominant source of orbital drag for satellites in low Earth orbit (LEO). Even though the atmosphere above 200 km is extremely tenuous, the high orbital velocities (~7.5 km/s) produce significant drag forces that cause orbital decay. Accurate density estimates are essential for:
- Orbit prediction: Drag is the largest non-gravitational perturbation for LEO satellites
- Satellite lifetime estimation: Determines how long a satellite remains in orbit
- Conjunction screening: Accurate propagation is needed for collision avoidance
- Re-entry prediction: Density uncertainty dominates re-entry timing errors
This tutorial demonstrates how to use satkit's nrlmsise00 function to query atmospheric density and explore how it varies with altitude, geographic location, and solar activity.
import satkit as sk
import numpy as np
import matplotlib.pyplot as plt
import scienceplots # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']
Basic Density Query¶
The sk.nrlmsise00 function returns a tuple of (density, temperature) where density is in kg/m^3 and temperature is in Kelvin.
The function takes altitude in kilometers as a positional argument, with optional keyword arguments for latitude, longitude, time, and whether to use space weather data.
When use_spaceweather=True (the default), the model uses historical or forecast F10.7 solar flux and geomagnetic Ap indices from satkit's built-in space weather database. When time or space weather data are not provided, it falls back to moderate default values (F10.7 = 150, Ap = 4).
# Query density at 400 km altitude (typical ISS orbit)
# at latitude 40 deg N, longitude -75 deg W
# at a specific time
tm = sk.time(2024, 6, 15, 12, 0, 0)
rho, temp = sk.nrlmsise00(
400.0,
latitude_deg=40.0,
longitude_deg=-75.0,
time=tm,
)
print("Altitude: 400 km")
print(f"Density: {rho:.4e} kg/m^3")
print(f"Temperature: {temp:.1f} K")
Altitude: 400 km Density: 4.6304e-12 kg/m^3 Temperature: 1180.1 K
# Without specifying location or time, the model uses defaults
# (lat=0, lon=0, moderate solar activity)
rho_default, temp_default = sk.nrlmsise00(400.0, use_spaceweather=False)
print(f"Default density at 400 km: {rho_default:.4e} kg/m^3")
print(f"Default temperature at 400 km: {temp_default:.1f} K")
Default density at 400 km: 2.6568e-12 kg/m^3 Default temperature at 400 km: 942.8 K
Density vs. Altitude¶
Atmospheric density drops roughly exponentially with altitude. The plot below shows density from 100 km (the Karman line) to 1000 km on a logarithmic scale. Note that the density spans many orders of magnitude across this range.
# Compute density at altitudes from 100 to 1000 km
altitudes = np.linspace(100, 1000, 200)
tm = sk.time(2024, 6, 15, 12, 0, 0)
densities = np.array([
sk.nrlmsise00(alt, latitude_deg=40.0, longitude_deg=-75.0, time=tm)[0]
for alt in altitudes
])
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogy(altitudes, densities)
ax.set_xlabel("Altitude (km)", fontsize=13)
ax.set_ylabel("Density (kg/m$^3$)", fontsize=13)
ax.set_title("NRLMSISE-00 Atmospheric Density vs. Altitude", fontsize=14)
ax.grid(True, which="both", alpha=0.3)
# Annotate key altitudes
for label, h in [("ISS (~400 km)", 400), ("Hubble (~540 km)", 540), ("Starlink (~550 km)", 550)]:
rho_ann = sk.nrlmsise00(float(h), latitude_deg=40.0, longitude_deg=-75.0, time=tm)[0]
ax.axhline(rho_ann, color="gray", linestyle="--", alpha=0.4)
ax.annotate(label, xy=(h, rho_ann), fontsize=10,
xytext=(h + 50, rho_ann * 3), arrowprops=dict(arrowstyle="->", color="gray"))
plt.tight_layout()
plt.show()
Solar Activity Effects¶
The upper atmosphere expands dramatically during periods of high solar activity. The F10.7 solar radio flux (measured in solar flux units, SFU) is the primary driver:
- Solar minimum: F10.7 ~ 70 SFU
- Solar maximum: F10.7 ~ 200+ SFU
At 400 km, density can vary by an order of magnitude between solar minimum and solar maximum. This makes solar activity forecasting critical for orbit prediction.
We demonstrate this by querying density at two different times corresponding to recent solar minimum and solar maximum conditions. With use_spaceweather=True, the model automatically looks up the historical F10.7 and Ap values for the given time.
# Solar minimum (late 2019) vs. solar maximum (late 2024)
tm_solar_min = sk.time(2019, 12, 15, 12, 0, 0)
tm_solar_max = sk.time(2024, 10, 15, 12, 0, 0)
altitudes = np.linspace(150, 800, 150)
rho_min = np.array([
sk.nrlmsise00(alt, latitude_deg=40.0, longitude_deg=0.0, time=tm_solar_min)[0]
for alt in altitudes
])
rho_max = np.array([
sk.nrlmsise00(alt, latitude_deg=40.0, longitude_deg=0.0, time=tm_solar_max)[0]
for alt in altitudes
])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
# Left panel: density profiles
ax1.semilogy(altitudes, rho_min, label="Solar Min (Dec 2019)", color="steelblue")
ax1.semilogy(altitudes, rho_max, label="Solar Max (Oct 2024)", color="firebrick")
ax1.set_xlabel("Altitude (km)", fontsize=13)
ax1.set_ylabel("Density (kg/m$^3$)", fontsize=13)
ax1.set_title("Density: Solar Min vs. Max", fontsize=14)
ax1.legend(fontsize=11)
ax1.grid(True, which="both", alpha=0.3)
# Right panel: ratio of solar max to solar min density
ratio = rho_max / rho_min
ax2.plot(altitudes, ratio, color="darkgreen")
ax2.set_xlabel("Altitude (km)", fontsize=13)
ax2.set_ylabel("Density Ratio (Max / Min)", fontsize=13)
ax2.set_title("Solar Max / Solar Min Density Ratio", fontsize=14)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Practical Impact: Orbital Decay at Different Altitudes¶
To see the practical effect of atmospheric drag, we propagate two LEO satellites at different altitudes and compare how their semi-major axes evolve over time. The higher satellite experiences far less drag due to the exponentially lower density.
We use satkit's high-precision orbit propagator with the satproperties object to set the satellite's ballistic coefficient (Cd * A / m).
# Earth parameters
R_EARTH = 6378.137e3 # meters
MU_EARTH = 3.986004418e14 # m^3/s^2
def circular_orbit_state(alt_km):
"""Return GCRF state vector [x, y, z, vx, vy, vz] for a circular orbit at given altitude."""
r = R_EARTH + alt_km * 1e3
v = np.sqrt(MU_EARTH / r)
# Place satellite at (r, 0, 0) with velocity (0, v, 0) for equatorial circular orbit
return np.array([r, 0.0, 0.0, 0.0, v, 0.0])
def semi_major_axis(state):
"""Compute semi-major axis from Cartesian state vector."""
r = np.linalg.norm(state[0:3])
v = np.linalg.norm(state[3:6])
return 1.0 / (2.0 / r - v**2 / MU_EARTH)
# Set up two satellites at 300 km and 500 km
# Typical ballistic coefficient: Cd*A/m ~ 0.01 m^2/kg
sp = sk.satproperties(cdaoverm=0.01)
t0 = sk.time(2024, 6, 15, 0, 0, 0)
duration_days = 2.0
t_end = t0 + sk.duration(days=duration_days)
settings = sk.propsettings(abs_error=1e-8, rel_error=1e-8)
settings.precompute_terms(t0, t_end)
altitudes_prop = [300, 500]
colors = ["firebrick", "steelblue"]
fig, ax = plt.subplots(figsize=(9, 5))
for alt, color in zip(altitudes_prop, colors):
state0 = circular_orbit_state(alt)
# Propagate with drag
result = sk.propagate(
state0, t0, end=t_end,
propsettings=settings,
satproperties=sp,
)
# Sample the trajectory at regular intervals
n_samples = 500
times = [t0 + sk.duration(days=d) for d in np.linspace(0, duration_days, n_samples)]
sma = np.array([semi_major_axis(result.interp(t)) for t in times])
sma_alt = (sma - R_EARTH) / 1e3 # Convert to altitude in km
days = np.linspace(0, duration_days, n_samples)
ax.plot(days, sma_alt, label=f"{alt} km initial altitude", color=color)
ax.set_xlabel("Time (days)", fontsize=13)
ax.set_ylabel("Semi-Major Axis Altitude (km)", fontsize=13)
ax.set_title("Orbital Decay Due to Atmospheric Drag", fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
The satellite at 300 km experiences substantially more drag than the one at 500 km, owing to the exponential density profile. This is why most operational LEO satellites orbit above 400 km -- lower altitudes require frequent orbit-raising maneuvers to maintain their orbit.