sgp4 propagator accuracy

[trufanov-nok]

Hi, I'm new to julia and tried to reproduce my python calculations with this package. Unfortunately the results don't match with python's skyfield package.
I have 2 TLE records (attached). test.tle.zip Julia code:

using SatelliteToolbox
tles = read_tle("/home/truf/tests/test.tle")
o = init_orbit_propagator(Val(:sgp4), tles[1])
r, v = propagate!(o, 33208.95199608803)

result: ([-4.809400178348191e6, 4.756700552513811e6, 330.22513111666956], [675.0195082213636, 681.3381457850203, 7617.332409906825])

Python code:

from skyfield.api import EarthSatellite, load

line1 = "1 41386U 16016A   23049.16798291  .00026516  00000-0  30339-3 0  9998"
line2 = "2 41386  97.1815 134.9260 0010580 260.1146  99.8918 15.62057491388243"

line3 = "1 41386U 16016A   23049.55234578  .00023500  00000-0  26893-3 0  9994"
line4 = "2 41386  97.1814 135.3158 0010590 258.1661 101.8410 15.62074590388300"

ts = load.timescale()

satellite1 = EarthSatellite(line1, line2, '', ts)
satellite2 = EarthSatellite(line3, line4, '', ts)

# number of seconds
satellite2.epoch.utc_datetime().timestamp() - satellite1.epoch.utc_datetime().timestamp()

# propogate
geopos = satellite1.at(satellite2.epoch)

geopos.position.m
geopos.velocity.m_per_s

Results are:

>>> satellite2.epoch.utc_datetime().timestamp() - satellite1.epoch.utc_datetime().timestamp()
33208.95199608803

>>> geopos.position.m
array([-4784717.55451921,  4781517.86743171,    10831.07255969])
>>> geopos.velocity.m_per_s
array([ 695.52038204,  678.07536636, 7615.78105353])

They don't match.
Am I missing something or sgp4 implementatons of SatelliteToolbox and skyfield produce dfferent results? If so, which implementation is correct?

[ronisbr]

Hi @trufanov-nok !

Am I missing something or sgp4 implementatons of SatelliteToolbox and skyfield produce dfferent results? If so, which implementation is correct?

Both :) There are two differences between your code in SatelliteToolbox and in skyfield.

First, skyfield uses the constants from WGS72 by default whereas SatelliteToolbox uses from WGS84.

Second, the function at in skyfield provides values in GCRF reference frame:

The simplest form in which you can generate a satellite position is to call its at() method, which will return an (x,y,z) position relative to the Earth’s center in the Geocentric Celestial Reference System. (GCRS coordinates are based on even more precise axes than those of the old J2000 system.)

whereas propagate! in SatelliteToolbox uses the same reference frame as used to represent in the orbital elements. In the case of a TLE, this reference frame is TEME.

The code in SatelliteToolbox.jl to match skyfield's is:

using SatelliteToolbox

tles = tle"""
1 41386U 16016A   23049.16798291  .00026516  00000-0  30339-3 0  9998
2 41386  97.1815 134.9260 0010580 260.1146  99.8918 15.62057491388243
1 41386U 16016A   23049.55234578  .00023500  00000-0  26893-3 0  9994
2 41386  97.1814 135.3158 0010590 258.1661 101.8410 15.62074590388300
"""

orbp = init_orbit_propagator(Val(:sgp4), tles[1]; sgp4_gc = sgp4_gc_wgs72)

# Get the position and velocity in TEME reference frame.
r_teme, v_teme = propagate_to_epoch!(orbp, tles[2].epoch)

# Convert to GCRF. Notice that we will need EOP data. If this is not wanted, we
# can use J2000.
eop = get_iers_eop()

# Obtain the DCM that converts TEME to GCRF.
D_gcrf_teme = r_eci_to_eci(TEME(), GCRF(), tles[2].epoch, eop)

# Convert the values.
r_gcrf = D_gcrf_teme * r_teme
v_gcrf = D_gcrf_teme * v_teme

The result is:

julia> r_gcrf
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
    -4.784717548676036e6
     4.78151787311928e6
 10831.140628953564

julia> v_gcrf
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
  695.5204459152549
  678.0753200830255
 7615.78105182469

The minor differences are probably some differences in the implementation.

[trufanov-nok]

Indeed! Thanks a lot for explanation and the code sample!

[ronisbr]

You're welcome!

[trufanov-nok]

It seems only sgp4 propagator contains api to read orbit from TLE directly while docs mention other propagators should support it to... Was this api overload removed for them?

[ronisbr]

Yes! It is a problem in the docs. A TLE is a set of mean elements for SGP4. I removed that interface so that we can avoid some errors.