2

I'm trying to calculate the perigee and apogee of a given satellite based on its TLE

from skyfield.api import Topos, EarthSatellite
ts = load.timescale()

# latest TLE as of this morning, epoch is epoch=2020-08-24T12:30:01Z
line=['0 ISS (ZARYA)',
      '1 25544U 98067A   20237.52084486  .00016717  00000-0  10270-3 0  9031', 
      '2 25544  51.6430  10.2947 0001353  63.6269 296.5020 15.49179055  2606']
satellite = EarthSatellite(line[1], line[2], line[0], ts)
t = ts.utc(2020, 8, 24, 12, range(30,123))  #epoch + a full orbit
geocentric = satellite.at(t)
subpoint = geocentric.subpoint()
print(f"max {max(subpoint.elevation.km)}")
print(f"min {min(subpoint.elevation.km)}")

This produces an apogee of 437.7 and perigee of 418.5. The perigee looks right, but apogee looks about 17km too high.

Thinking I was reading the docs wrong, I've also tried calculating the topocentric distance along the way and got identical results (to 7 places)

difference = satellite - bluffton
topocentric = difference.at(t)
alt, az, distance = topocentric.altaz()
print(f"max {max(distance.km)}")
print(f"min {min(distance.km)}")

Which produces identical results within 7 decimals.

Doing this more by hand, at the the TLE's epoch:

revs_per_day = 15.49179055
eccentricity = 0.0001353
earth_equatorial_radius = 6378.14
period_hrs = 24.0 / revs_per_day
range = (6028.9 * (period_hrs * 60))** (2 / 3)
perigee = range * (1 + eccentricity) - earth_equatorial_radius

produces an apogee of 420.02 km and perigee of 418.18 km, which is more in line with what I expected.

What am I doing wrong? Am I not understanding what Skyfield's distances represent?

rtphokie
  • 609
  • 1
  • 6
  • 14

1 Answers1

0

Correct, you have a slight misconception about the distance you are asking for; and then there’s some other discrepancy that I am less sure of.

  1. The “subpoint” beneath a body is the position on the Earth’s surface, where that surface is modeled as the WGS-84 oblate spheroid with several km more fatness at the equator than at the poles. So instead of using the subpoint.elevation.km which measures against that rising and falling surface as the ISS crosses from high to low latitude, try just geocentric.distance().km. You can always then subtract the Earth’s mean radius from that if you are less interested in the full shape of the orbit from the Earth’s center, than in just that component that’s above ground.
  2. I am not well enough versed in the math to be sure why your back-of-the-envelope calculation is off. But we can at least compare the maximum and minimum of geocentric.distance().km with, say, NASA HORIZONS output for comparison. I attach some of its output below: it shows the ISS rising and falling from 6788 km to 6803 km over an orbit, almost the same distances as the 6789–6802 km range returned from Skyfield — suggesting that the range in altitude is indeed more than 2 km for the current orbit.
*******************************************************************************
Ephemeris / WWW_USER Mon Aug 24 17:36:17 2020 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: International Space Station (spacecraft) (-125544) {source: iss}
Center body name: Earth (399)                     {source: DE431mx}
Center-site name: GEOCENTRIC
*******************************************************************************
Start time      : A.D. 2020-Aug-25 00:00:00.0000 UT      
Stop  time      : A.D. 2020-Aug-25 03:00:00.0000 UT      
Step-size       : 1 minutes
*******************************************************************************
...
Center pole/equ : High-precision EOP model        {East-longitude positive}
Center radii    : 6378.1 x 6378.1 x 6356.8 km     {Equator, meridian, pole}    
Target primary  : Earth
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE431mx}
Rel. light bend : Sun, EARTH                      {source: DE431mx}
Rel. lght bnd GM: 1.3271E+11, 3.9860E+05 km^3/s^2                              
Atmos refraction: NO (AIRLESS)
RA format       : HMS
Time format     : CAL 
EOP file        : eop.200824.p201115                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2020-AUG-24. PREDICTS-> 2020-NOV-14
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate (     0.0=NO )                           
****************************************************************************************************************
 Date__(UT)__HR:MN     R.A._____(ICRF)_____DEC   APmag   S-brt            delta      deldot    S-O-T /r    S-T-O
****************************************************************************************************************
$$SOE
 2020-Aug-25 00:00     11 24 02.95 +19 59 34.2    n.a.    n.a. 0.00004542337255   0.0095997  18.8151 /T 161.1844
...
 2020-Aug-25 00:25     16 34 19.41 -47 50 11.1    n.a.    n.a. 0.00004547799868   0.0003566 101.0263 /T  78.9709
 2020-Aug-25 00:26     16 56 12.33 -49 10 57.0    n.a.    n.a. 0.00004547808220   0.0000571 104.6886 /T  75.3086
 2020-Aug-25 00:27     17 19 10.96 -50 15 22.9    n.a.    n.a. 0.00004547804847  -0.0002273 108.3437 /T  71.6534
 2020-Aug-25 00:28     17 43 04.17 -51 02 05.5    n.a.    n.a. 0.00004547790303  -0.0004991 111.9896 /T  68.0075
...
 2020-Aug-25 01:13     05 07 59.71 +49 48 30.6    n.a.    n.a. 0.00004538023740  -0.0018603  73.4466 /L 106.5505
 2020-Aug-25 01:14     05 31 35.95 +50 43 36.4    n.a.    n.a. 0.00004537965836  -0.0010181  69.7785 /L 110.2186
 2020-Aug-25 01:15     05 55 59.08 +51 20 10.6    n.a.    n.a. 0.00004537942066  -0.0001615  66.1208 /L 113.8763
 2020-Aug-25 01:16     06 20 51.01 +51 37 19.0    n.a.    n.a. 0.00004537952843   0.0007013  62.4763 /L 117.5208
 2020-Aug-25 01:17     06 45 50.83 +51 34 35.0    n.a.    n.a. 0.00004537998258   0.0015622  58.8484 /L 121.1487
 2020-Aug-25 01:18     07 10 36.81 +51 12 03.0    n.a.    n.a. 0.00004538078075   0.0024132  55.2411 /L 124.7560
Brandon Rhodes
  • 83,755
  • 16
  • 106
  • 147
  • 1
    The above makes a lot of sense. The ISS was my starting point because it has an orbit I'm more familiar with. Ultimately I'm going to be running this on the entire space-track catalog (minus the deorbited stuff of course) for a rough census of sorts of satellites and debris.| With that use case in mind, geocentric.distance().km will do just fine. Thanks as always for the quick and thorough response. – rtphokie Aug 26 '20 at 03:22