Interesting Plots¶

A collection of visualizations demonstrating satkit's atmospheric density models, gravity field, and coordinate frame transformations.

Air Density as a Function of the Solar Cycle¶

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.

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).

In [1]:

import satkit as sk
import plotly.graph_objects as go
import numpy as np
import math as m

start = sk.time(1995, 1, 1)
end = sk.time(2022, 12, 31)
duration = end - start
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])

altitude = 400e3
rho_400, _temperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))

altitude = 500e3
rho_500, emperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))


fig = go.Figure()
fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_400,
                          mode='lines', name='Altitude = 400km', line=dict(color='black', width=1)))
fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_500,
                         mode='lines', name='Altitude = 500km', line=dict(color='black', width=1, dash='dash')))
fig.update_layout(legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.5
))
fig.update_yaxes(type="log", title_font=dict(size=8), tickfont=dict(size=8), title='Density [kg/m<sup>3</sup>]')
fig.update_xaxes(title='Year')
fig.update_layout(yaxis_tickformat=".2e", width=650, height=512,
                  title='Air Density Changes with Solar Cycle',)
fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_layout(xaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),
    yaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),)
fig.show()

import satkit as sk import plotly.graph_objects as go import numpy as np import math as m start = sk.time(1995, 1, 1) end = sk.time(2022, 12, 31) duration = end - start timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)]) altitude = 400e3 rho_400, _temperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray])) altitude = 500e3 rho_500, emperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray])) fig = go.Figure() fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_400, mode='lines', name='Altitude = 400km', line=dict(color='black', width=1))) fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_500, mode='lines', name='Altitude = 500km', line=dict(color='black', width=1, dash='dash'))) fig.update_layout(legend=dict( yanchor="top", y=0.99, xanchor="left", x=0.5 )) fig.update_yaxes(type="log", title_font=dict(size=8), tickfont=dict(size=8), title='Density [kg/m³]') fig.update_xaxes(title='Year') fig.update_layout(yaxis_tickformat=".2e", width=650, height=512, title='Air Density Changes with Solar Cycle',) fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True) fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True) fig.update_layout(xaxis=dict( gridcolor='#dddddd', gridwidth=1, ), yaxis=dict( gridcolor='#dddddd', gridwidth=1, ),) fig.show()

Forces Acting on a Satellite vs. Altitude¶

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:

• Gravity dominates at all altitudes
• J2 (Earth's oblateness) is the largest perturbation, ~1000x stronger than higher-order terms
• Atmospheric drag is significant below ~800 km and varies with the solar cycle (max vs. min curves)
• Solar radiation pressure is roughly constant with altitude
• Third-body effects (Sun, Moon) grow with altitude as the satellite moves further from the Earth

In [2]:

## Force components acting on satellite vs altitude
import numpy as np
import satkit as sk
import math as m
import plotly.graph_objects as go


N = 1000
range_arr = np.logspace(m.log10(6378.2e3 + 100e3), m.log10(50e6), N)
grav1v = np.array([sk.gravity(np.array([a, 0, 0]), order=1) for a in range_arr])
grav1 = np.linalg.norm(grav1v, axis=1)
grav2v = np.array([sk.gravity(np.array([a, 0, 0]), order=2) for a in range_arr])
grav2 = np.linalg.norm(grav2v-grav1v, axis=1)

grav6v = np.array([sk.gravity(np.array([a, 0, 0]), order=6) for a in range_arr])
grav5v = np.array([sk.gravity(np.array([a, 0, 0]), order=5) for a in range_arr])
grav6 = np.linalg.norm(grav6v-grav5v, axis=1)

aoverm = 0.01
Cd = 2.2
Cr = 1.0

# Drag force at maximum & minimum density
didx = np.argwhere(range_arr - sk.consts.earth_radius < 800e3).flatten()
tm_max = sk.time(2001, 12, 1)
tm_min = sk.time(1996, 12, 1)
rho_max = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_max)[0] for a in range_arr[didx]])
rho_min = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_min)[0] for a in range_arr[didx]])
varr = np.sqrt(sk.consts.mu_earth/(range_arr + sk.consts.earth_radius))
drag_max = 0.5 * rho_max * varr[didx]**2 * Cd * aoverm
drag_min = 0.5 * rho_min * varr[didx]**2 * Cd * aoverm

moon_range = np.linalg.norm(sk.jplephem.geocentric_pos(sk.solarsystem.Moon, sk.time(2023, 1, 1)))
moon = sk.consts.mu_moon*((moon_range-range_arr)**(-2) - moon_range**(-2))
sun = sk.consts.mu_sun*((sk.consts.au-range_arr)**(-2) - sk.consts.au**(-2))

a_radiation = 4.56e-6 * 0.5 * Cr * aoverm * np.ones(range_arr.shape)

def add_line(fig, x, y, text, frac=0.5, ax=-20, ay=-30):
    fig.add_scatter(x=x, y=y, mode='lines', name=text, line=dict(width=2, color='black'))
    idx = int(len(x)*frac)
    fig.add_annotation(
        x=np.log10(x[idx]),
        y=np.log10(y[idx]),
        text=text,
        showarrow=True,
        font=dict(
            size=12,
            color='black',
        ),
        align="center",
        ax=ax,
        ay=ay,
        arrowcolor="#636363",
        arrowsize=1,
        arrowwidth=2,
        arrowhead=6,
        borderwidth=2,
        borderpad=4,
        opacity=0.8
    )

fig = go.Figure()
add_line(fig, range_arr/1e3, grav1/1e3, 'Gravity')
add_line(fig, range_arr/1e3, grav2/1e3, 'J2', 0.2, 0, -10)
add_line(fig, range_arr/1e3, grav6/1e3, 'J6', 0.8, 0, -10)
add_line(fig, range_arr[didx]/1e3, drag_max/1e3, 'Drag Max', 0.7, 30, 0)
add_line(fig, range_arr[didx]/1e3, drag_min/1e3, 'Drag Min', 0.8, 10, 30)
add_line(fig, range_arr/1e3, moon/1e3, 'Moon', 0.8, -10, -10)
add_line(fig, range_arr/1e3, sun/1e3, 'Sun', 0.7, -10, 10)
add_line(fig, range_arr/1e3, a_radiation/1e3, 'Radiation\nPressure', 0.3, -10, 10)

fig.update_xaxes(type="log", title='Distance from Earth Origin [km]', range=np.log10([6378.1, 50e3]))
fig.update_yaxes(type="log", title='Acceleration   [km/s<sup>2</sup>]', tickformat=".1e")
fig.update_layout(title='Satellite Forces vs Altitude', width=650, height=650)
fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_layout(showlegend=False, xaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),
    yaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),)
fig.show()

## Force components acting on satellite vs altitude import numpy as np import satkit as sk import math as m import plotly.graph_objects as go N = 1000 range_arr = np.logspace(m.log10(6378.2e3 + 100e3), m.log10(50e6), N) grav1v = np.array([sk.gravity(np.array([a, 0, 0]), order=1) for a in range_arr]) grav1 = np.linalg.norm(grav1v, axis=1) grav2v = np.array([sk.gravity(np.array([a, 0, 0]), order=2) for a in range_arr]) grav2 = np.linalg.norm(grav2v-grav1v, axis=1) grav6v = np.array([sk.gravity(np.array([a, 0, 0]), order=6) for a in range_arr]) grav5v = np.array([sk.gravity(np.array([a, 0, 0]), order=5) for a in range_arr]) grav6 = np.linalg.norm(grav6v-grav5v, axis=1) aoverm = 0.01 Cd = 2.2 Cr = 1.0 # Drag force at maximum & minimum density didx = np.argwhere(range_arr - sk.consts.earth_radius < 800e3).flatten() tm_max = sk.time(2001, 12, 1) tm_min = sk.time(1996, 12, 1) rho_max = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_max)[0] for a in range_arr[didx]]) rho_min = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_min)[0] for a in range_arr[didx]]) varr = np.sqrt(sk.consts.mu_earth/(range_arr + sk.consts.earth_radius)) drag_max = 0.5 * rho_max * varr[didx]**2 * Cd * aoverm drag_min = 0.5 * rho_min * varr[didx]**2 * Cd * aoverm moon_range = np.linalg.norm(sk.jplephem.geocentric_pos(sk.solarsystem.Moon, sk.time(2023, 1, 1))) moon = sk.consts.mu_moon*((moon_range-range_arr)**(-2) - moon_range**(-2)) sun = sk.consts.mu_sun*((sk.consts.au-range_arr)**(-2) - sk.consts.au**(-2)) a_radiation = 4.56e-6 * 0.5 * Cr * aoverm * np.ones(range_arr.shape) def add_line(fig, x, y, text, frac=0.5, ax=-20, ay=-30): fig.add_scatter(x=x, y=y, mode='lines', name=text, line=dict(width=2, color='black')) idx = int(len(x)*frac) fig.add_annotation( x=np.log10(x[idx]), y=np.log10(y[idx]), text=text, showarrow=True, font=dict( size=12, color='black', ), align="center", ax=ax, ay=ay, arrowcolor="#636363", arrowsize=1, arrowwidth=2, arrowhead=6, borderwidth=2, borderpad=4, opacity=0.8 ) fig = go.Figure() add_line(fig, range_arr/1e3, grav1/1e3, 'Gravity') add_line(fig, range_arr/1e3, grav2/1e3, 'J2', 0.2, 0, -10) add_line(fig, range_arr/1e3, grav6/1e3, 'J6', 0.8, 0, -10) add_line(fig, range_arr[didx]/1e3, drag_max/1e3, 'Drag Max', 0.7, 30, 0) add_line(fig, range_arr[didx]/1e3, drag_min/1e3, 'Drag Min', 0.8, 10, 30) add_line(fig, range_arr/1e3, moon/1e3, 'Moon', 0.8, -10, -10) add_line(fig, range_arr/1e3, sun/1e3, 'Sun', 0.7, -10, 10) add_line(fig, range_arr/1e3, a_radiation/1e3, 'Radiation\nPressure', 0.3, -10, 10) fig.update_xaxes(type="log", title='Distance from Earth Origin [km]', range=np.log10([6378.1, 50e3])) fig.update_yaxes(type="log", title='Acceleration [km/s²]', tickformat=".1e") fig.update_layout(title='Satellite Forces vs Altitude', width=650, height=650) fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True) fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True) fig.update_layout(showlegend=False, xaxis=dict( gridcolor='#dddddd', gridwidth=1, ), yaxis=dict( gridcolor='#dddddd', gridwidth=1, ),) fig.show()

Difference in Angle Between IAU-2010 Reduction and Approximate Calculation¶

satkit provides both the full IAU-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.

In [3]:

import plotly.graph_objects as go
import numpy as np
import satkit as sk
import math as m

start = sk.time(2023, 1, 1)
end = sk.time(1970, 1, 1)
duration = end - start
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])
qexact = sk.frametransform.qgcrf2itrf(timearray)
qapprox = sk.frametransform.qgcrf2itrf_approx(timearray)
qdiff = np.array([q1 * q2.conj for q1, q2 in zip(qexact, qapprox)]) # type: ignore
theta = np.array([q.angle for q in qdiff])
fig = go.Figure()
fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=theta*180.0/m.pi*3600,
                         mode='lines', name='Angle Difference',
                         line=dict(width=1, color='black')))
fig.update_layout(title='Angle Error of Approximate ITRF to GCRF Rotation', width=650, height=512,
                  font=dict(size=14))
fig.update_yaxes(title='Error [arcsec]')
fig.update_xaxes(title='Year')

import plotly.graph_objects as go import numpy as np import satkit as sk import math as m start = sk.time(2023, 1, 1) end = sk.time(1970, 1, 1) duration = end - start timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)]) qexact = sk.frametransform.qgcrf2itrf(timearray) qapprox = sk.frametransform.qgcrf2itrf_approx(timearray) qdiff = np.array([q1 * q2.conj for q1, q2 in zip(qexact, qapprox)]) # type: ignore theta = np.array([q.angle for q in qdiff]) fig = go.Figure() fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=theta*180.0/m.pi*3600, mode='lines', name='Angle Difference', line=dict(width=1, color='black'))) fig.update_layout(title='Angle Error of Approximate ITRF to GCRF Rotation', width=650, height=512, font=dict(size=14)) fig.update_yaxes(title='Error [arcsec]') fig.update_xaxes(title='Year')

Polar Motion¶

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.

In [4]:

import plotly.graph_objects as go
import numpy as np
import satkit as sk
import math as m

start = sk.time(2023, 1, 1)
end = sk.time(1970, 1, 1)
duration = end - start
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])

dut1, xp, yp, lod, dX, dY = zip(*[sk.frametransform.earth_orientation_params(t) for t in timearray])
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=xp, y=yp, z=[t.datetime() for t in timearray],  mode='lines', name='Polar Motion',))
fig.update_scenes(xaxis_title='X [arcsec]', yaxis_title='Y [arcsec]', zaxis_title='Year')
fig.update_layout(title='Polar Motion', width=650, height=600, font=dict(size=12))
fig.show()

import plotly.graph_objects as go import numpy as np import satkit as sk import math as m start = sk.time(2023, 1, 1) end = sk.time(1970, 1, 1) duration = end - start timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)]) dut1, xp, yp, lod, dX, dY = zip(*[sk.frametransform.earth_orientation_params(t) for t in timearray]) fig = go.Figure() fig.add_trace(go.Scatter3d(x=xp, y=yp, z=[t.datetime() for t in timearray], mode='lines', name='Polar Motion',)) fig.update_scenes(xaxis_title='X [arcsec]', yaxis_title='Y [arcsec]', zaxis_title='Year') fig.update_layout(title='Polar Motion', width=650, height=600, font=dict(size=12)) fig.show()

Precession / Nutation¶

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.

In [5]:

import satkit as sk
import numpy as np
import math as m
import plotly.graph_objects as go

# Precession / Nutation is the difference between the IAU-2006 GCRS frame and the CIRS frame
start = sk.time(2000, 1, 1)
end = sk.time(2020, 1, 1)
duration = end - start
N = 1000
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, N)])
qarray = np.array([sk.frametransform.qcirs2gcrf(t) for t in timearray])
theta = np.fromiter((q.angle for q in qarray), float)
rots = np.array([q * np.array([0.0, 0.0, 1.0]) for q in qarray])
pline = np.linspace(0,1,N)
pv = np.polyfit(pline, rots[:,0], 1)
rx0 =rots[:,0] - np.polyval(pv, pline)

rad2asec = 180.0/m.pi * 3600
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=rx0*rad2asec, y=rots[:,1]*rad2asec, z=[t.datetime() for t in timearray], mode='lines',
                          name='Precession/Nutation', line=dict(width=2, color='black')))
fig.update_layout(title='Precession / Nutation', width=650, height=600, font=dict(size=12))

fig.show()

import satkit as sk import numpy as np import math as m import plotly.graph_objects as go # Precession / Nutation is the difference between the IAU-2006 GCRS frame and the CIRS frame start = sk.time(2000, 1, 1) end = sk.time(2020, 1, 1) duration = end - start N = 1000 timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, N)]) qarray = np.array([sk.frametransform.qcirs2gcrf(t) for t in timearray]) theta = np.fromiter((q.angle for q in qarray), float) rots = np.array([q * np.array([0.0, 0.0, 1.0]) for q in qarray]) pline = np.linspace(0,1,N) pv = np.polyfit(pline, rots[:,0], 1) rx0 =rots[:,0] - np.polyval(pv, pline) rad2asec = 180.0/m.pi * 3600 fig = go.Figure() fig.add_trace(go.Scatter3d(x=rx0*rad2asec, y=rots[:,1]*rad2asec, z=[t.datetime() for t in timearray], mode='lines', name='Precession/Nutation', line=dict(width=2, color='black'))) fig.update_layout(title='Precession / Nutation', width=650, height=600, font=dict(size=12)) fig.show()