{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Atmospheric Density with NRLMSISE-00\n",
"\n",
"The [NRLMSISE-00](https://en.wikipedia.org/wiki/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).\n",
"\n",
"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:\n",
"\n",
"- **Orbit prediction**: Drag is the largest non-gravitational perturbation for LEO satellites\n",
"- **Satellite lifetime estimation**: Determines how long a satellite remains in orbit\n",
"- **Conjunction screening**: Accurate propagation is needed for collision avoidance\n",
"- **Re-entry prediction**: Density uncertainty dominates re-entry timing errors\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Basic Density Query\n",
"\n",
"The `sk.nrlmsise00` function returns a tuple of `(density, temperature)` where density is in kg/m^3 and temperature is in Kelvin.\n",
"\n",
"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.\n",
"\n",
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Query density at 400 km altitude (typical ISS orbit)\n",
"# at latitude 40 deg N, longitude -75 deg W\n",
"# at a specific time\n",
"tm = sk.time(2024, 6, 15, 12, 0, 0)\n",
"\n",
"rho, temp = sk.nrlmsise00(\n",
" 400.0,\n",
" latitude_deg=40.0,\n",
" longitude_deg=-75.0,\n",
" time=tm,\n",
")\n",
"\n",
"print(\"Altitude: 400 km\")\n",
"print(f\"Density: {rho:.4e} kg/m^3\")\n",
"print(f\"Temperature: {temp:.1f} K\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Without specifying location or time, the model uses defaults\n",
"# (lat=0, lon=0, moderate solar activity)\n",
"rho_default, temp_default = sk.nrlmsise00(400.0, use_spaceweather=False)\n",
"\n",
"print(f\"Default density at 400 km: {rho_default:.4e} kg/m^3\")\n",
"print(f\"Default temperature at 400 km: {temp_default:.1f} K\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Density vs. Altitude\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute density at altitudes from 100 to 1000 km\n",
"altitudes = np.linspace(100, 1000, 200)\n",
"tm = sk.time(2024, 6, 15, 12, 0, 0)\n",
"\n",
"densities = np.array([\n",
" sk.nrlmsise00(alt, latitude_deg=40.0, longitude_deg=-75.0, time=tm)[0]\n",
" for alt in altitudes\n",
"])\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"ax.semilogy(altitudes, densities)\n",
"ax.set_xlabel(\"Altitude (km)\", fontsize=13)\n",
"ax.set_ylabel(\"Density (kg/m$^3$)\", fontsize=13)\n",
"ax.set_title(\"NRLMSISE-00 Atmospheric Density vs. Altitude\", fontsize=14)\n",
"ax.grid(True, which=\"both\", alpha=0.3)\n",
"\n",
"# Annotate key altitudes\n",
"for label, h in [(\"ISS (~400 km)\", 400), (\"Hubble (~540 km)\", 540), (\"Starlink (~550 km)\", 550)]:\n",
" rho_ann = sk.nrlmsise00(float(h), latitude_deg=40.0, longitude_deg=-75.0, time=tm)[0]\n",
" ax.axhline(rho_ann, color=\"gray\", linestyle=\"--\", alpha=0.4)\n",
" ax.annotate(label, xy=(h, rho_ann), fontsize=10,\n",
" xytext=(h + 50, rho_ann * 3), arrowprops=dict(arrowstyle=\"->\", color=\"gray\"))\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solar Activity Effects\n",
"\n",
"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:\n",
"\n",
"- **Solar minimum**: F10.7 ~ 70 SFU\n",
"- **Solar maximum**: F10.7 ~ 200+ SFU\n",
"\n",
"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.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Solar minimum (late 2019) vs. solar maximum (late 2024)\n",
"tm_solar_min = sk.time(2019, 12, 15, 12, 0, 0)\n",
"tm_solar_max = sk.time(2024, 10, 15, 12, 0, 0)\n",
"\n",
"altitudes = np.linspace(150, 800, 150)\n",
"\n",
"rho_min = np.array([\n",
" sk.nrlmsise00(alt, latitude_deg=40.0, longitude_deg=0.0, time=tm_solar_min)[0]\n",
" for alt in altitudes\n",
"])\n",
"\n",
"rho_max = np.array([\n",
" sk.nrlmsise00(alt, latitude_deg=40.0, longitude_deg=0.0, time=tm_solar_max)[0]\n",
" for alt in altitudes\n",
"])\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))\n",
"\n",
"# Left panel: density profiles\n",
"ax1.semilogy(altitudes, rho_min, label=\"Solar Min (Dec 2019)\", color=\"steelblue\")\n",
"ax1.semilogy(altitudes, rho_max, label=\"Solar Max (Oct 2024)\", color=\"firebrick\")\n",
"ax1.set_xlabel(\"Altitude (km)\", fontsize=13)\n",
"ax1.set_ylabel(\"Density (kg/m$^3$)\", fontsize=13)\n",
"ax1.set_title(\"Density: Solar Min vs. Max\", fontsize=14)\n",
"ax1.legend(fontsize=11)\n",
"ax1.grid(True, which=\"both\", alpha=0.3)\n",
"\n",
"# Right panel: ratio of solar max to solar min density\n",
"ratio = rho_max / rho_min\n",
"ax2.plot(altitudes, ratio, color=\"darkgreen\")\n",
"ax2.set_xlabel(\"Altitude (km)\", fontsize=13)\n",
"ax2.set_ylabel(\"Density Ratio (Max / Min)\", fontsize=13)\n",
"ax2.set_title(\"Solar Max / Solar Min Density Ratio\", fontsize=14)\n",
"ax2.grid(True, alpha=0.3)\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Practical Impact: Orbital Decay at Different Altitudes\n",
"\n",
"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.\n",
"\n",
"We use satkit's high-precision orbit propagator with the `satproperties` object to set the satellite's ballistic coefficient (Cd * A / m)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Earth parameters\n",
"R_EARTH = 6378.137e3 # meters\n",
"MU_EARTH = 3.986004418e14 # m^3/s^2\n",
"\n",
"def circular_orbit_state(alt_km):\n",
" \"\"\"Return GCRF state vector [x, y, z, vx, vy, vz] for a circular orbit at given altitude.\"\"\"\n",
" r = R_EARTH + alt_km * 1e3\n",
" v = np.sqrt(MU_EARTH / r)\n",
" # Place satellite at (r, 0, 0) with velocity (0, v, 0) for equatorial circular orbit\n",
" return np.array([r, 0.0, 0.0, 0.0, v, 0.0])\n",
"\n",
"def semi_major_axis(state):\n",
" \"\"\"Compute semi-major axis from Cartesian state vector.\"\"\"\n",
" r = np.linalg.norm(state[0:3])\n",
" v = np.linalg.norm(state[3:6])\n",
" return 1.0 / (2.0 / r - v**2 / MU_EARTH)\n",
"\n",
"\n",
"# Set up two satellites at 300 km and 500 km\n",
"# Typical ballistic coefficient: Cd*A/m ~ 0.01 m^2/kg\n",
"sp = sk.satproperties(cdaoverm=0.01)\n",
"\n",
"t0 = sk.time(2024, 6, 15, 0, 0, 0)\n",
"duration_days = 2.0\n",
"t_end = t0 + sk.duration(days=duration_days)\n",
"\n",
"settings = sk.propsettings(abs_error=1e-8, rel_error=1e-8)\n",
"settings.precompute_terms(t0, t_end)\n",
"\n",
"altitudes_prop = [300, 500]\n",
"colors = [\"firebrick\", \"steelblue\"]\n",
"\n",
"fig, ax = plt.subplots(figsize=(9, 5))\n",
"\n",
"for alt, color in zip(altitudes_prop, colors):\n",
" state0 = circular_orbit_state(alt)\n",
"\n",
" # Propagate with drag\n",
" result = sk.propagate(\n",
" state0, t0, end=t_end,\n",
" propsettings=settings,\n",
" satproperties=sp,\n",
" )\n",
"\n",
" # Sample the trajectory at regular intervals\n",
" n_samples = 500\n",
" times = [t0 + sk.duration(days=d) for d in np.linspace(0, duration_days, n_samples)]\n",
" sma = np.array([semi_major_axis(result.interp(t)) for t in times])\n",
" sma_alt = (sma - R_EARTH) / 1e3 # Convert to altitude in km\n",
" days = np.linspace(0, duration_days, n_samples)\n",
"\n",
" ax.plot(days, sma_alt, label=f\"{alt} km initial altitude\", color=color)\n",
"\n",
"ax.set_xlabel(\"Time (days)\", fontsize=13)\n",
"ax.set_ylabel(\"Semi-Major Axis Altitude (km)\", fontsize=13)\n",
"ax.set_title(\"Orbital Decay Due to Atmospheric Drag\", fontsize=14)\n",
"ax.legend(fontsize=11)\n",
"ax.grid(True, alpha=0.3)\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}