Failing at getting apogee and perigee from TLE

Question

At Celestrak, the current orbit for ISS is currently listed as 423x417 km, but if I'm getting the TLE data from Celestrak and use orb.get_lonlatalt(now) from pyorbital by iterating the orbit for the next period (I'm doing that by reading the altitudine every 30 seconds for the next 93 minutes in the future), I get a maximum altitude of 436 km instead of 423 km (the perigee is not off, at 417 km, exactly as specified on Celestrak website). What am I doing wrong and why my (apogee) data is different than on Celestrak?

Secondary question: what's the cleanest, fastest, most elegant way to get apogee and perigee (in Python) from TLE data?

Answer 1

As usual when we discuss TLE headaches (here, here, here, here, and here, just to name a few), the core issue is keeping track of whether you are working with mean elements or osculating ones. Within the SGP4 ecosystem, apogee and perigee are not osculating quantities, they are mean quantities, so their values cannot easily be derived from osculating altitude.

Instead, the proper comparison to make is to mean semi-major axis and mean eccentricity. Doing that, I find over the course of these four days (centered around the epoch 21356.70730882, to match the graphs in uhoh's answer, rather than 21356.62544795) apogee slowly and smoothly falls from 423.186 km to 423.020 km, and perigee slowly and smoothly falls from 417.009 km to 416.852 km, so they round to the same integers shown by Celestrak for the whole time.

I used the default, out of the box, U.S. Space Force SGP4 Python wrapper. Lines 308-310 of the demonstration script Sgp4Prop.py calls

     self.Sgp4Prop.Sgp4GetPropOut(key, self.PROPOUT_OSC_ELEM, oscElem)
     self.Sgp4Prop.Sgp4GetPropOut(key, self.PROPOUT_MEAN_ELEM, meanElem)
     self.Sgp4Prop.Sgp4GetPropOut(key, self.PROPOUT_NODAL_AP_PER, nodalApPer)

If you take $a$ and $e$ from oscElem[0:2], you get results that bounce around. If instead you take $a$ and $e$ from meanElem[0:2], you get results that are very slowly and smoothly changing, exactly as they were designed, when the averaging was done analytically, decades ago, when the theoretical model of which perturbations would be subtracted was chosen.

To compare among the variables, you could attempt to make the huge effort of reconstructing exactly what those perturbations are and reversing them, or you could just ask SGP4 to tell you mean elements instead of osculating ones.

nodalApPer[0:3] is mean motion, then apogee, then perigee. For the data set in question, using the mean element values, a(1+e) equals apogee and a(1-e) always equals perigee, to within 1 part in 2^52.

Answer 2

To answer my own question, I've finally found an elegant way of calculating apogee and perigee from TLE:

from sgp4.api import Satrec
line=['0 ISS (ZARYA)',
'1 25544U 98067A   21356.62544795  .00006800  00000-0  13125-3 0  9998',
'2 25544  51.6428 130.9420 0004657 342.5227  11.5462 15.49048823317794']
sat = Satrec.twoline2rv(line[1],line[2])
print((sat.alta * sat.radiusearthkm), (sat.altp * sat.radiusearthkm))

Answer 3

I'm happy to see OP has used Brandon Rhodes' sgp4 python library to obtain periapsis and apoapsis from a TLE. I've always used the skyfield python library for this which is from the same developer and builds on the sgp4 library but has way more features.

From https://celestrak.org/NORAD/elements/table.php?tleFile=stations&title=Space%20Stations&orbits=1&pointsPerRev=90&frame=1 as of right now:

Int. Des.  Cat No.    name   T (min) i (deg) Apo x Peri (km)  ecc  TLE age (day)
1998-067A   25544 ISS (ZARYA) 92.96   51.64  423 x 416    0.0004540  0.33

which cites the most recent TLE

ISS (ZARYA)             
1 25544U 98067A   21356.70730882  .00006423  00000-0  12443-3 0  9994
2 25544  51.6431 130.5342 0004540 343.5826 107.2903 15.49048054317816

What I'll do is propagate the satellite for +/- 2 days around it's epoch (time of best accuracy) and ask for the apo and peri from sgp4 and also plot the 'altitude' and the 'height' above Earth's surface at 30 second intervals to sample at a fine enough step size so that we can just use a simple numerical extrema detector to get pretty close to the perigees and apogees.

What's interesting is that these do not at all coincide with the maxima and minima of elevation above the Earth's surface! Since the orbit is very very close to circular, the Earth's equatorial bulge tends to force a second minimum altitude near the ascending and descending nodes, separate from the perigee.

What is even more interesting is that I can not reproduce Celestrak's 425 x 416!

For orbital altitudes that go into things like 423 x 416 km, we take the distances from Earth's geocenter and subtract 6378.137 kilometers; Earth's geoid's equatorial radius.

This is simply a convention.

• How is the altitude of a satellite defined, given that the Earth is not spherical?

We can also see the eccentricity is slowly increasing over a four day period centered on the TLE's epoch.

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load, EarthSatellite, Loader, wgs84
import datetime

loaddata = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
ts = loaddata.timescale() # include builtin=True if you want to use older files (you may miss some leap-seconds)
eph = loaddata('de421.bsp') # a small one, fine for this

satname, L1, L2 = """ISS (ZARYA)             
1 25544U 98067A   21356.70730882  .00006423  00000-0  12443-3 0  9994
2 25544  51.6431 130.5342 0004540 343.5826 107.2903 15.49048054317816""".splitlines()

sat = EarthSatellite(L1, L2, satname)

hw_days, dt_sec = 2.0, 10.

daysec = 24 * 3600 # seconds

earth_eq_radius = 6378.137 # km

days = np.arange(-hw_days * daysec, hw_days * daysec + 1, dt_sec) / daysec

times = ts.tt_jd(sat.epoch.tt + days)

g = sat.at(times) # geocentric position object

p = g.position.km  # geocentric positions (km)

d = g.distance().km

h = d - earth_eq_radius

el = wgs84.subpoint(g).elevation.km # height above the wgs84 geoid

def get_maxes_and_mins(a):
    """simple imperfect way to find points above or below both of those on either side
       and misses those where two are exactly equal as well as possible endpoints"""
    maxes = np.where((a[1:-1] > a[2:]) * (a[1:-1] > a[:-2]))[0]
    mins  = np.where((a[1:-1] < a[2:]) * (a[1:-1] < a[:-2]))[0]
    return maxes, mins

h_maxes, h_mins = get_maxes_and_mins(h)
el_maxes, el_mins = get_maxes_and_mins(el)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=[12, 6])

ax1.plot(days, h, label='height')
ax1.plot(days, el, label='elevation')

ax1.plot(days[h_maxes], h[h_maxes], '.k')
ax1.plot(days[h_mins], h[h_mins], '.k')

ax1.plot(days[el_maxes], el[el_maxes], '.k')
ax1.plot(days[el_mins], el[el_mins], '.k')

ax1.set_xlim(-0.5, 0.5)
ax1.set_ylim(410, 440)

ax1.legend()

ax2.plot(days[h_maxes], h[h_maxes], '-')
ax2.plot(days[h_mins], h[h_mins], '-')

ax2.plot(days[el_maxes], el[el_maxes], '-')
ax2.plot(days[el_mins], el[el_mins], '-')

ax2.set_xlim(-hw_days, hw_days)
ax2.set_ylim(412, 437)

epoch_string = datetime.datetime(*sat.epoch.utc[:5]).isoformat().replace('T', ' ') + 'UTC'
ax2.set_xlabel('days relative to TLE epoch')

plt.suptitle('ISS TLE epoch: ' + epoch_string)

plt.show()