GPS Orbit Fit Example¶

This example walks through a full GPS-class orbit determination workflow: starting from precise SP3 ephemeris truth data, fit an initial state plus drag/SRP parameters and a small empirical acceleration, then validate the propagated orbit against the truth and compare ODE integrators on the same problem.

Overview¶

1. Load reference data — read GPS satellite positions from an SP3 file of precise orbit determination results (positions only, no velocity).
2. Rotate to inertial frame — transform the ITRF positions into GCRF.
3. Fit initial conditions — use a linearized least-squares loop over the state transition matrix to recover the 6-DOF initial state, wrapped in a non-linear solve for the solar radiation pressure coefficient ($C_r A/m$) and a constant empirical acceleration in the NTW frame.
4. Validate — compare the propagated positions against SP3 truth and plot the residuals.
5. Fit vs. free-run — fit over day 1 only, then propagate through day 2 with no new measurements to see how errors grow outside the fit window.
6. Compare integrators — benchmark several Runge-Kutta integrators and the Gauss-Jackson 8 multistep predictor-corrector on the same problem.

The SP3 file contains a full 24 hours of satellite positions, recorded once every 5 minutes.

In [1]:

# Imports
import gzip
import os
import time
import urllib.request

import matplotlib.pyplot as plt
import numpy as np
import scienceplots  # noqa: F401
from scipy.optimize import minimize

import satkit as sk

plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']

# Imports import gzip import os import time import urllib.request import matplotlib.pyplot as plt import numpy as np import scienceplots # noqa: F401 from scipy.optimize import minimize import satkit as sk plt.style.use(["science", "no-latex", "../satkit.mplstyle"]) %config InlineBackend.figure_formats = ['svg']

Setup and SP3 File Reader¶

The SP3 file format (Hilla 2016, SP3-d) contains precise satellite positions (and optionally velocities) at regular intervals, typically generated by post-processing networks of ground stations. We read the file and extract ITRF positions for a single GPS satellite.

In [2]:

# Function to read in the SP3 file


def read_sp3file(fname, satnum=20):
    """
    Read SP3 file
    (file containing "true" GPS ephemerides)
    and output UTC time and position in ITRF frame
    """

    # Read in the test vectors
    with open(fname, "r") as fd:
        lines = fd.readlines()

    def line2date(lines):
        for line in lines:
            year = int(line[3:7])
            month = int(line[8:10])
            day = int(line[11:13])
            hour = int(line[14:16])
            minute = int(line[17:19])
            sec = float(line[20:32])
            yield sk.time(year, month, day, hour, minute, sec)

    def line2pos(lines):
        for line in lines:
            lvals = line.split()
            yield np.array([float(lvals[1]), float(lvals[2]), float(lvals[3])])

    datelines = list(filter(lambda x: x[0] == "*", lines))
    match = f"PG{satnum:02d}"
    satlines = list(filter(lambda x: x[0:4] == match, lines))
    dates = np.fromiter(line2date(datelines), sk.time)
    pitrf = np.stack(np.fromiter(line2pos(satlines), list), axis=0) * 1.0e3  # type: ignore

    return (pitrf, dates)

# Function to read in the SP3 file def read_sp3file(fname, satnum=20): """ Read SP3 file (file containing "true" GPS ephemerides) and output UTC time and position in ITRF frame """ # Read in the test vectors with open(fname, "r") as fd: lines = fd.readlines() def line2date(lines): for line in lines: year = int(line[3:7]) month = int(line[8:10]) day = int(line[11:13]) hour = int(line[14:16]) minute = int(line[17:19]) sec = float(line[20:32]) yield sk.time(year, month, day, hour, minute, sec) def line2pos(lines): for line in lines: lvals = line.split() yield np.array([float(lvals[1]), float(lvals[2]), float(lvals[3])]) datelines = list(filter(lambda x: x[0] == "*", lines)) match = f"PG{satnum:02d}" satlines = list(filter(lambda x: x[0:4] == match, lines)) dates = np.fromiter(line2date(datelines), sk.time) pitrf = np.stack(np.fromiter(line2pos(satlines), list), axis=0) * 1.0e3 # type: ignore return (pitrf, dates)

Load SP3 Data and Seed the Initial State¶

Download the day-1 SP3 file, rotate the ITRF positions into GCRF, and build a crude initial state from the first two points. The velocity seed is just a finite difference across 5 minutes — coarse, but good enough for the fit loop to converge.

In [3]:

# Download SP3 file if not present
fname = "./ESA0OPSFIN_20233640000_01D_05M_ORB.SP3"
url = "http://navigation-office.esa.int/products/gnss-products/2294/ESA0OPSFIN_20233640000_01D_05M_ORB.SP3.gz"

if not os.path.exists(fname):
    with urllib.request.urlopen(url) as response:
        with open(fname, "wb") as out_file:
            with gzip.GzipFile(fileobj=response) as uncompressed:
                out_file.write(uncompressed.read())

# Read in the SP3 file
[pitrf, timearr] = read_sp3file(fname)

# Rotate positions to the GCRF frame
pgcrf = np.stack(
    np.fromiter(
        (q * p for q, p in zip([sk.frametransform.rotation(sk.frame.ITRF, sk.frame.GCRF, t) for t in timearr], pitrf)),
        list,  # type: ignore
    ),
    axis=0,
)  # type: ignore

# Crude estimation of initial velocity from the first two position states.
# The 5-minute spacing makes this rough, but the fit loop will refine it.
vgcrf = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds

# Initial 6-DOF state (position + velocity) in GCRF
state0 = np.concatenate((pgcrf[0, :], vgcrf))

# Propagator settings shared across fits.
# Available integrators: sk.integrator.rkv98 (default), rkv98_nointerp, rkv87, rkv65, rkts54
# Available gravity models: sk.gravmodel.egm96 (default), jgm3, jgm2, itugrace16
settings = sk.propsettings(
    abs_error=1e-12,
    rel_error=1e-12,
    gravity_degree=10,
)
# Only compute sun, moon positions and earth rotation vectors once for all propagations
settings.precompute_terms(timearr[0], timearr[-1])

# Download SP3 file if not present fname = "./ESA0OPSFIN_20233640000_01D_05M_ORB.SP3" url = "http://navigation-office.esa.int/products/gnss-products/2294/ESA0OPSFIN_20233640000_01D_05M_ORB.SP3.gz" if not os.path.exists(fname): with urllib.request.urlopen(url) as response: with open(fname, "wb") as out_file: with gzip.GzipFile(fileobj=response) as uncompressed: out_file.write(uncompressed.read()) # Read in the SP3 file [pitrf, timearr] = read_sp3file(fname) # Rotate positions to the GCRF frame pgcrf = np.stack( np.fromiter( (q * p for q, p in zip([sk.frametransform.rotation(sk.frame.ITRF, sk.frame.GCRF, t) for t in timearr], pitrf)), list, # type: ignore ), axis=0, ) # type: ignore # Crude estimation of initial velocity from the first two position states. # The 5-minute spacing makes this rough, but the fit loop will refine it. vgcrf = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds # Initial 6-DOF state (position + velocity) in GCRF state0 = np.concatenate((pgcrf[0, :], vgcrf)) # Propagator settings shared across fits. # Available integrators: sk.integrator.rkv98 (default), rkv98_nointerp, rkv87, rkv65, rkts54 # Available gravity models: sk.gravmodel.egm96 (default), jgm3, jgm2, itugrace16 settings = sk.propsettings( abs_error=1e-12, rel_error=1e-12, gravity_degree=10, ) # Only compute sun, moon positions and earth rotation vectors once for all propagations settings.precompute_terms(timearr[0], timearr[-1])

Fit Initial State and Validate¶

We fit the orbit in two nested loops:

• Inner loop (linearized least squares) — iteratively solves for the 6-element initial state correction $\Delta s_0$ using the position block of the state transition matrix $\Phi_p(t; s_0)$:

$$\min_{\Delta s_0} \; \sum_t \| p_\text{gps}(t) - \hat{p}(t; s_0) - \Phi_p(t; s_0)\,\Delta s_0 \|^2$$

where $p_\text{gps}(t)$ is the SP3 position in GCRF and $\hat{p}(t; s_0)$ is the propagated position from initial state $s_0$. Converges in about 5 iterations.

• Outer loop (Nelder-Mead) — non-linearly searches for parameters the STM doesn't cover: the solar radiation pressure coefficient $C_r A/m$ and a constant empirical acceleration in the velocity-aligned NTW frame $[a_N, a_T, a_W]$.

The empirical acceleration is the interesting physical piece. It is not a real thruster — it is a catch-all for any constant bias left over after the modelled forces are applied. satkit's default force model includes static spherical-harmonic Earth gravity, Sun and Moon third-body point-mass gravity, atmospheric drag (NRLMSISE-00), solar radiation pressure with eclipse modelling, solid Earth tides (IERS 2010 §6.2.1 Step 1), and the general-relativistic correction of IERS 2010 §10.3 Eq. 10.12 (Schwarzschild, geodesic precession and Lense–Thirring; Petit & Luzum 2010). The empirical absorbs what's not covered: ocean tides, higher-order SRP effects (box-wing geometry, antenna thrust, thermal re-radiation), Earth albedo radiation pressure, and any residual Earth-orientation errors. Because we use the NTW frame, the tangential component $a_T$ acts purely along the velocity vector and therefore maps one-to-one to orbital energy (and hence semi-major-axis drift) — it is the most effective single parameter for absorbing along-track residual growth, which is the dominant error mode in GNSS orbit determination. In operational orbit determination these are called "empirical accelerations" and fitting them is standard practice — e.g. the CODE 5-parameter SRP model (Beutler et al. 1994). The batch least-squares formulation is that of Tapley, Schutz & Born (2004), §4.3; the outer loop is the simplex method of Nelder & Mead (1965).

In [4]:

# Fit helpers: satproperties builder and linearized LSQ loop


def make_satprops(craoverm, ntw_accel):
    """Build a satproperties with a constant NTW-frame thrust acceleration.

    NTW axes (Vallado §3.3): T = velocity unit vector, W = orbit-normal
    (r × v), N = W × T (in-plane, normal to velocity).  Strict along-velocity
    semantics even for eccentric orbits — the T component acts purely on
    orbital energy (semi-major axis).
    """
    thrusts = [
        sk.thrust.constant(
            list(ntw_accel), timearr[0], timearr[-1], frame=sk.frame.NTW
        )
    ]
    return sk.satproperties(craoverm=craoverm, thrusts=thrusts)


def linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp, iters=5):
    """
    Linearized least squares fit of initial state to SP3 truth data, holding
    CrAoverM and any thrust terms fixed.
    """
    state0_s = state0.copy()
    for idx in range(iters):
        # Propagate state and state transition matrix over times of interest
        res = sk.propagate(
            state0_s,
            timearr[0],
            timearr[-1],
            output_phi=True,
            propsettings=settings,
            satproperties=sp,
        )

        # Get state and state transition matrix at times of GPS truth data
        statearr, phiarr = zip(*[res.interp(t, output_phi=True) for t in timearr])
        phiarr = np.array(phiarr)
        statearr = np.array(statearr)

        # Linearized least squares solve for state0 update
        H = np.sum([p[0:3, :].T @ p[0:3, :] for p in phiarr], axis=0)
        b = np.sum(
            [
                p[0:3, :].T @ (pgcrf[i, :] - statearr[i, 0:3]).T
                for i, p in enumerate(phiarr)
            ],
            axis=0,
        )
        dstate0 = np.linalg.solve(H, b)
        state0_s = state0_s + dstate0

    perr = np.zeros((len(timearr), 3))
    for i in range(len(timearr)):
        perr[i, :] = res.interp(timearr[i])[0:3] - pgcrf[i, :]

    return state0_s, res, perr

# Fit helpers: satproperties builder and linearized LSQ loop def make_satprops(craoverm, ntw_accel): """Build a satproperties with a constant NTW-frame thrust acceleration. NTW axes (Vallado §3.3): T = velocity unit vector, W = orbit-normal (r × v), N = W × T (in-plane, normal to velocity). Strict along-velocity semantics even for eccentric orbits — the T component acts purely on orbital energy (semi-major axis). """ thrusts = [ sk.thrust.constant( list(ntw_accel), timearr[0], timearr[-1], frame=sk.frame.NTW ) ] return sk.satproperties(craoverm=craoverm, thrusts=thrusts) def linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp, iters=5): """ Linearized least squares fit of initial state to SP3 truth data, holding CrAoverM and any thrust terms fixed. """ state0_s = state0.copy() for idx in range(iters): # Propagate state and state transition matrix over times of interest res = sk.propagate( state0_s, timearr[0], timearr[-1], output_phi=True, propsettings=settings, satproperties=sp, ) # Get state and state transition matrix at times of GPS truth data statearr, phiarr = zip(*[res.interp(t, output_phi=True) for t in timearr]) phiarr = np.array(phiarr) statearr = np.array(statearr) # Linearized least squares solve for state0 update H = np.sum([p[0:3, :].T @ p[0:3, :] for p in phiarr], axis=0) b = np.sum( [ p[0:3, :].T @ (pgcrf[i, :] - statearr[i, 0:3]).T for i, p in enumerate(phiarr) ], axis=0, ) dstate0 = np.linalg.solve(H, b) state0_s = state0_s + dstate0 perr = np.zeros((len(timearr), 3)) for i in range(len(timearr)): perr[i, :] = res.interp(timearr[i])[0:3] - pgcrf[i, :] return state0_s, res, perr

Outer loop: fit $C_r A/m$ and empirical NTW acceleration¶

Nelder-Mead searches over four scalars: $C_r A/m$ and $[a_N, a_T, a_W]$. At each outer step the inner linearized LSQ re-fits the 6-element initial state (the STM only covers position/velocity, so these acceleration parameters must be handled by the outer non-linear loop). The NTW accel components are scaled by ACCEL_SCALE so all parameters are $O(0.01)$ and Nelder-Mead behaves.

In [5]:

# Nelder-Mead outer optimization over CrA/m and constant NTW empirical acceleration
ACCEL_SCALE = 1e-9  # m/s^2 per unit of scaled parameter


def minfunc_full(params, state0, timearr, pgcrf, settings):
    craoverm = params[0]
    ntw_accel = np.array(params[1:4]) * ACCEL_SCALE
    if craoverm < 0 or craoverm > 1:
        return 1e30
    sp = make_satprops(craoverm, ntw_accel)
    _, _, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp)
    return float(np.sum(perr**2))


x0 = np.array([0.01, 0.0, 0.0, 0.0])
r = minimize(
    minfunc_full,
    x0,
    args=(state0, timearr, pgcrf, settings),
    method="Nelder-Mead",
    options={"xatol": 1e-5, "fatol": 1e-4, "maxiter": 400},
)

craoverm_fit = r.x[0]
ntw_accel_fit = np.array(r.x[1:4]) * ACCEL_SCALE

# Final least-squares fit with optimized CrA/m and NTW empirical acceleration
sp = make_satprops(craoverm_fit, ntw_accel_fit)
state0, res, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp)

print(f"Satellite radiation pressure susceptibility, Cr A over M: {craoverm_fit:.4f} m^2/kg")
print(
    f"Fitted empirical NTW acceleration (m/s^2): "
    f"N={ntw_accel_fit[0]: .3e}, T={ntw_accel_fit[1]: .3e}, W={ntw_accel_fit[2]: .3e}"
)
print(
    f"Mean position error of fit over 24 hours: {np.mean(np.linalg.norm(perr, axis=1)):.3f} meters"
)

# Nelder-Mead outer optimization over CrA/m and constant NTW empirical acceleration ACCEL_SCALE = 1e-9 # m/s^2 per unit of scaled parameter def minfunc_full(params, state0, timearr, pgcrf, settings): craoverm = params[0] ntw_accel = np.array(params[1:4]) * ACCEL_SCALE if craoverm < 0 or craoverm > 1: return 1e30 sp = make_satprops(craoverm, ntw_accel) _, _, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp) return float(np.sum(perr**2)) x0 = np.array([0.01, 0.0, 0.0, 0.0]) r = minimize( minfunc_full, x0, args=(state0, timearr, pgcrf, settings), method="Nelder-Mead", options={"xatol": 1e-5, "fatol": 1e-4, "maxiter": 400}, ) craoverm_fit = r.x[0] ntw_accel_fit = np.array(r.x[1:4]) * ACCEL_SCALE # Final least-squares fit with optimized CrA/m and NTW empirical acceleration sp = make_satprops(craoverm_fit, ntw_accel_fit) state0, res, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp) print(f"Satellite radiation pressure susceptibility, Cr A over M: {craoverm_fit:.4f} m^2/kg") print( f"Fitted empirical NTW acceleration (m/s^2): " f"N={ntw_accel_fit[0]: .3e}, T={ntw_accel_fit[1]: .3e}, W={ntw_accel_fit[2]: .3e}" ) print( f"Mean position error of fit over 24 hours: {np.mean(np.linalg.norm(perr, axis=1)):.3f} meters" )

Satellite radiation pressure susceptibility, Cr A over M: 0.0232 m^2/kg
Fitted empirical NTW acceleration (m/s^2): N=-1.532e-10, T= 5.100e-10, W= 6.639e-12
Mean position error of fit over 24 hours: 0.405 meters

Validate: residuals vs SP3¶

Per-component position errors relative to the SP3 truth over the full 24-hour fit window.

In [6]:

# Plot position error
fig, ax = plt.subplots()
dates = [t.datetime() for t in timearr]
ax.plot(dates, perr[:, 0], label="$X$")
ax.plot(dates, perr[:, 1], label="$Y$")
ax.plot(dates, perr[:, 2], label="$Z$")
ax.set_xlabel("Time")
ax.set_ylabel("Position Error (m)")
ax.set_title("Propagation Error vs SP3 Truth for GPS Satellite")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.25), ncol=3)
fig.autofmt_xdate()
plt.tight_layout()
plt.show()

# Plot position error fig, ax = plt.subplots() dates = [t.datetime() for t in timearr] ax.plot(dates, perr[:, 0], label="$X$") ax.plot(dates, perr[:, 1], label="$Y$") ax.plot(dates, perr[:, 2], label="$Z$") ax.set_xlabel("Time") ax.set_ylabel("Position Error (m)") ax.set_title("Propagation Error vs SP3 Truth for GPS Satellite") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.25), ncol=3) fig.autofmt_xdate() plt.tight_layout() plt.show()

Fit vs. Free-Run: How Errors Grow¶

The least-squares fit above stayed at the meter level over its 24-hour window. But what happens if we hold the fitted initial state fixed and propagate beyond the fit interval?

Below we download a second day of SP3 truth data and propagate the already-fitted state0 across both days. The day-1 window is the fit (errors at the meter level); day 2 is pure prediction with no new measurements, and shows how quickly errors grow once we leave the fit interval.

In [7]:

# Download a second SP3 file (day after the first) and concatenate to get 2 days of truth data.
# Day 1 is 2023-364 (GPS week 2294). Day 2 is 2023-365 = Dec 31 2023 (Sunday), which starts GPS week 2295.
fname2 = "./ESA0OPSFIN_20233650000_01D_05M_ORB.SP3"
url2 = "http://navigation-office.esa.int/products/gnss-products/2295/ESA0OPSFIN_20233650000_01D_05M_ORB.SP3.gz"

if not os.path.exists(fname2):
    with urllib.request.urlopen(url2) as response:
        with open(fname2, "wb") as out_file:
            with gzip.GzipFile(fileobj=response) as uncompressed:
                out_file.write(uncompressed.read())

# Read day 2 and rotate to GCRF
pitrf2, timearr2 = read_sp3file(fname2)
pgcrf2 = np.stack(
    np.fromiter(
        (q * p for q, p in zip([sk.frametransform.rotation(sk.frame.ITRF, sk.frame.GCRF, t) for t in timearr2], pitrf2)),
        list,  # type: ignore
    ),
    axis=0,
)  # type: ignore

# Concatenate day 1 + day 2 for full 2-day truth.  The last point of day 1 and
# first point of day 2 are typically the same epoch, so drop the duplicate.
if (timearr2[0] - timearr[-1]).seconds == 0:
    timearr_full = np.concatenate((timearr, timearr2[1:]))
    pgcrf_full = np.concatenate((pgcrf, pgcrf2[1:, :]), axis=0)
else:
    timearr_full = np.concatenate((timearr, timearr2))
    pgcrf_full = np.concatenate((pgcrf, pgcrf2), axis=0)

# Propagation settings covering both days
settings2 = sk.propsettings(
    abs_error=1e-12,
    rel_error=1e-12,
    gravity_degree=10,
)
settings2.precompute_terms(timearr_full[0], timearr_full[-1])

# Propagate the already-fitted state0 (and sp) across the full 2-day span.
# Day 1 is the fit window; day 2 is pure prediction — no new measurements used.
res_full = sk.propagate(
    state0,
    timearr_full[0],
    timearr_full[-1],
    propsettings=settings2,
    satproperties=sp,
)

perr_full = np.zeros((len(timearr_full), 3))
for i, t in enumerate(timearr_full):
    perr_full[i, :] = res_full.interp(t)[0:3] - pgcrf_full[i, :]

perr_full_norm = np.linalg.norm(perr_full, axis=1)

# Split day 1 (fit window) vs day 2 (free run) for reporting
fit_end = timearr[-1]
is_fit = np.array([(t - fit_end).seconds <= 0 for t in timearr_full])
print(f"Mean |err| over day 1 (fit window):  {np.mean(perr_full_norm[is_fit]):8.3f} m")
print(f"Max  |err| over day 1 (fit window):  {np.max(perr_full_norm[is_fit]):8.3f} m")
print(f"Mean |err| over day 2 (free run):    {np.mean(perr_full_norm[~is_fit]):8.3f} m")
print(f"Max  |err| over day 2 (free run):    {np.max(perr_full_norm[~is_fit]):8.3f} m")

# Plot: components + error magnitude, shading the fit window vs free-run region
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 5), sharex=True)
dates_full = [t.datetime() for t in timearr_full]

ax1.plot(dates_full, perr_full[:, 0], label="$X$")
ax1.plot(dates_full, perr_full[:, 1], label="$Y$")
ax1.plot(dates_full, perr_full[:, 2], label="$Z$")
ax1.axvline(fit_end.datetime(), color="k", linestyle="--", linewidth=1)
ax1.axvspan(
    dates_full[0], fit_end.datetime(), alpha=0.1, color="tab:green", label="Fit window"
)
ax1.axvspan(
    fit_end.datetime(), dates_full[-1], alpha=0.1, color="tab:red", label="Free run"
)
ax1.set_ylabel("Position Error (m)")
ax1.set_title("Fitted-State Propagation Error vs SP3 Truth (Day 1 fit, Day 2 free-run)")
ax1.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=8)

ax2.semilogy(dates_full, perr_full_norm, color="tab:purple")
ax2.axvline(fit_end.datetime(), color="k", linestyle="--", linewidth=1)
ax2.axvspan(dates_full[0], fit_end.datetime(), alpha=0.1, color="tab:green")
ax2.axvspan(fit_end.datetime(), dates_full[-1], alpha=0.1, color="tab:red")
ax2.set_xlabel("Time")
ax2.set_ylabel("|Error| (m, log)")

fig.autofmt_xdate()
plt.tight_layout()
plt.show()

# Download a second SP3 file (day after the first) and concatenate to get 2 days of truth data. # Day 1 is 2023-364 (GPS week 2294). Day 2 is 2023-365 = Dec 31 2023 (Sunday), which starts GPS week 2295. fname2 = "./ESA0OPSFIN_20233650000_01D_05M_ORB.SP3" url2 = "http://navigation-office.esa.int/products/gnss-products/2295/ESA0OPSFIN_20233650000_01D_05M_ORB.SP3.gz" if not os.path.exists(fname2): with urllib.request.urlopen(url2) as response: with open(fname2, "wb") as out_file: with gzip.GzipFile(fileobj=response) as uncompressed: out_file.write(uncompressed.read()) # Read day 2 and rotate to GCRF pitrf2, timearr2 = read_sp3file(fname2) pgcrf2 = np.stack( np.fromiter( (q * p for q, p in zip([sk.frametransform.rotation(sk.frame.ITRF, sk.frame.GCRF, t) for t in timearr2], pitrf2)), list, # type: ignore ), axis=0, ) # type: ignore # Concatenate day 1 + day 2 for full 2-day truth. The last point of day 1 and # first point of day 2 are typically the same epoch, so drop the duplicate. if (timearr2[0] - timearr[-1]).seconds == 0: timearr_full = np.concatenate((timearr, timearr2[1:])) pgcrf_full = np.concatenate((pgcrf, pgcrf2[1:, :]), axis=0) else: timearr_full = np.concatenate((timearr, timearr2)) pgcrf_full = np.concatenate((pgcrf, pgcrf2), axis=0) # Propagation settings covering both days settings2 = sk.propsettings( abs_error=1e-12, rel_error=1e-12, gravity_degree=10, ) settings2.precompute_terms(timearr_full[0], timearr_full[-1]) # Propagate the already-fitted state0 (and sp) across the full 2-day span. # Day 1 is the fit window; day 2 is pure prediction — no new measurements used. res_full = sk.propagate( state0, timearr_full[0], timearr_full[-1], propsettings=settings2, satproperties=sp, ) perr_full = np.zeros((len(timearr_full), 3)) for i, t in enumerate(timearr_full): perr_full[i, :] = res_full.interp(t)[0:3] - pgcrf_full[i, :] perr_full_norm = np.linalg.norm(perr_full, axis=1) # Split day 1 (fit window) vs day 2 (free run) for reporting fit_end = timearr[-1] is_fit = np.array([(t - fit_end).seconds <= 0 for t in timearr_full]) print(f"Mean |err| over day 1 (fit window): {np.mean(perr_full_norm[is_fit]):8.3f} m") print(f"Max |err| over day 1 (fit window): {np.max(perr_full_norm[is_fit]):8.3f} m") print(f"Mean |err| over day 2 (free run): {np.mean(perr_full_norm[~is_fit]):8.3f} m") print(f"Max |err| over day 2 (free run): {np.max(perr_full_norm[~is_fit]):8.3f} m") # Plot: components + error magnitude, shading the fit window vs free-run region fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 5), sharex=True) dates_full = [t.datetime() for t in timearr_full] ax1.plot(dates_full, perr_full[:, 0], label="$X$") ax1.plot(dates_full, perr_full[:, 1], label="$Y$") ax1.plot(dates_full, perr_full[:, 2], label="$Z$") ax1.axvline(fit_end.datetime(), color="k", linestyle="--", linewidth=1) ax1.axvspan( dates_full[0], fit_end.datetime(), alpha=0.1, color="tab:green", label="Fit window" ) ax1.axvspan( fit_end.datetime(), dates_full[-1], alpha=0.1, color="tab:red", label="Free run" ) ax1.set_ylabel("Position Error (m)") ax1.set_title("Fitted-State Propagation Error vs SP3 Truth (Day 1 fit, Day 2 free-run)") ax1.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=8) ax2.semilogy(dates_full, perr_full_norm, color="tab:purple") ax2.axvline(fit_end.datetime(), color="k", linestyle="--", linewidth=1) ax2.axvspan(dates_full[0], fit_end.datetime(), alpha=0.1, color="tab:green") ax2.axvspan(fit_end.datetime(), dates_full[-1], alpha=0.1, color="tab:red") ax2.set_xlabel("Time") ax2.set_ylabel("|Error| (m, log)") fig.autofmt_xdate() plt.tight_layout() plt.show()

Mean |err| over day 1 (fit window):     0.405 m
Max  |err| over day 1 (fit window):     0.807 m
Mean |err| over day 2 (free run):       2.107 m
Max  |err| over day 2 (free run):       6.538 m

Comparing Integrators¶

The propagator supports several Runge-Kutta integrators of different orders, plus the Gauss-Jackson 8 multistep predictor-corrector. Higher-order Runge-Kutta methods can take larger steps for the same accuracy, so they often require fewer total function evaluations despite more stages per step. Gauss-Jackson 8 is specialised for smooth 2nd-order orbit ODEs and typically needs 3–10× fewer force evaluations than the adaptive methods at comparable accuracy — at the cost of a user-chosen fixed step size and no state-transition-matrix support.

Below we compare performance and accuracy of these integrators on the same problem, using the fitted initial state from above.

In [8]:

# Integrators to compare (from lowest to highest order).  Gauss-Jackson 8 is a
# fixed-step multistep predictor-corrector — it uses gj_step_seconds from
# propsettings rather than abs_error / rel_error.  120 s is a reasonable step
# for GPS (MEO) altitudes.
integrators = [
    ("rkts54         - Tsitouras 5(4)", sk.integrator.rkts54),
    ("rkv65          - Verner 6(5)", sk.integrator.rkv65),
    ("rkv87          - Verner 8(7)", sk.integrator.rkv87),
    ("rkv98          - Verner 9(8)", sk.integrator.rkv98),
    ("gauss_jackson8 - GJ8 (120 s)", sk.integrator.gauss_jackson8),
]

print(
    f"{'Integrator':<35} {'Steps':>6} {'Evals':>6} {'Time (ms)':>10} {'Max Err (m)':>12}"
)
print("-" * 75)

for name, integ in integrators:
    s = sk.propsettings(
        abs_error=1e-12,
        rel_error=1e-12,
        gravity_degree=10,
        integrator=integ,
        gj_step_seconds=120.0,  # only used by gauss_jackson8
    )
    s.precompute_terms(timearr[0], timearr[-1])

    t0 = time.perf_counter()
    result = sk.propagate(
        state0,
        timearr[0],
        timearr[-1],
        propsettings=s,
        satproperties=sp,
    )
    elapsed = (time.perf_counter() - t0) * 1e3

    # Compute max position error vs SP3 truth
    max_err = max(
        np.linalg.norm(result.interp(t)[0:3] - pgcrf[i, :])
        for i, t in enumerate(timearr)
    )

    stats = result.stats
    print(
        f"{name:<35} {stats.num_accept:>6} {stats.num_eval:>6} {elapsed:>10.1f} {max_err:>12.3f}"
    )

# Integrators to compare (from lowest to highest order). Gauss-Jackson 8 is a # fixed-step multistep predictor-corrector — it uses gj_step_seconds from # propsettings rather than abs_error / rel_error. 120 s is a reasonable step # for GPS (MEO) altitudes. integrators = [ ("rkts54 - Tsitouras 5(4)", sk.integrator.rkts54), ("rkv65 - Verner 6(5)", sk.integrator.rkv65), ("rkv87 - Verner 8(7)", sk.integrator.rkv87), ("rkv98 - Verner 9(8)", sk.integrator.rkv98), ("gauss_jackson8 - GJ8 (120 s)", sk.integrator.gauss_jackson8), ] print( f"{'Integrator':<35} {'Steps':>6} {'Evals':>6} {'Time (ms)':>10} {'Max Err (m)':>12}" ) print("-" * 75) for name, integ in integrators: s = sk.propsettings( abs_error=1e-12, rel_error=1e-12, gravity_degree=10, integrator=integ, gj_step_seconds=120.0, # only used by gauss_jackson8 ) s.precompute_terms(timearr[0], timearr[-1]) t0 = time.perf_counter() result = sk.propagate( state0, timearr[0], timearr[-1], propsettings=s, satproperties=sp, ) elapsed = (time.perf_counter() - t0) * 1e3 # Compute max position error vs SP3 truth max_err = max( np.linalg.norm(result.interp(t)[0:3] - pgcrf[i, :]) for i, t in enumerate(timearr) ) stats = result.stats print( f"{name:<35} {stats.num_accept:>6} {stats.num_eval:>6} {elapsed:>10.1f} {max_err:>12.3f}" )

Integrator                           Steps  Evals  Time (ms)  Max Err (m)
---------------------------------------------------------------------------
rkts54         - Tsitouras 5(4)       1429   8619       21.5        0.807
rkv65          - Verner 6(5)           627   6272       15.5        0.807
rkv87          - Verner 8(7)           160   2722        6.7        0.807
rkv98          - Verner 9(8)           112   2354        5.8        0.807
gauss_jackson8 - GJ8 (120 s)           720    777        1.9        0.807