TLE Fitting

Fit a TLE to a set of satellite states (position and velocity)

This uses Levenberg-Marquart non-linear least-squares fitting to tune TLE parameters to minimize the difference between the state positions and the SGP4-computed positions

Note that a TLE represents states in the TEME frame. Inputs are rotated into TEME frame from GCRF

# Imports
import satkit as sk
import numpy as np
import math as m

# Create a high-precision state
# Altitude for circular orbit
altitude = 450e3

# Radius & velocity
r0 = altitude + sk.consts.earth_radius
v0 = m.sqrt(sk.consts.mu_earth / r0)

# Inclination
inclination = 15 * m.pi / 180.0

# Create the state (3D position in meters, 3D velocity in meters / second)
state0 = np.array([r0, 0, 0, 0, v0 * m.cos(inclination), v0 * m.sin(inclination)])
# Make up an epoch
time0 = sk.time(2024, 3, 15, 13, 0, 0)

# Propagate the state forward by a day with high-precision propagator
res = sk.propagate(state0, time0, time0 + sk.duration(days=1.0))

# Get interpolated states every 10 minutes
times = [time0 + sk.duration(minutes=i) for i in range(0, 1440, 10)]
states = [res.interp(t) for t in times]

# Fit the TLE
(tle, fitresults) = sk.TLE.fit_from_states(states, times, time0 + sk.duration(days=0.5)) # type: ignore

# Print the result
print(tle)
print(fitresults['success'])

            TLE: none
                         NORAD ID: 00000,
                      Launch Year: 2000,
                            Epoch: 2024-03-16T01:00:00.000000Z,
                  Mean Motion Dot: 0 revs / day^2,
              Mean Motion Dot Dot: 0 revs / day^3,
                             Drag: -0.0001532958087389173,
                      Inclination: 14.991373023845723 deg,
                             RAAN: 355.99763657441315 deg,
                            eccen: 0.0013613550403152985,
                   Arg of Perigee: 199.31611943912216 deg,
                     Mean Anomaly: 62.01096415449389 deg,
                      Mean Motion: 15.408691923872402 revs / day
                            Rev #: 0

Convergence in parameter value

# Compute position errors (differences between TLE & state)

# Get the positions from sgp4
(pteme, vteme) = sk.sgp4(tle, times)
# Rotate positions from TEME to GCRF frame
pgcrf = [sk.frametransform.qteme2gcrf(t) * p for t, p in zip(times, pteme)]
# Take difference between state vector and SGP4 positions, and compute norm
pdiff = [p - s[0:3] for p, s in zip(pgcrf, states)]
pdiff = np.array([np.linalg.norm(p) for p in pdiff])


# Plot position errors
import plotly.graph_objects as go

fig = go.Figure()
fig.add_trace(go.Scatter(x=[t.datetime() for t in times], y=pdiff, mode='lines', name='Position Error',
                         line=dict(color='black', width=2)))
fig.update_layout(title='TLE Fitting Position Errors',
                  xaxis_title='Time',
                  yaxis_title='Position Error (m)')
fig.update_xaxes(showline=True, linewidth=2, linecolor="black", mirror=True)
fig.update_yaxes(showline=True, linewidth=2, linecolor="black", mirror=True)
fig.update_layout(
    xaxis=dict(
        gridcolor="#dddddd",
        gridwidth=1,
    ),
    yaxis=dict(
        gridcolor="#dddddd",
        gridwidth=1,
    ),
)

fig.show()