Skip to content

Instantly share code, notes, and snippets.

@cbassa
Created October 31, 2019 21:46
Show Gist options
  • Save cbassa/6c7c51c07ed662bfe0014044f05e1d63 to your computer and use it in GitHub Desktop.
Save cbassa/6c7c51c07ed662bfe0014044f05e1d63 to your computer and use it in GitHub Desktop.
sgp4 comparison
#!/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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment