Compare PyEphem and Skyfield modelled position to actual position

orbit_compare.py

# Sentinel 1a position data taken from https://qc.sentinel1.eo.esa.int 
# POD Precise Orbit Ephemerides	AUX_POEORB	(only last 366 days)
# 
# S1A_OPER_AUX_POEORB_OPOD_20151028T122426_V20151007T225943_20151009T005943
# 
# !! Coordinate system unknow, indicates fixed but not which one assume ITRF
#
# <OSV>
#   <TAI>TAI=2015-10-08T21:51:09.000000</TAI>
#   <UTC>UTC=2015-10-08T21:50:33.000000</UTC>
#   <UT1>UT1=2015-10-08T21:50:33.218710</UT1>
#   <Absolute_Orbit>+08067</Absolute_Orbit>
#   <X unit="m">3803127.125157</X>
#   <Y unit="m">-5968232.460514</Y>
#   <Z unit="m">22940.818265</Z>
#   <VX unit="m/s">-1353.493961</VX>
#   <VY unit="m/s">-823.241903</VY>
#   <VZ unit="m/s">7430.107470</VZ>
#   <Quality>NOMINAL</Quality>
# </OSV>

import numpy as np
import dateutil.parser as dtparse
from math import sqrt, sin, cos, radians

# Unit m
rObs = np.array([3803127.125157, -5968232.460514, 22940.818265])
# Unit m/s
vObs = np.array([-1353.493961, -823.241903, 7430.107470])
# UTC Time
tObs = dtparse.parse("2015-10-08T21:50:33.000000")

# Sentinel 1a TLE from https://www.space-track.org/#/tle
tle = """Sentinel 1a
1 39634U 14016A   15281.91006962 -.00000329  00000-0 -60181-4 0  9993
2 39634  98.1821 287.3784 0001263  95.2589 264.8781 14.59197669 80679
"""

# ephem calculations
import ephem
tlelines = tle.split("\n")
sentinel1a = ephem.readtle(tlelines[0], tlelines[1], tlelines[2])
tdiff = sentinel1a._epoch.datetime() - tObs
print "-"*80
print "Time diff between position measurement and TLE epoch (s): %s " % tdiff.total_seconds()
print "-"*80
# Compte SGP4 at time of obs
sentinel1a.compute(tObs)

# Get radius and angles needed to convert from sphere to cartesian
r = sentinel1a.elevation + ephem.earth_radius
theta = sentinel1a.sublong
phi = radians(90.0) - sentinel1a.sublat

x = r * cos(theta) * sin(phi)
y = r * sin(theta) * sin(phi)
z = r * cos(phi)

rMod = np.array([x, y, z])
rDiff = sqrt(np.sum((rObs - rMod)**2, axis=0))

print " pyephem ".center(80, '-')
print "Total error in position sqrt(xd^2+yd^2+zd^2)(m): %s " % rDiff
print "Difference between observered and computed radius (m): %s" % (sqrt(sum(rObs ** 2)) - sqrt(sum(rMod ** 2)))
print "-"*80


# Skyfield 
from skyfield import api, constants, timelib, sgp4lib
eph = api.load('de421.bsp')
earth = eph['earth']

sat = earth.satellite(tle)

jd = timelib.JulianDate(utc=tObs.replace(tzinfo=api.utc))
mod_gcrs = sat.gcrs(jd)
rGCRS = mod_gcrs.position.au * constants.AU_M
rDiffGCRS = sqrt(np.sum((rObs - rGCRS)**2, axis=0))

print " skyfield ".center(80, '-')
print "Total error in position in GCRS sqrt(xd^2+yd^2+zd^2)(m): %s " % rDiffGCRS
print "Difference between observered and computed radius (m): %s" % (sqrt(sum(rObs ** 2)) - sqrt(sum(rGCRS ** 2)))

rTEME, vTEME, error = sat._position_and_velocity_TEME_km(jd)
rTEME /= constants.AU_KM
vTEME /= constants.AU_KM
vTEME *= constants.DAY_S
rITRF, vITRF = sgp4lib.TEME_to_ITRF(jd.ut1, rTEME, vTEME)
rITRF = rITRF * constants.AU_M
rDiffITRF = sqrt(np.sum((rObs - rITRF)**2, axis=0))

print "Total error in position in ITRF sqrt(xd^2+yd^2+zd^2)(m): %s " % rDiffITRF
print "Difference between observered and computed radius (m): %s" % (sqrt(sum(rObs ** 2)) - sqrt(sum(rITRF ** 2)))

# Include corrections for polar moitions
arcminute = constants.DEG2RAD / 60.0
arcsecond = arcminute / 60.0
second = 1.0 / (24.0 * 60.0 * 60.0)

# Polar motion parameters taken from http://www.celestrak.com/SpaceData/eop19620101.txt
# 2015 10 08 57303  0.206685  0.303262  0.2200140  0.0013979 -0.098966 -0.013020  0.000007 -0.000001  36

xp = 0.206685 * arcsecond
yp = 0.303262 * arcsecond

rITRF, vITRF = sgp4lib.TEME_to_ITRF(jd.ut1, rTEME, vTEME, xp=xp, yp=yp)
rITRF = rITRF * constants.AU_M
rDiffITRF = sqrt(np.sum((rObs - rITRF)**2, axis=0))

print "Total error in position in ITRF + polar sqrt(xd^2+yd^2+zd^2)(m): %s " % rDiffITRF
print "Difference between observered and computed radius (m): %s" % (sqrt(sum(rObs ** 2)) - sqrt(sum(rITRF ** 2)))
print " skyfield ".center(80, '-')