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 Inerital Frame: Rotate the GPS satellite positions from the ITRF (Earth-fixed) to the GCRF (Inertial) frame

  3. Fits Initial Conditions: Uses numerical optimization to determine the initial velocity and solar radiation pressure coefficient (Cr·A/m) that best match the reference trajectory by integrating the initial state and minimizing the difference between the predicted and reported position

  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

[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
import plotly.graph_objects as go

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

[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
vgcrf = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds

# Initial state for non-linear least squares is initial velocity
# and susceptibility to radiation pressuer : Cr A / m
# v0 = np.array([vgcrf[0], vgcrf[1], vgcrf[2], 0.02])

# Skip ahead to the right answer... this takes a long time to run
v0 = np.array([2.47130555e03, 2.94682777e03, -5.34171918e02, 2.13018578e-02])


# Function to minimize ... takes the sum squared error between
# propagated position and SP3 truth position
def minfunc(v):
    tstart = timearr[0]
    tend = timearr[-1]
    settings = sk.propsettings()
    settings.gravity_order = 10 # type: ignore
    settings.enable_interp = True # type: ignore
    satprops = sk.satproperties_static(craoverm = v[3])

    res = sk.propagate(
        np.concatenate((pgcrf[0, :], v[0:3])),
        begin=tstart,
        end=tend,
        propsettings=settings,
        satproperties=satprops,
    )
    pest = np.zeros((len(timearr), 3))
    for i in range(len(timearr)):
        pest[i, :] = res.interp(timearr[i])[0:3]

    return np.sum(np.sum((pest - pgcrf) ** 2, axis=1), axis=0)


# Minimize difference between true satellite state and propagated state
# by tuning initial velocity and CrAoverM
r = minimize(minfunc, v0, method="Nelder-Mead")

# Propagate with optimized initial conditions
settings = sk.propsettings()
satprops = sk.satproperties_static(craoverm = r.x[3])
res = sk.propagate(
    np.concatenate((pgcrf[0, :], r.x[0:3])),
    begin=timearr[0],
    end=timearr[-1],
    propsettings=settings,
    satproperties=satprops,
)
perr = np.zeros((len(timearr), 3))
for i in range(len(timearr)):
    perr[i, :] = res.interp(timearr[i])[0:3] - pgcrf[i, :]
# print the mean error
print("Mean position error (m): ", np.mean(np.linalg.norm(perr, axis=1)))

# Plot position error
fig = go.Figure()

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",
    xaxis_title="Time",
    yaxis_title="Position Error (m)",
)
fig.show()

Mean position error (m):  0.872184186566293