19
$\begingroup$

I have obtained a Two Line Element (TLE) of a satellite in Earth orbit from Celestrak at https://celestrak.org/NORAD/elements/ and I would like to use it to calculate an orbit.

How can I calculate the orbit from the TLE, and then plot it in 3D, using Python and the Skyfield package?

$\endgroup$

3 Answers 3

6
$\begingroup$

Code updated for python3*

def makecubelimits(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 = """1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993
2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
L1, L2 = TLE.splitlines()

from skyfield.api import Loader, EarthSatellite
from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')
data = load('de421.bsp')
ts   = load.timescale()

planets = load('de421.bsp')
earth   = planets['earth']

Roadster = EarthSatellite(L1, L2)

print(Roadster.epoch.tt)
hours = np.arange(0, 3, 0.01)

time = ts.utc(2018, 2, 7, hours)

Rpos    = Roadster.at(time).position.km
Rposecl = Roadster.at(time).ecliptic_position().km

print(Rpos.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:    
    fig = plt.figure(figsize=[10, 8])  # [12, 10]

    ax  = fig.add_subplot(1, 1, 1, projection='3d')

    x, y, z = Rpos
    ax.plot(x, y, z)
    for x, y, z in lons:
        ax.plot(x, y, z, '-k')
    for x, y, z in lats:
        ax.plot(x, y, z, '-k')

    centers, hw = makecubelimits(ax)

    print("centers are: ", centers)
    print("hw is:       ", hw)

    plt.show()

r_Roadster = np.sqrt((Rpos**2).sum(axis=0))
alt_roadster = r_Roadster - re

if True:
    plt.figure()
    plt.plot(hours, r_Roadster)
    plt.plot(hours, alt_roadster)
    plt.xlabel('hours', fontsize=14)
    plt.ylabel('Geocenter radius or altitude (km)', fontsize=14)
    plt.show()
$\endgroup$
0
18
$\begingroup$

SGP4 is the standard procedure that TLEs are intended to work with. All of the following are extremely helpful, but the most important point would be use a standard, recent SGP4 package that is recommended, do not try to use the elements in a TLE yourself. This is becuase the TLE and the SGP4 package are built for each other.

One interesting point is that for high altitude orbits with periods longer than 225 minutes, a modern SGP4 implementation will switch to an algorithm that is called SDP4. See the question “Deep space” corrections in SGP4; how does it account for the Sun's and Moon's gravity? for more on that.


The easiest to access SGP4 that I know is in the Python package Skyfield. You can find SGP4 routines in many languages in many places. I would recommend you choose something that is used widely and well maintained, and not just any random code on the internet.

Skyfield's SGP4 is based on :https://github.com/brandon-rhodes/python-sgp4

Here's a plot of the low earth orbit of the SpaceX F9 upper stage known as Roadster generated using Skyfield and Roadster's TLE, of course before the second burn that put it in heliocentric orbit.

Roadster's temporary LEO

def makecubelimits(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 = """1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993
2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
L1, L2 = TLE.splitlines()

from skyfield.api import Loader, EarthSatellite
from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')
data = load('de421.bsp')
ts   = load.timescale()

planets = load('de421.bsp')
earth   = planets['earth']

Roadster = EarthSatellite(L1, L2)

print(Roadster.epoch.tt)
hours = np.arange(0, 3, 0.01)

time = ts.utc(2018, 2, 7, hours)

Rpos    = Roadster.at(time).position.km
Rposecl = Roadster.at(time).ecliptic_position().km

print(Rpos.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:    
    fig = plt.figure(figsize=[10, 8])  # [12, 10]

    ax  = fig.add_subplot(1, 1, 1, projection='3d')

    x, y, z = Rpos
    ax.plot(x, y, z)
    for x, y, z in lons:
        ax.plot(x, y, z, '-k')
    for x, y, z in lats:
        ax.plot(x, y, z, '-k')

    centers, hw = makecubelimits(ax)

    print("centers are: ", centers)
    print("hw is:       ", hw)

    plt.show()

r_Roadster = np.sqrt((Rpos**2).sum(axis=0))
alt_roadster = r_Roadster - re

if True:
    plt.figure()
    plt.plot(hours, r_Roadster)
    plt.plot(hours, alt_roadster)
    plt.xlabel('hours', fontsize=14)
    plt.ylabel('Geocenter radius or altitude (km)', fontsize=14)
    plt.show()
$\endgroup$
9
  • 2
    $\begingroup$ Excellent answer showing the incredible power of Python and its packages! $\endgroup$
    – Uwe
    Mar 10, 2018 at 10:21
  • 1
    $\begingroup$ If you're interested, here's a modified and updated version for Python3 : pastebin.com/iuMPSkJq Thanks for the script! $\endgroup$ May 14, 2018 at 13:59
  • $\begingroup$ @EricDuminil that's great! It looks like you can delete one of the two occurrences of load() in your pastebin script, and that you've omitted the definition of load load = Loader(path) Also I'm not sure how to use the comments at the top, but they look intriguing. Usually all my Skyfield scripts need adding the extra parenthesis around tuples like here[f*np.pi for f in 0.5, 1, 2] Did anything else actually require changing for 2 to 3? $\endgroup$
    – uhoh
    May 14, 2018 at 14:11
  • 1
    $\begingroup$ @uhoh: My goal was to get a working Python3 script with the least amount of changes to your code. The comments at first are just indications on how to install the required packages in an Anaconda environment (anaconda.com/download). There's also a link to download the JPL ephemeris. load is defined directly in from skyfield.api import load. The only required changes for Python3 were parens around print arguments and tuples. $\endgroup$ May 14, 2018 at 15:17
  • $\begingroup$ @EricDuminil okay thanks! I should just update all my posted scripts in various Stack Exchange sites. I have some that were written with ancient (early) versions of Skyfield and it's been revamped significantly by the time 1.0 was released. $\endgroup$
    – uhoh
    May 14, 2018 at 21:57
4
$\begingroup$

The answer given before seems to work great, I am just here to give a different answer if you're more interested in learning about the software that goes into numerically propagation trajectories and 3D plotting

The general procedure to go about solving this problem would be:

  1. Read TLE. It will give you inclination, RAAN, eccentricity, argument of periapsis, mean anomaly, mean motion (revs/day), and some others.
  2. Calculate orbital period (T) in seconds (where n is mean motion)

$T=\frac{24.0 * 3600.0}{n}$

  1. Calculate semi-major axis

$sma = (\frac{T^2 \mu}{4\pi^2})^\frac{1}{3}$

  1. Calculate eccentric anomaly via Newton's root solving method (doesn't have to be Newton, but thats what I use). Since this is a transcendental equation, there is no algebraic solution for eccentric anomaly E, thus it is solved iteratively:

$M_{e} = E - esin(E) $

  1. Calculate true anomaly $\theta$:

$\theta = 2\space arctan(\sqrt{\frac{1+e}{1-e}}tan(\frac{E}{2})) $

  1. Propagate orbit for set amount of time using whichever perturbations you want
  2. Plot 3D

I show how to do steps 1-5 in a video on TLEs, using an ISS TLE example.

Step 6 is a bit more involved, and requires a bit of knowledge on ordinary differential equations (ODEs) and ODE solvers, which I cover in this series: https://www.youtube.com/playlist?list=PLOIRBaljOV8gn074rWFWYP1dCr2dJqWab

Another extension of this type of analysis would be to also calculate and plot the groundtracks with respect to Earth's surface. This is again more involved when it comes to software but you'll learn tons if you go through the process of whats going on under the hood of large Python packages: https://www.youtube.com/playlist?list=PLOIRBaljOV8ghaN9XcC4ubv-QLPOQByYQ

$\endgroup$
1
  • 1
    $\begingroup$ Thank you for the info, since I'm new to here (as you saw). I just updated the comment with the algorithm for calculating the COEs $\endgroup$ Mar 5, 2021 at 17:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.