sgp4 comparison

sgp4_comparison.py

#!/usr/bin/env python3
"""
This test will compare the output of ephem, skyfield and sgp4

Reference information from cbassa/sattools skymap
Time: 2019-10-31T18:22:00
Location: 52.8344N, 6.3785E, 10m
TLE:
USA 276
1 42689U 17022A   19300.77210760 0.00002000  00000-0  26008-4 0    09
2 42689  49.9898 242.6774 0014628 311.5070  48.5401 15.57936661    09

Results:
x:   +3209.1737 km; vx: +5.50221 km/s
y:   -2960.1539 km; vy: +5.33375 km/s
z:   +5174.5299 km; vz: -0.34825 km/s
r:     603.47 km; v:  -3.246 km/s
l:   2.05; b:  49.84; h:     392.17 km

R.A.:  19:11:13.88  Decl.:  10:59:40.6
R.A.:  19:10:17.79  Decl.:  10:57:39.8 (J2000)
Azi.: 225.6 Alt.: 40.1

Phase: 104.19
Magnitude: 5.38
"""
import ephem
import datetime
from skyfield.api import Topos, load
from skyfield.api import EarthSatellite
from sgp4.earth_gravity import wgs84
from sgp4.io import twoline2rv
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, Angle, CartesianRepresentation
from astropy.coordinates import GCRS, PrecessedGeocentric, ICRS

if __name__ == "__main__":
    # Set observer
    observer = ephem.Observer()
    observer.lon = "6.3785"
    observer.lat = "52.8344"
    observer.elevation = 10
    observer.date = ephem.date("2019/10/31 18:22:00")
    
    tle = ["USA 276",
           "1 42689U 17022A   19300.77210760 0.00002000  00000-0  26008-4 0    09",
           "2 42689  49.9898 242.6774 0014628 311.5070  48.5401 15.57936661    09"]
    satellite = ephem.readtle(tle[0], tle[1], tle[2])
    satellite.compute(observer)

    print("Ephem: ", satellite.ra, satellite.dec, satellite.range)

    satellite = EarthSatellite(tle[1], tle[2], tle[0])
    ts = load.timescale()
    t = ts.utc(2019, 10, 31, 18, 22, 0)
    location = Topos(latitude_degrees=52.8344, longitude_degrees=6.3785, elevation_m=0.0)
    ra, dec, dist = (satellite-location).at(t).radec(epoch='date')
    print("Skyfield: ", ra, dec, dist.km)
    print("GCRS: ", satellite.at(t).position.km)

    satellite = twoline2rv(tle[1], tle[2], wgs84)
    p, v = satellite.propagate(2019, 10, 31, 18, 22, 0)
    print("sgp4: ", p)

    t = Time("2019-10-31T18:22:00", format="isot", scale="utc")
    pos = SkyCoord(frame=PrecessedGeocentric(equinox=t, obstime=t),
                   x=p[0], y=p[1], z=p[2], unit=u.km, representation_type='cartesian')
    print("sgp4 + astropy: ", pos.gcrs.cartesian)