6
$\begingroup$

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 different I can do 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

enter image description here

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – uhoh Jan 31 at 1:36
  • $\begingroup$ It would also be interesting, given two (or a set of) TLEs, to determine if an orbital maneuver had been performed. $\endgroup$ – rickhg12hs Jan 31 at 16:09
  • 1
    $\begingroup$ @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. $\endgroup$ – sudeep Jan 31 at 18:56
  • 1
    $\begingroup$ @rickhg12hs yes, the orbital maneuver is also one of the interesting points to determine from TLEs. $\endgroup$ – sudeep Jan 31 at 18:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.