{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Interesting Plots\n",
"\n",
"A collection of visualizations demonstrating satkit's atmospheric density models, gravity field, and coordinate frame transformations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Air Density as a Function of the Solar Cycle\n",
"\n",
"Atmospheric density at orbital altitudes varies by more than an order of magnitude over the ~11-year solar cycle. During solar maximum, increased UV and EUV radiation heats the upper atmosphere, causing it to expand and increasing drag on satellites. This effect dominates orbit lifetime predictions for low-Earth orbit spacecraft.\n",
"\n",
"The plot below uses the NRLMSISE-00 density model at 400 km and 500 km altitude, showing density variations from 1995 to 2022 (spanning roughly two solar cycles)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\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']\n",
"import numpy as np\n",
"import math as m\n",
"\n",
"start = sk.time(1995, 1, 1)\n",
"end = sk.time(2022, 12, 31)\n",
"duration = end - start\n",
"timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])\n",
"\n",
"altitude = 400e3\n",
"rho_400, _temperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))\n",
"\n",
"altitude = 500e3\n",
"rho_500, emperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))\n",
"\n",
"\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"dates = [t.as_datetime() for t in timearray]\n",
"ax.semilogy(dates, rho_400, \"k-\", linewidth=1, label=\"Altitude = 400km\")\n",
"ax.semilogy(dates, rho_500, \"k--\", linewidth=1, label=\"Altitude = 500km\")\n",
"ax.set_xlabel(\"Year\")\n",
"ax.set_ylabel(\"Density [kg/m$^3$]\")\n",
"ax.set_title(\"Air Density Changes with Solar Cycle\")\n",
"ax.legend()\n",
"fig.autofmt_xdate()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forces Acting on a Satellite vs. Altitude\n",
"\n",
"This plot (inspired by Montenbruck & Gill, *Satellite Orbits*) shows the magnitude of various perturbation accelerations as a function of distance from the Earth's center. Key observations:\n",
"\n",
"- **Gravity** dominates at all altitudes\n",
"- **J2** (Earth's oblateness) is the largest perturbation, ~1000x stronger than higher-order terms\n",
"- **Atmospheric drag** is significant below ~800 km and varies with the solar cycle (max vs. min curves)\n",
"- **Solar radiation pressure** is roughly constant with altitude\n",
"- **Third-body** effects (Sun, Moon) grow with altitude as the satellite moves further from the Earth"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"## Force components acting on satellite vs altitude\n",
"import numpy as np\n",
"import satkit as sk\n",
"import math as m\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']\n",
"\n",
"\n",
"N = 1000\n",
"range_arr = np.logspace(m.log10(6378.2e3 + 100e3), m.log10(50e6), N)\n",
"grav1v = np.array([sk.gravity(np.array([a, 0, 0]), order=1) for a in range_arr])\n",
"grav1 = np.linalg.norm(grav1v, axis=1)\n",
"grav2v = np.array([sk.gravity(np.array([a, 0, 0]), order=2) for a in range_arr])\n",
"grav2 = np.linalg.norm(grav2v-grav1v, axis=1)\n",
"\n",
"grav6v = np.array([sk.gravity(np.array([a, 0, 0]), order=6) for a in range_arr])\n",
"grav5v = np.array([sk.gravity(np.array([a, 0, 0]), order=5) for a in range_arr])\n",
"grav6 = np.linalg.norm(grav6v-grav5v, axis=1)\n",
"\n",
"aoverm = 0.01\n",
"Cd = 2.2\n",
"Cr = 1.0\n",
"\n",
"# Drag force at maximum & minimum density\n",
"didx = np.argwhere(range_arr - sk.consts.earth_radius < 800e3).flatten()\n",
"tm_max = sk.time(2001, 12, 1)\n",
"tm_min = sk.time(1996, 12, 1)\n",
"rho_max = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_max)[0] for a in range_arr[didx]])\n",
"rho_min = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_min)[0] for a in range_arr[didx]])\n",
"varr = np.sqrt(sk.consts.mu_earth/(range_arr + sk.consts.earth_radius))\n",
"drag_max = 0.5 * rho_max * varr[didx]**2 * Cd * aoverm\n",
"drag_min = 0.5 * rho_min * varr[didx]**2 * Cd * aoverm\n",
"\n",
"moon_range = np.linalg.norm(sk.jplephem.geocentric_pos(sk.solarsystem.Moon, sk.time(2023, 1, 1)))\n",
"moon = sk.consts.mu_moon*((moon_range-range_arr)**(-2) - moon_range**(-2))\n",
"sun = sk.consts.mu_sun*((sk.consts.au-range_arr)**(-2) - sk.consts.au**(-2))\n",
"\n",
"a_radiation = 4.56e-6 * 0.5 * Cr * aoverm * np.ones(range_arr.shape)\n",
"\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 8))\n",
"def add_line(ax, x, y, text, frac=0.5, dx=-20, dy=-20):\n",
" ax.loglog(x, y, \"k-\", linewidth=1.5)\n",
" idx = int(len(x)*frac)\n",
" ax.annotate(text, xy=(x[idx], y[idx]), fontsize=9,\n",
" xytext=(dx, dy), textcoords=\"offset points\",\n",
" arrowprops=dict(arrowstyle=\"->\", color=\"gray\"))\n",
"\n",
"add_line(ax, range_arr/1e3, grav1/1e3, \"Gravity\")\n",
"add_line(ax, range_arr/1e3, grav2/1e3, \"J2\", 0.2, 0, -15)\n",
"add_line(ax, range_arr/1e3, grav6/1e3, \"J6\", 0.8, 0, -15)\n",
"add_line(ax, range_arr[didx]/1e3, drag_max/1e3, \"Drag Max\", 0.7, 30, 0)\n",
"add_line(ax, range_arr[didx]/1e3, drag_min/1e3, \"Drag Min\", 0.8, 10, 30)\n",
"add_line(ax, range_arr/1e3, moon/1e3, \"Moon\", 0.8, -10, -15)\n",
"add_line(ax, range_arr/1e3, sun/1e3, \"Sun\", 0.7, -10, 15)\n",
"add_line(ax, range_arr/1e3, a_radiation/1e3, \"Radiation\\nPressure\", 0.3, -10, 15)\n",
"ax.set_xlabel(\"Distance from Earth Origin [km]\")\n",
"ax.set_ylabel(\"Acceleration [km/s$^2$]\")\n",
"ax.set_title(\"Satellite Forces vs Altitude\")\n",
"ax.set_xlim(6378.1, 50e3)\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": "## Difference in Angle Between IERS 2010 Reduction and Approximate Calculation\n\n`satkit` provides both the full IERS 2010 precession/nutation model (`qgcrf2itrf`) and a faster approximate version (`qgcrf2itrf_approx`). This plot shows the angular difference between the two over several decades. The error remains below a few arcseconds, making the approximate model suitable for many applications where speed matters more than sub-arcsecond accuracy."
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"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']\n",
"import numpy as np\n",
"import satkit as sk\n",
"import math as m\n",
"\n",
"start = sk.time(2023, 1, 1)\n",
"end = sk.time(1970, 1, 1)\n",
"duration = end - start\n",
"timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])\n",
"qexact = sk.frametransform.qgcrf2itrf(timearray)\n",
"qapprox = sk.frametransform.qgcrf2itrf_approx(timearray)\n",
"qdiff = np.array([q1 * q2.conj for q1, q2 in zip(qexact, qapprox)]) # type: ignore\n",
"theta = np.array([min(q.angle, 2*m.pi - q.angle) for q in qdiff])\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"ax.plot([t.as_datetime() for t in timearray], theta*180.0/m.pi*3600, \"k-\", linewidth=1)\n",
"ax.set_xlabel(\"Year\")\n",
"ax.set_ylabel(\"Error [arcsec]\")\n",
"ax.set_title(\"Angle Error of Approximate ITRF to GCRF Rotation\")\n",
"fig.autofmt_xdate()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Polar Motion\n",
"\n",
"The Earth's rotational axis wobbles slightly relative to its crust, a phenomenon known as polar motion. The plot below shows the x and y components of polar motion (in arcseconds) over several decades. The dominant signal is the ~14-month Chandler wobble superimposed on an annual cycle, with a slow secular drift."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"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']\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"import matplotlib.dates as mdates\n",
"import numpy as np\n",
"import satkit as sk\n",
"import math as m\n",
"\n",
"start = sk.time(2023, 1, 1)\n",
"end = sk.time(1970, 1, 1)\n",
"duration = end - start\n",
"timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])\n",
"\n",
"dut1, xp, yp, lod, dX, dY = zip(*[sk.frametransform.earth_orientation_params(t) for t in timearray])\n",
"fig = plt.figure(figsize=(8, 7))\n",
"ax = fig.add_subplot(111, projection=\"3d\")\n",
"dates_num = mdates.date2num([t.as_datetime() for t in timearray])\n",
"ax.plot(list(xp), list(yp), dates_num, linewidth=0.5)\n",
"ax.set_xlabel(\"X [arcsec]\")\n",
"ax.set_ylabel(\"Y [arcsec]\")\n",
"ax.set_zlabel(\"Year\")\n",
"ax.set_title(\"Polar Motion\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Precession / Nutation\n",
"\n",
"Precession and nutation describe the long-term and short-term wobble of the Earth's rotational axis in inertial space. This is captured by the rotation between the GCRS (Geocentric Celestial Reference System) and CIRS (Celestial Intermediate Reference System) frames. The plot shows the nutation components (with the linear precession trend removed from x) over a 20-year period; the dominant ~18.6-year and semi-annual periods are clearly visible."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": "import satkit as sk\nimport numpy as np\nimport math as m\nimport matplotlib.pyplot as plt\nimport scienceplots # noqa: F401\nplt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n%config InlineBackend.figure_formats = ['svg']\nfrom mpl_toolkits.mplot3d import Axes3D\nimport matplotlib.dates as mdates\n\n# Precession / Nutation is the difference between the IERS 2010 GCRS frame and the CIRS frame\nstart = sk.time(2000, 1, 1)\nend = sk.time(2020, 1, 1)\nduration = end - start\nN = 1000\ntimearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, N)])\nqarray = np.array([sk.frametransform.qcirs2gcrf(t) for t in timearray])\ntheta = np.fromiter((q.angle for q in qarray), float)\nrots = np.array([q * np.array([0.0, 0.0, 1.0]) for q in qarray])\npline = np.linspace(0,1,N)\npv = np.polyfit(pline, rots[:,0], 1)\nrx0 =rots[:,0] - np.polyval(pv, pline)\n\nrad2asec = 180.0/m.pi * 3600\nfig = plt.figure(figsize=(8, 7))\nax = fig.add_subplot(111, projection=\"3d\")\ndates_num = mdates.date2num([t.as_datetime() for t in timearray])\nax.plot(rx0*rad2asec, rots[:,1]*rad2asec, dates_num, \"k-\", linewidth=0.5)\nax.set_xlabel(\"X [arcsec]\")\nax.set_ylabel(\"Y [arcsec]\")\nax.set_zlabel(\"Year\")\nax.set_title(\"Precession / Nutation\")\nplt.tight_layout()\nplt.show()"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}