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 bounded, non-linear solve for solar radiation pressure ($CrA/m$) around the least-squares fit to improve accuracy. This adds complexity but highlights the precision achievable with the high-fidelity model.

4. Validates Results: Compares the propagated positions against the SP3 truth data and visualizes the position errors

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 math as m
import numpy.typing as npt
from scipy.optimize import minimize_scalar
import plotly.graph_objects as go

# Necessary imports import satkit as sk import numpy as np import math as m import numpy.typing as npt from scipy.optimize import minimize_scalar import plotly.graph_objects as go

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
if not os.path.exists(fname):
    import urllib.request
    import gzip
    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
settings = sk.propsettings()
# Only compute sun, moon positions and earth rotation vectors once for all propagations
settings.precompute_terms(timearr[0], timearr[-1])
settings.abs_error = 1e-10
settings.rel_error = 1e-12
settings.gravity_order = 10

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

# Wrap linearized solve inside a bounded non-linear optimization of CrAoverM
def minfunc(v, state0, timearr, pgcrf, settings):
    sp = sk.satproperties_static(craoverm = v)
    _, _, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp)
    return np.sum(np.sum(perr ** 2, axis=1), axis=0)

# Non-linear minimization of CrAoverM, with linearized least squares fit of initial state at each step
v0 = 0.01
r = minimize_scalar(
    lambda v: minfunc(v, state0, timearr, pgcrf, settings),
    bounds=(0, 1),
    method='bounded',
  )

# Final least squares fit with optimized CrAoverM
sp = sk.satproperties_static(craoverm = r.x)
state0, res, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp)

# print the mean error
print(f"Satellite radiation pressure susceptibility, Cr A over M: {r.x:.4f} m^2/kg")
print(f"Mean position error of fit over 24 hours: {np.mean(np.linalg.norm(perr, axis=1)):.3f} meters")

# Plot position error
fig = go.Figure()
# Add markers for each component of the position error
fig.add_trace(
    go.Scatter(x=[t.datetime() for t in timearr], y=perr[:, 0], mode="lines", name="X")
)
fig.add_trace(
    go.Scatter(x=[t.datetime() for t in timearr], y=perr[:, 1], mode="lines", name="Y")
)
fig.add_trace(
    go.Scatter(x=[t.datetime() for t in timearr], y=perr[:, 2], mode="lines", name="Z")
)
fig.update_layout(
    title="Propagation Error vs SP3 Truth for GPS Satellite",
    title_font_size=18,
    xaxis_title="Time",
    yaxis_title="Position Error (m)",
    xaxis_title_font=dict(size=16),
    yaxis_title_font=dict(size=16),
    legend_font=dict(size=14),
    font=dict(size=14),
)
fig.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 if not os.path.exists(fname): import urllib.request import gzip 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 settings = sk.propsettings() # Only compute sun, moon positions and earth rotation vectors once for all propagations settings.precompute_terms(timearr[0], timearr[-1]) settings.abs_error = 1e-10 settings.rel_error = 1e-12 settings.gravity_order = 10 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 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 # Wrap linearized solve inside a bounded non-linear optimization of CrAoverM def minfunc(v, state0, timearr, pgcrf, settings): sp = sk.satproperties_static(craoverm = v) _, _, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp) return np.sum(np.sum(perr ** 2, axis=1), axis=0) # Non-linear minimization of CrAoverM, with linearized least squares fit of initial state at each step v0 = 0.01 r = minimize_scalar( lambda v: minfunc(v, state0, timearr, pgcrf, settings), bounds=(0, 1), method='bounded', ) # Final least squares fit with optimized CrAoverM sp = sk.satproperties_static(craoverm = r.x) state0, res, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp) # print the mean error print(f"Satellite radiation pressure susceptibility, Cr A over M: {r.x:.4f} m^2/kg") print(f"Mean position error of fit over 24 hours: {np.mean(np.linalg.norm(perr, axis=1)):.3f} meters") # Plot position error fig = go.Figure() # Add markers for each component of the position error fig.add_trace( go.Scatter(x=[t.datetime() for t in timearr], y=perr[:, 0], mode="lines", name="X") ) fig.add_trace( go.Scatter(x=[t.datetime() for t in timearr], y=perr[:, 1], mode="lines", name="Y") ) fig.add_trace( go.Scatter(x=[t.datetime() for t in timearr], y=perr[:, 2], mode="lines", name="Z") ) fig.update_layout( title="Propagation Error vs SP3 Truth for GPS Satellite", title_font_size=18, xaxis_title="Time", yaxis_title="Position Error (m)", xaxis_title_font=dict(size=16), yaxis_title_font=dict(size=16), legend_font=dict(size=14), font=dict(size=14), ) fig.show()

Satellite radiation pressure susceptibility, Cr A over M: 0.0233 m^2/kg
Mean position error of fit over 24 hours: 0.538 meters

Fit Initial State and Validate¶

Download the SP3 file, rotate positions from ITRF to GCRF, then iteratively fit the initial state (position + velocity) using linearized least squares with the state transition matrix. A bounded scalar optimization over the solar radiation pressure coefficient ($C_r A/m$) wraps the linear solve to improve accuracy.

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