How can I plot satellite's trajectory from three different TLEs to detect any deviation on path with time?

Question

I have obtained TLE information of a satellite from space-track.org for three different times. I am trying to study any deviation of the satellite trajectory from its usual path. Here I am comparing with three-time frames (current, 6 months ago, 1 year ago) with TLEs.

When I plot, I am getting what looks like a single path for three TLEs overlapped.

How can I study any deviation in these three TLEs visually and in more detail? Any suggestions on what I can do differently for the same.

I am using Python and the Skyfield package.

from skyfield.api import load, EarthSatellite
from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 
def plotaxislimit(axis, centers=None, hw=None):
    lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
    if centers == None:
        centers = [0.5*sum(pair) for pair in lims]
 
    if hw == None:
        widths  = [pair[1] - pair[0] for pair in lims]
        hw      = 0.5*max(widths)
        ax.set_xlim(centers[0]-hw, centers[0]+hw)
        ax.set_ylim(centers[1]-hw, centers[1]+hw)
        ax.set_zlim(centers[2]-hw, centers[2]+hw)
        print("hw was None so set to:", hw)
    else:
        try:
            hwx, hwy, hwz = hw
            print("ok hw requested: ", hwx, hwy, hwz)
            ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
            ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
            ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
        except:
            print("nope hw requested: ", hw)
            ax.set_xlim(centers[0]-hw, centers[0]+hw)
            ax.set_ylim(centers[1]-hw, centers[1]+hw)
            ax.set_zlim(centers[2]-hw, centers[2]+hw)
 
    return centers, hw
 
TLE_SAT01 = """1 28884U 05041A   20030.12392372  .00000082  00000-0  00000+0 0  9993
2 28884   0.0587 269.2229 0001764  37.9111  93.2791  1.00270805 52263"""
L1Sat01, L2Sat01 = TLE_SAT01.splitlines()

TLE_SAT02 = """1 28884U 05041A   19210.42688337  .00000078  00000-0  00000+0 0  9998
2 28884   0.0595 269.5353 0001937 216.2932 201.6080  1.00272640 50407"""
L1Sat02, L2Sat02 = TLE_SAT02.splitlines()

TLE_SAT03 = """1 28884U 05041A   19029.94520572  .00000068  00000-0  00000+0 0  9998
2 28884   0.0411 270.7112 0001828  43.2457  22.1983  1.00273963 48594"""
L1Sat03, L2Sat03 = TLE_SAT03.splitlines()

 
halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads = 180/pi, pi/180
 
data = load('de421.bsp')
ts   = load.timescale()
 
planets = load('de421.bsp')
earth   = planets['earth']
 
Sat01 = EarthSatellite(L1Sat01, L2Sat01, name ='Sat01')
Sat02 = EarthSatellite(L1Sat02, L2Sat02, name ='Sat02')
Sat03 = EarthSatellite(L1Sat03, L2Sat03, name ='Sat03')

print("Satellites Epoch Details::")
print("SAT01 MJD::", Sat01.epoch.tt)
print("SAT01 Date::",Sat01.epoch.utc_jpl())
print("SAT02 MJD::", Sat02.epoch.tt)
print("SAT02 Date::", Sat02.epoch.utc_jpl())
print("SAT03 MJD::", Sat03.epoch.tt)
print("SAT03 Date::", Sat03.epoch.utc_jpl())

hours = np.arange(0, 24, 0.30)
time1 = ts.utc(2020, 1, 30, hours)
time2 = ts.utc(2019, 7, 29, hours)
time3 = ts.utc(2019, 1, 29, hours)
 
Sat01pos    = Sat01.at(time1).position.km
Sat01posecl = Sat01.at(time1).ecliptic_position().km
print("Satellite 01 Position Shape Details::")
print(Sat01pos.shape)

Sat02pos    = Sat02.at(time2).position.km
Sat02posecl = Sat02.at(time2).ecliptic_position().km
print("Satellite 02 Position Shape Details::")
print(Sat02pos.shape)

Sat03pos    = Sat03.at(time3).position.km
Sat03posecl = Sat03.at(time3).ecliptic_position().km
print("Satellite 02 Position Shape Details::")
print(Sat03pos.shape)
 
re = 6378.
 
theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
lon0 = re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
                     lon0[1]*cph + lon0[0]*sph,
                     lon0[2]) )
    lons.append(lon)
 
lat0 = re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
    lats.append(lat)
 
if True:    
    figPlot = plt.figure(figsize=[12, 10])
    figPlot.suptitle('SATELLITE PROJECTION', fontsize=14, fontweight='bold')
    axdet  = figPlot.add_subplot(1, 1, 1, projection='3d')
 
    x, y, z = Sat01pos
    axdet.plot(x, y, z)
    #ax.text(8500, 500, 5000, Sat01.name)
    for x, y, z in lons:
        axdet.plot(x, y, z, '-k')
    for x, y, z in lats:
        axdet.plot(x, y, z, '-k')
        
    x, y, z = Sat02pos
    axdet.plot(x, y, z)
    #ax.text(5500, 500, 6500, Sat02.name)
    for x, y, z in lons:
        axdet.plot(x, y, z, '-k')
    for x, y, z in lats:
        axdet.plot(x, y, z, '-k')
        
    x, y, z = Sat03pos
    axdet.plot(x, y, z)
    #ax.text(5500, 500, 8500, Sat03.name)
    for x, y, z in lons:
        axdet.plot(x, y, z, '-k')
    for x, y, z in lats:
        axdet.plot(x, y, z, '-k')
 
    centers, hw = plotaxislimit(axdet)
 
    print("centers are: ", centers)
    print("hw is:       ", hw)
 
    plt.show()

Output from the code:

Satellites Epoch Details::
SAT01 MJD:: 2458878.6247244608
SAT01 Date:: A.D. 2020-Jan-30 02:58:27.0094 UT
SAT02 MJD:: 2458693.927684111
SAT02 Date:: A.D. 2019-Jul-29 10:14:42.7232 UT
SAT03 MJD:: 2458513.446006461
SAT03 Date:: A.D. 2019-Jan-29 22:41:05.7742 UT
Satellite 01 Position Shape Details::
(3, 80)
Satellite 02 Position Shape Details::
(3, 80)
Satellite 03 Position Shape Details::
(3, 80)
hw was None so set to: 46380.64771570716
centers are:  [-4.247959478401754, 6.076372600262403, 0.0]
hw is:        46380.64771570716

Answer 1

You're plotting the orbit of a geosynchronous satellite in active commercial use. At the scale you've drawn it, it has to overlap exactly, because variation smaller than you can see at this magnification would be big enough to cause customer service complaints. If the orbit weren't actively maintained, it would start to drift, but operational satellites maneuver regularly (several times per month) to stay in the little box assigned to them.

What is your end goal? Are you just doing this to satisfy your curiosity, or do you have a specific reason, like getting a degree in the field? If so, by all means, carry on! If not, you may wish to consider how much time and money you are willing to invest in getting an answer of sufficient quality.

I would advise subtracting, and looking at the whole curve in each component, not just an average, because many of the deviations you might see will be oscillating. First look for big things that ought to be conserved in Keplerian systems, like angular momentum and eccentricity vector. Then look at the classical orbital elements, osculating and mean separately. Save position and velocity for last, since you need to develop a sense of how fast those change even when everything is fine. However, assigning a particular meaning to any observed difference is tricky, because you're trying to do this with very sparse, low-quality data.

NORAD catalog number 28884 (aka 2005-041A) is Galaxy-15, an operational C-band telecom satellite in geosynchronous orbit. Such vehicles are naturally constrained -- synchrony demands exactly one specific period (86164 seconds), and therefore exactly one semi-major axis (exact value at the single-digit kilometer level depends on which perturbations you turned on). Due both to limitations on the precision with which parameters can be controlled or known, and to natural change in those parameters due to higher order terms in the gravity potential, tides, other massive bodies, and so on, all satellites' orbital parameters are constantly changing. If they don't regularly make small maneuvers, they drift out of their assigned orbit box into someone else's, which operators try very hard to avoid. If the inclination or eccentricity are too high, they'll drift out of their box for part of every day. If that happens north-south (inclination too high), it will only annoy their own customers; if it happens east-west (eccentricity too high) it will annoy their neighbors, too. If both are too high, it counts mainly as north-south, because if it crosses the longitude boundary only after it has gotten too high in latitude, then there probably isn't anyone there to interfere with anyway, though some commercial satellites do indeed exhibit significant daily motion.

Propagating a state for six months or a year is not easy, and demands good input data and a good solver. SGP4+TLE doesn't qualify. The reason that it is so popular is that it is standardized, reliable, easily available and TLEs are regularly generated for (nearly) all satellites in Earth orbit. Anyone who can get anything else uses that instead, because essentially everything else is better. SGP4 was obsolete decades ago, but it is still used because the inertia of government bureaucracy is astonishingly large. Also, of course, the US government actually uses something else internally -- TLEs are only for outsiders. I was very pleased to hear over a year ago that they had released an updated version with many advertised improvements, yet to my dismay neither I nor anyone else I know has been able to get their hands on the new "Type 4" TLEs that SGP4-XP needs as input.

Even the best propagator in the world can't recover facts you don't know about the maneuver schedule the vehicle followed during that year, unless you also have a large number of real raw measurements (because you are the vehicle's owner or operator, you have paid a commercial space situational awareness provider a sizable chunk of money, or you own your own telescope and time dedicated to the task -- which includes amateur astronomers, universities, and governments). With those, a good differential correction orbit determination engine can indeed calculate what maneuvers were performed at what times, but it is not easy or cheap; prices for these things start out around $30,000 and go up from there.