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

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.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)]

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

• This is a really interesting question! You have three trajectories and you'd like to somehow visualize how they differ. You can't just subtract because the times don't line up. There are probably several ways, I'll try to post one later today. There are also things called Poincaré maps. Another thing you can look at is just plotting some orbital parameters from the raw TLEs as a function of time (especially if you have many more TLEs, similar to what I did here
– uhoh
Jan 31, 2020 at 1:36
• It would also be interesting, given two (or a set of) TLEs, to determine if an orbital maneuver had been performed. Jan 31, 2020 at 16:09
• @uhoh - thanks for the input, and very insightful work. I am also trying and working on it. I will share if something works out. Looking forward to your post. Jan 31, 2020 at 18:56
• @rickhg12hs yes, the orbital maneuver is also one of the interesting points to determine from TLEs. Jan 31, 2020 at 18:57
• @sudeep this is indeed really interesting; please could you elaborate on what you mean by "When I plot, I am getting what looks like a single path for three TLEs overlapped. " and what is it you are actually trying to achieve with the code and numrical/graphical outputs? i.e. what is the problem and what do you want to do. Nov 22, 2020 at 23:38

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.

• I like your answer; I made some small adjustments but please feel free to roll back or edit further. I think that SGP4 + TLEs is so popular because the TLEs are the only satellite data so readily available. Unless I'm mistaken, getting any better data on state vectors of specific satellites requires some combination of money, permission, and access to tracking or telemetry. Some scientific Earth observation satellites might be the exception; in their data products they may indeed have precise times and positions.
– uhoh
Dec 4, 2021 at 22:27
• – uhoh
Dec 4, 2021 at 22:41
• @uhoh I like your changes. The TLE database is a uniquely valuable resource -- I just wish they would actually let us use the updates they've been advertising. Dec 5, 2021 at 0:56
• – uhoh
Dec 5, 2021 at 0:59
• @uhoh Yes, those updates. There's now a link to version 8.2 on the space-track.org documentation page, but the link doesn't seem to work, so I'm holding off on recommending it until they fix it. Dec 5, 2021 at 1:03