comparison with SatTrack SGP4/SDP4 solution?

[jmfriedt]

I am looking into understanding the behavior of a geostationary satellite around its equilibrium point. As a great fan of SatTrack [1] for decades, this has been my first obvious choice, but the lack of resolution and the suspicious "*" symbols when computing the solution led me to look for an alternative (I have always used SatTrack successfully for LEO satellites, first time I am looking into geostationary). I have hence compared this Python library with the SatTrack output and observe some significant difference. Obviously I am using the TLE well beyond their expiry date, but am wondering why, as both libraries seem to be relying on the same SPG4 (or SPD4 for geostationary bodies) solver, the solutions are diverging. Has any benchmark been assessed between solver implementations? FYI, the TLE is from Telstar11N downloaded from Celestrak at two dates, namely March 27 2022 and April 18 2022 and am comparing how quickly the solutions diverge. Thanks and the associated software for generating the satellite range/azimuth/elevation with the Python library

#1 34111U 09009A   22086.18243766 -.00000270  00000+0  00000+0 0  9991
#2 34111   0.0189 162.3771 0001463 203.1699 207.1412  1.00270119 47943
TLE = """1 34111U 09009A   22086.18243766 -.00000270  00000+0  00000+0 0  9991
2 34111   0.0189 162.3771 0001463 203.1699 207.1412  1.00270119 47943"""

L1, L2 = TLE.splitlines()

from skyfield.api import Loader, EarthSatellite
from skyfield.api import N,S,E,W, wgs84
from skyfield.timelib import Time
import numpy as np
import csv

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')
data = load('de421.bsp')
ts   = load.timescale()

planets = load('de421.bsp')
earth   = planets['earth']
lab     = wgs84.latlon(+47.0*N, +6.0*E,elevation_m=143)

telstar11n = EarthSatellite(L1, L2)
hours = np.arange(0, 750, 0.1)
time = ts.utc(2022, 4, 12, hours)
diff=telstar11n-lab
Rpos    = diff.at(time).position.km
Rposecl = diff.at(time).ecliptic_position().km
topocentric=diff.at(time)
alt, az, distance = topocentric.altaz()
print(distance.km)

re = 6378.

r_Roadster = np.sqrt((Rpos**2).sum(axis=0))
alt_roadster = r_Roadster - re
with open('dist2703.csv', 'w', newline='') as csvfile:
    o=csv.writer(csvfile, delimiter=' ', quotechar='|', quoting=csv.QUOTE_MINIMAL)
    o.writerow(distance.km)
with open('alt2703.csv', 'w', newline='') as csvfile:
    o=csv.writer(csvfile, delimiter=' ', quotechar='|', quoting=csv.QUOTE_MINIMAL)
    o.writerow(alt.degrees)
with open('az2703.csv', 'w', newline='') as csvfile:
    o=csv.writer(csvfile, delimiter=' ', quotechar='|', quoting=csv.QUOTE_MINIMAL)
    o.writerow(az.degrees)

[1] https://ieeexplore.ieee.org/document/499660/similar whose source code I still have from https://sources.debian.org/src/sattrack/3.1.6-9/ after applying a Y2K bug patch

[brandon-rhodes]

@jmfriedt — No, I'm not aware that anyone has tried to sit down and figure out the difference between the results from the official SGP4 routine, and the results from SatTrack, which I had not myself heard of before. Its results look to be somewhat discrete; I see stair-steps in your graph rather than smooth curves?

Currently the focus of the sgp4 library is agreement with the official test suite of SGP4 itself; I'm not aware of anyone trying to do a thorough comparison between its results and those of any third-party libraries.

Did the folks who wrote SatTrack do their own propagation? The source code—

https://sources.debian.org/src/sattrack/3.1.6-9/src/sattrack/satprop.c/

—doesn't read like the SGP4 source code itself, which is unusual; most attempts at satellite propagation seem to have simply cut-and-pasted the SGP4 code, sometimes with modifications, but the SatTrack code looks in several ways different. Then again, I guess sgp4.c from PyEphem also looks different, so maybe there are more independent implementations of SGP4 out there than I realize!

[brandon-rhodes]

@jmfriedt — I wanted to check back and see if you learned any more by comparing the two satellite solutions; or whether you maybe compared both Skyfield and SatTrack to some other, third solution to see if it agreed more closely with one of them or the other. Let me know if you run across something Skyfield should be doing differently. Thanks!

[CelesTrak]

@jmfriedt and @brandon-rhodes - CelesTrak is actually revisiting this issue as we are in the process of converting our SGP4 code to work with GP data instead of TLEs. Details of new GP data are located at: https://celestrak.org/NORAD/documentation/gp-data-formats.php

In our research and review of various SGP4 code, we have found that nearly all very in some fashion or another. This is an issue with data standards and limitations of the legacy TLE format. Thus, the CelesTrak initiatives to migrate to GP format. Deviations or and/or errors often occur related to the conversation process needed to work within TLE data limitations. This can often generate errors in handling Leap Seconds, Julian Dates, etc.

When comparing two different SGP4 processes it is more important that they "agree", not that they nessesarily "match". For comparison CelesTrak recommends the use of distance measurements (e.g. mm, m, km) as opposed to other metrics such as percentages. If the two processes vary on the order of kilometers, something is wrong with the way one process is propagating. For example, a validation we recently performed between two different SGP4 propagators was performed across the entire satellite catalogue and the overall variation was millimeters. A paper describing our processes is located at: https://celestrak.org/publications/AIAA/2006-6753/

Our Best,

CelesTrak https://celestrak.org 501(c)3 non-profit

[rmathar]

It seems that the wiggles in the plots from the skyfielders/Python code are calculated, not just a plotting artefact, because the 2nd LTE lines tells us that the satellite has 1.00270 revolutions per day and there are approximately 10 wiggles in 10 days in the plots.

The major difference in the two codes is that the StaTrack solutions start in a sort of quadratic/parabolic deviation from their starting date as a function of time, whereas the Python code has a linear contribution. In the TLE the second derivative and drag term of the mean motion are zero, so the first derivative ought to dominate.

The thing to focus on is the deviation of the alt/az coordinates of the order of 0.2 deg right from the start for the data of 27 Mar 22, no propagation in time yet. (This is much larger than the tens of arcseconds expected from difference in handling atmospheric refraction.) So there seems to be disagreement in the codes what 'azimuth' and 'elevation' actually mean.

Another observation: for the tracking starting at 27 mar 22 the positions agree after 27 days (sidereal month). This might be related to the fact that in the sample skyfielders code the JPL ephemerides 'earth' is loaded, which includes the Earth's oscillation around the earth-moon barycenter, where other codes might look at the TLE's as just describing motion around the Earth's barycenter.

[brandon-rhodes]

@CelesTrak — Thank you for your July comment pointing out that different propagators can still be expected to agree to within millimeters, and that the difference in distances reported by @jmfriedt indicate a serious divergence in opinion about where a satellite is.

@rmathar — It doesn't look like the earth object is used anywhere in @jmfriedt’s code after being created, so I suspect it was just cut-and-pasted in as part of another example? I don't think it can affect the outcome here since that object isn't used in the calculation. But I agree with your other points, and they make good sense of several of the details in the plots.

I have the feeling that there's lots of detail missing because the SatTrack elevation and azimuth are being rounded off to tenths of degrees, which erases any little details like the nearly-diurnal waviness of the SGP4 solution that @rmathar points out. It could be there in the SatTrack data too, but the rounding to tenths of degrees turns whatever curve was there into a stair-step.

As @jmfriedt hasn't replied since April, and since Skyfield is happy to use the official SGP4 solution for satellites, and since there's no volunteers have popped up here in the issue who want to go through SatTrack's code line-by-line figuring out how it works and why its results are different from the official routine, I'm going to go ahead and close this issue. I welcome any further comments on the closed issue and will plan on responding to them when I get the chance, but I don't guess at this point that any Skyfield contributors are going to be jumping in to do a detailed analysis of the differences between SGP4 and SatTrack.

[jmfriedt]

As @jmfriedt hasn't replied since April, and since Skyfield is happy to use the official SGP4 solution for satellites, and since there's no volunteers have popped up here in the issue who want to go through SatTrack's code line-by-line figuring out how it works and why its results are different from the official routine, I'm going to go ahead and close this issue.

This exchange has been most instructive and although I did not participate for lack of technical knowledge I was most interested in learning about such subtleties in geostationary satellite path prediction. On a positive note, I did the same for GPS satellite position prediction and the result is amazingly consistent with online tools and observations. Thanks a lot for looking into the issue, greatly appreciated.