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 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 Schwarzschild general-relativistic correction (IERS 2010 §10.3 Eq. 10.12). 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, Lense-Thirring / de Sitter relativistic terms, 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, JPL's GPS orbits, etc.).

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= 4.097e-12, T= 5.355e-10, W=-1.577e-10
Mean position error of fit over 24 hours: 0.376 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.376 m
Max  |err| over day 1 (fit window):     0.777 m
Mean |err| over day 2 (free run):       2.041 m
Max  |err| over day 2 (free run):       6.577 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       20.8        0.777
rkv65          - Verner 6(5)           628   6282       15.1        0.777
rkv87          - Verner 8(7)           160   2722        6.5        0.777
rkv98          - Verner 9(8)           111   2333        5.6        0.777
gauss_jackson8 - GJ8 (120 s)           720    777        1.9        0.777