SGP4 vs Numerical Propagation¶

SGP4 and numerical propagation are two fundamentally different approaches to predicting satellite orbits:

• SGP4 is an analytical model that uses simplified perturbation theory (J2--J4 zonal harmonics + atmospheric drag via B). It is fast and standardized, but limited in accuracy — TLEs are designed to be used *only with SGP4, and their accuracy degrades over days.
• Numerical propagation integrates the full equations of motion with configurable force models (high-order gravity, lunisolar perturbations, drag with NRLMSISE-00, solar radiation pressure). It is slower but can achieve meter-level accuracy.

This tutorial propagates the same orbit both ways and shows how and why they diverge.

Setup: Common Initial Conditions¶

We start from an ISS TLE. SGP4 propagates the TLE directly, outputting position in the TEME frame. For the numerical propagator, we convert the SGP4 state at TLE epoch from TEME to GCRF.

In [1]:

import satkit as sk
import numpy as np

# ISS TLE
tle = sk.TLE.from_lines([
    "ISS (ZARYA)",
    "1 25544U 98067A   24100.50000000  .00016717  00000-0  30233-3 0  9993",
    "2 25544  51.6416 247.4627 0006703 130.5360 229.6208 15.49438695449041"
])

epoch = tle.epoch
print(f"TLE epoch: {epoch}")
print(f"Mean motion: {tle.mean_motion:.8f} rev/day")
print(f"Inclination: {tle.inclination:.4f} deg")
print(f"B*: {tle.bstar:.6e}")

# SGP4 state at epoch (TEME frame)
pos_teme, vel_teme = sk.sgp4(tle, epoch)
print("\nSGP4 at epoch (TEME):")
print(f"  pos = [{pos_teme[0]/1e3:.3f}, {pos_teme[1]/1e3:.3f}, {pos_teme[2]/1e3:.3f}] km")
print(f"  vel = [{vel_teme[0]/1e3:.6f}, {vel_teme[1]/1e3:.6f}, {vel_teme[2]/1e3:.6f}] km/s")
print(f"  altitude ≈ {(np.linalg.norm(pos_teme) - sk.consts.earth_radius) / 1e3:.1f} km")

# Convert to GCRF for numerical propagator
q_teme2gcrf = sk.frametransform.qteme2gcrf(epoch)
pos_gcrf = q_teme2gcrf * pos_teme
vel_gcrf = q_teme2gcrf * vel_teme
state0 = np.concatenate([pos_gcrf, vel_gcrf])
print("\nInitial state (GCRF):")
print(f"  pos = [{pos_gcrf[0]/1e3:.3f}, {pos_gcrf[1]/1e3:.3f}, {pos_gcrf[2]/1e3:.3f}] km")

import satkit as sk import numpy as np # ISS TLE tle = sk.TLE.from_lines([ "ISS (ZARYA)", "1 25544U 98067A 24100.50000000 .00016717 00000-0 30233-3 0 9993", "2 25544 51.6416 247.4627 0006703 130.5360 229.6208 15.49438695449041" ]) epoch = tle.epoch print(f"TLE epoch: {epoch}") print(f"Mean motion: {tle.mean_motion:.8f} rev/day") print(f"Inclination: {tle.inclination:.4f} deg") print(f"B*: {tle.bstar:.6e}") # SGP4 state at epoch (TEME frame) pos_teme, vel_teme = sk.sgp4(tle, epoch) print("\nSGP4 at epoch (TEME):") print(f" pos = [{pos_teme[0]/1e3:.3f}, {pos_teme[1]/1e3:.3f}, {pos_teme[2]/1e3:.3f}] km") print(f" vel = [{vel_teme[0]/1e3:.6f}, {vel_teme[1]/1e3:.6f}, {vel_teme[2]/1e3:.6f}] km/s") print(f" altitude ≈ {(np.linalg.norm(pos_teme) - sk.consts.earth_radius) / 1e3:.1f} km") # Convert to GCRF for numerical propagator q_teme2gcrf = sk.frametransform.qteme2gcrf(epoch) pos_gcrf = q_teme2gcrf * pos_teme vel_gcrf = q_teme2gcrf * vel_teme state0 = np.concatenate([pos_gcrf, vel_gcrf]) print("\nInitial state (GCRF):") print(f" pos = [{pos_gcrf[0]/1e3:.3f}, {pos_gcrf[1]/1e3:.3f}, {pos_gcrf[2]/1e3:.3f}] km")

TLE epoch: 2024-04-09T12:00:00.000000Z
Mean motion: 15.49438695 rev/day
Inclination: 51.6416 deg
B*: 3.023300e-04

SGP4 at epoch (TEME):
  pos = [-2606.449, -6280.854, -0.050] km
  vel = [4.391336, -1.811001, 6.006375] km/s
  altitude ≈ 422.1 km

Initial state (GCRF):
  pos = [-2640.489, -6266.616, 6.392] km

Propagate Both for 3 Days¶

We propagate both models forward and compare positions in the GCRF frame.

In [2]:

import satkit as sk
import numpy as np

# Time array: 3 days, sampled every 2 minutes
n_days = 3
dt_minutes = 2
n_points = int(n_days * 24 * 60 / dt_minutes) + 1
times = [epoch + sk.duration(minutes=i * dt_minutes) for i in range(n_points)]

# --- SGP4 propagation ---
sgp4_teme_pos, sgp4_teme_vel = sk.sgp4(tle, times)

# Convert all SGP4 results from TEME to GCRF
q_arr = [sk.frametransform.qteme2gcrf(t) for t in times]
sgp4_gcrf = np.array([q * p for q, p in zip(q_arr, sgp4_teme_pos)])

# --- Numerical propagation ---
# High-fidelity force model
settings = sk.propsettings(
    gravity_model=sk.gravmodel.jgm3,
    gravity_degree=20,
    abs_error=1e-10,
    rel_error=1e-10,
)

result = sk.propagate(
    state0,
    epoch,
    end=epoch + sk.duration(days=n_days),
    propsettings=settings,
)

# Interpolate numerical solution at the same times
num_gcrf = np.array([result.interp(t)[:3] for t in times])

# --- Position difference ---
pos_diff = np.linalg.norm(sgp4_gcrf - num_gcrf, axis=1)
hours = np.array([(t - epoch).hours for t in times])

print("Position difference:")
print(f"  After  1 hour:  {pos_diff[30]:.1f} m")
print(f"  After  6 hours: {pos_diff[180]:.1f} m")
print(f"  After 24 hours: {pos_diff[720]:.0f} m")
print(f"  After 72 hours: {pos_diff[-1]:.0f} m")

import satkit as sk import numpy as np # Time array: 3 days, sampled every 2 minutes n_days = 3 dt_minutes = 2 n_points = int(n_days * 24 * 60 / dt_minutes) + 1 times = [epoch + sk.duration(minutes=i * dt_minutes) for i in range(n_points)] # --- SGP4 propagation --- sgp4_teme_pos, sgp4_teme_vel = sk.sgp4(tle, times) # Convert all SGP4 results from TEME to GCRF q_arr = [sk.frametransform.qteme2gcrf(t) for t in times] sgp4_gcrf = np.array([q * p for q, p in zip(q_arr, sgp4_teme_pos)]) # --- Numerical propagation --- # High-fidelity force model settings = sk.propsettings( gravity_model=sk.gravmodel.jgm3, gravity_degree=20, abs_error=1e-10, rel_error=1e-10, ) result = sk.propagate( state0, epoch, end=epoch + sk.duration(days=n_days), propsettings=settings, ) # Interpolate numerical solution at the same times num_gcrf = np.array([result.interp(t)[:3] for t in times]) # --- Position difference --- pos_diff = np.linalg.norm(sgp4_gcrf - num_gcrf, axis=1) hours = np.array([(t - epoch).hours for t in times]) print("Position difference:") print(f" After 1 hour: {pos_diff[30]:.1f} m") print(f" After 6 hours: {pos_diff[180]:.1f} m") print(f" After 24 hours: {pos_diff[720]:.0f} m") print(f" After 72 hours: {pos_diff[-1]:.0f} m")

Position difference:
  After  1 hour:  270.8 m
  After  6 hours: 817.0 m
  After 24 hours: 10426 m
  After 72 hours: 73699 m

Visualize the Divergence¶

In [3]:

import matplotlib.pyplot as plt
import scienceplots  # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']

fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

ax = axes[0]
ax.plot(hours, pos_diff / 1e3, linewidth=1.5)
ax.set_ylabel("Position Difference [km]")
ax.set_title("SGP4 vs Numerical Propagation (ISS)")

# Component-wise differences
ax = axes[1]
diff = sgp4_gcrf - num_gcrf
ax.plot(hours, diff[:, 0] / 1e3, label="X", linewidth=0.8)
ax.plot(hours, diff[:, 1] / 1e3, label="Y", linewidth=0.8)
ax.plot(hours, diff[:, 2] / 1e3, label="Z", linewidth=0.8)
ax.set_xlabel("Time Since Epoch [hours]")
ax.set_ylabel("Component Difference [km]")
ax.set_title("GCRF Component Differences")
ax.legend()

plt.tight_layout()
plt.show()

import matplotlib.pyplot as plt import scienceplots # noqa: F401 plt.style.use(["science", "no-latex", "../satkit.mplstyle"]) %config InlineBackend.figure_formats = ['svg'] fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True) ax = axes[0] ax.plot(hours, pos_diff / 1e3, linewidth=1.5) ax.set_ylabel("Position Difference [km]") ax.set_title("SGP4 vs Numerical Propagation (ISS)") # Component-wise differences ax = axes[1] diff = sgp4_gcrf - num_gcrf ax.plot(hours, diff[:, 0] / 1e3, label="X", linewidth=0.8) ax.plot(hours, diff[:, 1] / 1e3, label="Y", linewidth=0.8) ax.plot(hours, diff[:, 2] / 1e3, label="Z", linewidth=0.8) ax.set_xlabel("Time Since Epoch [hours]") ax.set_ylabel("Component Difference [km]") ax.set_title("GCRF Component Differences") ax.legend() plt.tight_layout() plt.show()

Effect of Gravity Model Fidelity¶

SGP4 only models J2--J4 (zonal harmonics up to degree 4). The numerical propagator can use much higher degree/order. Let's compare the divergence with different gravity truncation levels to see how much of the error comes from higher-order gravity terms.

In [4]:

import satkit as sk
import numpy as np

gravity_degrees = [4, 10, 20, 40]
results_by_degree = {}

for deg in gravity_degrees:
    settings = sk.propsettings(
        gravity_model=sk.gravmodel.jgm3,
        gravity_degree=deg,
        abs_error=1e-10,
        rel_error=1e-10,
    )
    res = sk.propagate(state0, epoch,
                       end=epoch + sk.duration(days=n_days),
                       propsettings=settings)
    num_pos = np.array([res.interp(t)[:3] for t in times])
    results_by_degree[deg] = np.linalg.norm(sgp4_gcrf - num_pos, axis=1)

fig, ax = plt.subplots(figsize=(10, 5))
for deg in gravity_degrees:
    ax.plot(hours, results_by_degree[deg] / 1e3, linewidth=1.2, label=f"Degree {deg}")
ax.set_xlabel("Time Since Epoch [hours]")
ax.set_ylabel("Position Difference [km]")
ax.set_title("SGP4 vs Numerical: Effect of Gravity Model Order")
ax.legend()
plt.tight_layout()
plt.show()

import satkit as sk import numpy as np gravity_degrees = [4, 10, 20, 40] results_by_degree = {} for deg in gravity_degrees: settings = sk.propsettings( gravity_model=sk.gravmodel.jgm3, gravity_degree=deg, abs_error=1e-10, rel_error=1e-10, ) res = sk.propagate(state0, epoch, end=epoch + sk.duration(days=n_days), propsettings=settings) num_pos = np.array([res.interp(t)[:3] for t in times]) results_by_degree[deg] = np.linalg.norm(sgp4_gcrf - num_pos, axis=1) fig, ax = plt.subplots(figsize=(10, 5)) for deg in gravity_degrees: ax.plot(hours, results_by_degree[deg] / 1e3, linewidth=1.2, label=f"Degree {deg}") ax.set_xlabel("Time Since Epoch [hours]") ax.set_ylabel("Position Difference [km]") ax.set_title("SGP4 vs Numerical: Effect of Gravity Model Order") ax.legend() plt.tight_layout() plt.show()

Altitude Comparison¶

Another way to see the divergence: compare the altitude profiles. SGP4's simplified drag model (B*) produces slightly different orbital decay than the NRLMSISE-00 density model used by the numerical propagator.

In [5]:

import satkit as sk
import numpy as np

# SGP4 altitudes (via ITRF -> geodetic)
sgp4_alt = []
for t, p_teme in zip(times, sgp4_teme_pos):
    p_itrf = sk.frametransform.qteme2itrf(t) * p_teme
    sgp4_alt.append(sk.itrfcoord(p_itrf).altitude / 1e3)

# Numerical altitudes (GCRF -> ITRF -> geodetic)
num_alt = []
for t in times:
    state = result.interp(t)
    p_itrf = sk.frametransform.qgcrf2itrf(t) * state[:3]
    num_alt.append(sk.itrfcoord(p_itrf).altitude / 1e3)

fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)

ax = axes[0]
ax.plot(hours, sgp4_alt, linewidth=0.5, alpha=0.7, label="SGP4")
ax.plot(hours, num_alt, linewidth=0.5, alpha=0.7, label="Numerical")
ax.set_ylabel("Altitude [km]")
ax.set_title("Altitude Comparison")
ax.legend()

ax = axes[1]
alt_diff = np.array(sgp4_alt) - np.array(num_alt)
ax.plot(hours, alt_diff * 1e3, linewidth=0.8)
ax.set_xlabel("Time Since Epoch [hours]")
ax.set_ylabel("Altitude Difference [m]")
ax.set_title("SGP4 $-$ Numerical Altitude Difference")

plt.tight_layout()
plt.show()

import satkit as sk import numpy as np # SGP4 altitudes (via ITRF -> geodetic) sgp4_alt = [] for t, p_teme in zip(times, sgp4_teme_pos): p_itrf = sk.frametransform.qteme2itrf(t) * p_teme sgp4_alt.append(sk.itrfcoord(p_itrf).altitude / 1e3) # Numerical altitudes (GCRF -> ITRF -> geodetic) num_alt = [] for t in times: state = result.interp(t) p_itrf = sk.frametransform.qgcrf2itrf(t) * state[:3] num_alt.append(sk.itrfcoord(p_itrf).altitude / 1e3) fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True) ax = axes[0] ax.plot(hours, sgp4_alt, linewidth=0.5, alpha=0.7, label="SGP4") ax.plot(hours, num_alt, linewidth=0.5, alpha=0.7, label="Numerical") ax.set_ylabel("Altitude [km]") ax.set_title("Altitude Comparison") ax.legend() ax = axes[1] alt_diff = np.array(sgp4_alt) - np.array(num_alt) ax.plot(hours, alt_diff * 1e3, linewidth=0.8) ax.set_xlabel("Time Since Epoch [hours]") ax.set_ylabel("Altitude Difference [m]") ax.set_title("SGP4 $-$ Numerical Altitude Difference") plt.tight_layout() plt.show()