Actual accuracy of Skyfield+SGP4 in computing satellite paths ?

[alxobs]

I recently used Skyfield to compare my telescopic observations of artificial satellite paths with computed ephemeris from TLEs. My goal being to estimate the practical accuracy of the SGP4 inversion.
The recorded images, after plate solving, (12’’ telescope + CMOS camera) allowed satellite paths to be located (in position and orientation) with a precision of a couple of arc seconds in the sky.
My first comparisons, by using standard procedure as described in the Skyfield documentation, were quite unsatisfying, since they lead to discrepancies of order of several arc minutes. Despite the fact that I was using TLE elements with optimal epoch.
I then found that much better results (with an overall position accuracy always better than 20 arc seconds) could be obtained by changing two things:
i) by selecting the WGS72 reference ellipsoid in the SGP4 calculation (the twoline2rv routine taking the WGS84 one by default),
ii) by changing the geographic altitude of my observing site by the orthometric one (i.e. the one above the ellipsoid and not the sea level).
I wonder if other people have done same findings and whether my procedure is correct or not.
alx.

[brandon-rhodes]

First, thank you for opening an issue based on real-world observations — it's a rare treat to get to compare Skyfield's output to what's really seen through a telescope!

Could you share the exact details of how you are calling Skyfield? The default is intended to be WGS72. From what I can see, Skyfield calls it like this:

python-skyfield/skyfield/sgp4lib.py

Line 96 in eff6582

satrec = Satrec.twoline2rv(line1, line2)

When no argument is provided, the sgp4 library is written to use WGS72, unless I did something wrong? Here's how it works in C:

https://github.com/brandon-rhodes/python-sgp4/blob/0f5546db4348641ed245aa97aeefc3df805fcf46/extension/wrapper.cpp#L105

And how it works in the Python fallback code:

https://github.com/brandon-rhodes/python-sgp4/blob/0f5546db4348641ed245aa97aeefc3df805fcf46/sgp4/model.py#L54

But maybe these defaults aren't happening correctly in your case? Let me know how you are calling Skyfield's routines, and we can work to get those defaults actually working properly!

[alxobs]

Thanks for your reply. I agree with most of your remarks.
Following is my very simple Skyfield’s usage for computing small ephemeris :

code.txt

As you can see, the use of WGS72 is forced at the first statement, but this is likely not mandatory.
I also compared observing site definitions by exchanging geoids (72 vs 84), but without significant effect.
I got some better calculations by taking site altitude above the ellipsoid rather than above the sea level (can you confirm this by looking at the sources ?).
An instructive exemple of observations vs. calculations is given by the following figure. The image is a composite of a few camera images taken while tracking ZTF comet a few months ago. During the observation (~3000 images at 1 fps, so in about one hour), five different satellite paths could be recorded at five distant times. The composite image size is 25’x25’, increasing RA on x-axis, Dec on y-axis.
The superimposed red lines are the satellite paths as computed by Skyfield.

The agreement is not too bad, but the distance between the observed and calculated tracks still can reach up to ~45", depending on the satellite. I have no idea whether this result is consistent or not with the expected SGP4 reconstruction accuracy.
Nevertheless, a slighly better fit can be obtained by empirically shifting as a whole the calculated tracks by dRA = 25’’, dDec = -10’’, as shown below.

Let me know if my approach is correct and if you think an improvement is possible.