High-Precision Propagation¶

This example demonstrates high-precision orbit propagation using satkit's advanced numerical integration capabilities. It propagates a GPS satellite orbit over a 24-hour period and validates the results against high-fidelity SP3 ephemeris data from the European Space Agency (ESA).

Overview¶

The example performs the following steps:

1. Loads Reference Data: Reads GPS satellite positions from an SP3 file containing precise orbit determination results (note: velocity is not included)
2. Rotate to Inertial Frame: Rotate the GPS satellite positions from the ITRF (Earth-fixed) to the GCRF (inertial) frame
3. Fits Initial Conditions: Use the high-precision propagator to output the state and state transition matrix at the SP3 timestamps. Use that matrix to linearize the system and solve for corrections to the initial state (position and velocity) that minimize errors over the 1-day data set. The least-squares fit minimizes: $||p_{gps}(t) - \hat{p}(t; s_0) - \Phi_p(t; s_0) \, \Delta s_0||^2$ where $p_{gps}(t)$ is the SP3 position in GCRF, $\hat{p}(t; s_0)$ is the propagated position from initial state $s_0$, $\Phi_p(t; s_0)$ is the position block of the state transition matrix, and $\Delta s_0$ is the correction to the initial state, which is what we are solving for.

Note: We wrap a non-linear solve around the least-squares fit to also recover the solar radiation pressure coefficient ($C_r A/m$) and — more interestingly — a constant empirical acceleration in the velocity-aligned NTW frame. This empirical acceleration is a catch-all for any un-modelled constant bias in the force model (mismodelled SRP shape, thermal re-radiation, antenna thrust, gravity truncation residuals, Earth-orientation errors, etc.). Fitting it is standard practice in operational GNSS orbit determination and significantly tightens the 24-hour residuals.

4. Validates Results: Compares the propagated positions against the SP3 truth data and visualizes the position errors
5. Fit vs. Free-Run: Downloads a second day of SP3 truth, fits the initial state using only day 1, then propagates across both days with the fit held fixed. The comparison against day 2 (pure prediction, no new measurements) shows how error grows once we leave the fit window.
6. Compares Integrators: Benchmarks several Runge-Kutta integrators and the Gauss-Jackson 8 multistep predictor-corrector on the same problem, reporting steps, force evaluations, wall time, and max position error vs SP3 truth.

This example showcases satkit's ability to achieve meter-level accuracy over extended propagation periods when properly configured with high-fidelity force models and optimized initial conditions.

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

In [1]:

# Necessary imports
import satkit as sk
import numpy as np
import matplotlib.pyplot as plt
import scienceplots  # noqa: F401

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

# Necessary imports import satkit as sk import numpy as np import matplotlib.pyplot as plt import scienceplots # noqa: F401 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)

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"
import os
import urllib.request
import gzip

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.qitrf2gcrf(timearr), pitrf)),
        list,  # type: ignore
    ),
    axis=0,
)  # type: ignore
# Crude estimation of initial velocity is just the difference between the 1st two position states divided by the difference
# in time.  This will have a very high error value, since the two points are 5 minutes apart.
vgcrf = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds

# Create the initial state
state0 = np.concatenate((pgcrf[0, :], vgcrf))

# Settings to use in propagations
# The integrator and gravity model can be set via constructor kwargs or properties.
# 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])


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 non-linear optimization over CrA/m plus a constant NTW-frame acceleration
# [aN, aT, aW].  Physically, this is a catch-all "empirical acceleration" that
# absorbs any *un-modelled* constant bias in the force model (mismodelled SRP
# shape, antenna thrust, thermal re-radiation, gravity-model truncation
# residuals, etc.).  The along-track (T) term is the most important — it maps
# directly to orbital energy / semi-major-axis drift.  At each outer step we
# do a linearized least-squares re-fit of 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.
from scipy.optimize import minimize

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 the mean error and fitted parameters
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"
)

# Plot position error
fig, ax = plt.subplots(figsize=(10, 5))
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()
fig.autofmt_xdate()
plt.tight_layout()
plt.show()

# 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" import os import urllib.request import gzip 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.qitrf2gcrf(timearr), pitrf)), list, # type: ignore ), axis=0, ) # type: ignore # Crude estimation of initial velocity is just the difference between the 1st two position states divided by the difference # in time. This will have a very high error value, since the two points are 5 minutes apart. vgcrf = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds # Create the initial state state0 = np.concatenate((pgcrf[0, :], vgcrf)) # Settings to use in propagations # The integrator and gravity model can be set via constructor kwargs or properties. # 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]) 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 non-linear optimization over CrA/m plus a constant NTW-frame acceleration # [aN, aT, aW]. Physically, this is a catch-all "empirical acceleration" that # absorbs any *un-modelled* constant bias in the force model (mismodelled SRP # shape, antenna thrust, thermal re-radiation, gravity-model truncation # residuals, etc.). The along-track (T) term is the most important — it maps # directly to orbital energy / semi-major-axis drift. At each outer step we # do a linearized least-squares re-fit of 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. from scipy.optimize import minimize 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 the mean error and fitted parameters 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" ) # Plot position error fig, ax = plt.subplots(figsize=(10, 5)) 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() fig.autofmt_xdate() plt.tight_layout() plt.show()

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

Fit Initial State and Validate¶

Download the SP3 file, rotate positions from ITRF to GCRF, then 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)$. This assumes the 6-DOF state is the only thing being estimated and converges in ~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 "un-modelled acceleration" that absorbs any constant bias left over after the gravity / third-body / SRP model is applied. Such biases come from things like mismodelled SRP box-wing geometry, antenna thrust, thermal re-radiation, truncation of the spherical-harmonic gravity expansion, or errors in Earth-orientation data. 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.).

The final plot shows per-component position errors relative to the SP3 truth data over the full 24-hour period.

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 [4]:

import time

# 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}"
    )

import time # 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       18.3        0.766
rkv65          - Verner 6(5)           628   6282       13.1        0.766
rkv87          - Verner 8(7)           160   2722        5.7        0.766
rkv98          - Verner 9(8)           112   2354        4.9        0.766
gauss_jackson8 - GJ8 (120 s)           720    777        1.7        0.766

Fit vs. Free-Run: How Errors Grow¶

The least-squares fit above was performed over a single 24-hour SP3 file, and errors stayed at the meter level over that 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, fit the initial state using only day 1, then propagate through both days. We compare the propagated trajectory to the SP3 truth across both days to visualize how errors grow once we leave the fit window.

In [5]:

# 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):
    import urllib.request
    import gzip

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

# New 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])

# Fit the initial state using ONLY day 1 data (timearr, pgcrf), with CrAoverM from the earlier fit.
# Re-seed velocity from the first two day-1 points.
v0_seed = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds
state0_day1 = np.concatenate((pgcrf[0, :], v0_seed))
state0_fit, _, perr_fit = linearized_least_squares_fit(
    state0_day1, timearr, pgcrf, settings, sp
)

# Now "free run" the fitted state over 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_fit,
    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=(10, 7), 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 left", ncol=2)

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): import urllib.request import gzip 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.qitrf2gcrf(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) # New 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]) # Fit the initial state using ONLY day 1 data (timearr, pgcrf), with CrAoverM from the earlier fit. # Re-seed velocity from the first two day-1 points. v0_seed = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds state0_day1 = np.concatenate((pgcrf[0, :], v0_seed)) state0_fit, _, perr_fit = linearized_least_squares_fit( state0_day1, timearr, pgcrf, settings, sp ) # Now "free run" the fitted state over 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_fit, 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=(10, 7), 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 left", ncol=2) 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.385 m
Max  |err| over day 1 (fit window):     0.766 m
Mean |err| over day 2 (free run):       2.007 m
Max  |err| over day 2 (free run):       6.437 m