Comparing the results of sgp4 and cysgp4

compare_sgp4_cysgp4.py

#!/usr/bin/env python3
import numpy as np
from sgp4.earth_gravity import wgs84
from sgp4.io import twoline2rv
from astropy.time import Time
import astropy.units as u
from cysgp4 import PyTle, PyObserver, propagate_many, Satellite, PyDateTime

if __name__ == "__main__":
    # TLE
    tle = ["USA 259",
           "1 40344U 14081A   19289.95767003 0.00000526  00000-0  00000-0 0    01",
           "2 40344  63.4408 358.0868 7184845 269.5721  15.5423  2.00643946    03"]

    # Time range, from epoch in 5 day steps until 80 days 
    dt = np.arange(0, 115201.0, 5*1440)
    t = Time(58772.95767003, format="mjd") + dt * u.min

    # Use pythonic SGP4
    satellite = twoline2rv(tle[1], tle[2], wgs84)
    p = []
    v = []
    for tm in t:
        pm, vm = satellite.propagate(tm.datetime.year,
                                     tm.datetime.month,
                                     tm.datetime.day,
                                     tm.datetime.hour,
                                     tm.datetime.minute,
                                     tm.datetime.second +
                                     tm.datetime.microsecond/1000000.0)
        p.append(pm)
        v.append(vm)

    pteme_sgp4 = np.asarray(p)
    vteme_sgp4 = np.asarray(v)

    print("pysgp4")
    for i in range(len(dt)):
        print("%17.8f %16.8f %16.8f %16.8f %14.8f %14.8f %14.8f" % (dt[i],
                                                                    pteme_sgp4[i, 0],
                                                                    pteme_sgp4[i, 1],
                                                                    pteme_sgp4[i, 2],
                                                                    vteme_sgp4[i, 0],
                                                                    vteme_sgp4[i, 1],
                                                                    vteme_sgp4[i, 2]))

    # Use cysgp4 propagate_many
    observer = np.array([PyObserver(0.0, 0.0, 0.0)])[np.newaxis, :, np.newaxis]
    tles = np.array([PyTle(*tle)])[np.newaxis, np.newaxis, :]
    mjds = t.mjd[:, np.newaxis, np.newaxis]

    result = propagate_many(mjds, tles, observer, on_error="coerce_to_nan", do_eci_pos=True, do_eci_vel=True, do_geo=False, do_topo=True)

    pteme_cysgp4 = result['eci_pos'].squeeze()
    vteme_cysgp4 = result['eci_vel'].squeeze()

    print("cysgp4 propagate_many")
    for i in range(len(dt)):
        print("%17.8f %16.8f %16.8f %16.8f %14.8f %14.8f %14.8f" % (dt[i],
                                                                    pteme_cysgp4[i, 0],
                                                                    pteme_cysgp4[i, 1],
                                                                    pteme_cysgp4[i, 2],
                                                                    vteme_cysgp4[i, 0],
                                                                    vteme_cysgp4[i, 1],
                                                                    vteme_cysgp4[i, 2]))

    # Use cysgp4 Satellite
    pytle = PyTle(*tle)
    pyobs = PyObserver(0.0, 0.0, 0.0)
    pydt = PyDateTime.from_mjd(t[-1].mjd)
    sat = Satellite(pytle, pyobs, pydt)

    pteme_cysgp4_single = np.asarray(sat.eci_pos().loc)
    vteme_cysgp4_single = np.asarray(sat.eci_pos().vel)

    print("cysgp4 Satellite")
    print("%17.8f %16.8f %16.8f %16.8f %14.8f %14.8f %14.8f" % (80 * 1440,
                                                                pteme_cysgp4_single[0],
                                                                pteme_cysgp4_single[1],
                                                                pteme_cysgp4_single[2],
                                                                vteme_cysgp4_single[0],
                                                                vteme_cysgp4_single[1],
                                                                vteme_cysgp4_single[2]))