2
$\begingroup$

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?

$\endgroup$
3
  • 2
    $\begingroup$ There are at least two possible things going on. 1) Right now "At Celestrak, the current orbit is..." is unclear. Where exactly in Celestrak are you finding those numbers, can you link to the exact page? It's a big site! 2) The Earth's gravity is not exactly spherical so orbits can not be exactly elliptical. So perigee/apogees estimated from elliptical orbits inferred from mean orbital elements will not give correct numbers. 3) Note that those are pseudo-altitudes above Earth's surface, not distances from Earth's center, and the Earth's surface is not a sphere! $\endgroup$
    – uhoh
    Dec 20, 2021 at 22:27
  • 1
    $\begingroup$ Items 2 and 3 have been discussed in previous questions and answers already, so this may end up being partially a duplicate question, but right now it will be very helpful if you give more details on exactly where your exact numbers are coming from, because your question might turn out to be unique and in need of a new answer. $\endgroup$
    – uhoh
    Dec 20, 2021 at 22:29
  • 1
    $\begingroup$ I've added a link to the page I mentioned. $\endgroup$
    – Claudiu
    Dec 21, 2021 at 4:57

3 Answers 3

3
$\begingroup$

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.

$\endgroup$
1
  • 2
    $\begingroup$ @uhoh yes, i can reproduce them; that's the last paragraph. The way those terms are defined by SGP4, as mean rather than osculating quantities, 417 and 423 km, which the library outputs when you call the function that computes perigee and apogee, is the correct answer. No amount of trying to average the osculating altitude by hand will reproduce the same assumptions that went into the definitions of mean apogee and mean perigee. You have to average the same way NORAD does, just like Revisiting Space-Track 3 says. $\endgroup$
    – Ryan C
    Dec 24, 2021 at 15:24
0
$\begingroup$

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))
$\endgroup$
6
  • 1
    $\begingroup$ Shouldn't that be print((sat.alta + sat.radiusearthkm), (sat.altp + sat.radiusearthkm)) , assuming all referenced values are in kilometers? $\endgroup$
    – notovny
    Dec 22, 2021 at 20:46
  • 1
    $\begingroup$ It would be great if you added some more details to explain your syntax for users who are not familiar with this. $\endgroup$ Dec 22, 2021 at 21:38
  • 2
    $\begingroup$ @notovny No, it's correct with * rather than +. The docs at pypi.org/project/sgp4 say "Derived Orbit Properties These are computed when the satellite is first loaded, as a convenience for callers who might be interested in them. They aren’t used by the SGP4 propagator itself. a — Semi-major axis (earth radii). altp — Altitude of the satellite at perigee (earth radii, assuming a spherical Earth). alta — Altitude of the satellite at apogee (earth radii, assuming a spherical Earth)." $\endgroup$
    – Ryan C
    Dec 22, 2021 at 21:50
  • 2
    $\begingroup$ You could argue that there should be a +1 before the multiplication, but the original target stated was 417 and 423 km, not 6795 and 6801 km. $\endgroup$
    – Ryan C
    Dec 22, 2021 at 21:59
  • 1
    $\begingroup$ @RyanC Fair enough, objection withdrawn. $\endgroup$
    – notovny
    Dec 22, 2021 at 22:16
0
$\begingroup$

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.com/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.

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

plot of ISS various heights and altitudes

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()
$\endgroup$

Your Answer

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

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