Allow change of gravity model in Satrec

[astrojuanlu]

The new Satrec class uses the WGS-72 gravity model:

https://github.com/brandon-rhodes/python-sgp4/blob/10e8076f18bb148d5e2a2434ae20daf9dcfc6fd0/extension/wrapper.cpp#L113

https://github.com/brandon-rhodes/python-sgp4/blob/10e8076f18bb148d5e2a2434ae20daf9dcfc6fd0/sgp4/model.py#L47

There is some discussion in Vallado et al. "Revisiting Spacetrack Report #3: Rev 1" but I haven't fully understood it yet. Do you have any supporting documentation that says that WGS-72 should always be used, or do you think it could be left as a parameter?

[brandon-rhodes]

In making that simplification I relied on VI.B. in the paper:

However, because many operational sites may still have embedded software containing a version of SGP4 using WGS-72, and the fact that the accuracy of the theory would not really be impacted, AFSPC may well have chosen to retain the older set of constants to better maintain interoperability with its internal resources. We use WGS-72 as the default value.

The simplification provided by omitting the parameter was part of making the C wrapper code simple enough that I was able to sit down and get it written and released to users. I figured that folks needing the non-default value could still use the old API, which continues to work with the pure Python code.

Do you have a use case where the difference between them is crucial? I would be happy to look at the difference in output.

[astrojuanlu]

Thanks for the super quick response!

I would be happy to look at the difference in output.

I knew you would ask this :) I'll try to produce some numbers to see the actual difference.

[astrojuanlu]

Here is the difference in output between old API with WGS-84 and 72 and the new one (72):

$ cat test_sgp4.py 
import datetime as dt

from sgp4.api import Satrec, jday
from sgp4.io import twoline2rv
from sgp4.earth_gravity import wgs72, wgs84

TLE_LINES = (
    "1 43204U 18015K   18339.11168986  .00000941  00000-0  42148-4 0  9999",
    "2 43204  97.3719 104.7825 0016180 271.1347 174.4597 15.23621941 46156",
)

REF_DATE = dt.datetime(2019, 1, 1)

print(
    twoline2rv(TLE_LINES[0], TLE_LINES[1], wgs84).propagate(
        REF_DATE.year,
        REF_DATE.month,
        REF_DATE.day,
        REF_DATE.hour,
        REF_DATE.minute,
        REF_DATE.second + REF_DATE.microsecond * 1e-6,
    )
)
print(
    twoline2rv(TLE_LINES[0], TLE_LINES[1], wgs72).propagate(
        REF_DATE.year,
        REF_DATE.month,
        REF_DATE.day,
        REF_DATE.hour,
        REF_DATE.minute,
        REF_DATE.second + REF_DATE.microsecond * 1e-6,
    )
)
print(
    Satrec.twoline2rv(TLE_LINES[0], TLE_LINES[1]).sgp4(*jday(
        REF_DATE.year,
        REF_DATE.month,
        REF_DATE.day,
        REF_DATE.hour,
        REF_DATE.minute,
        REF_DATE.second + REF_DATE.microsecond * 1e-6,
    ))[1:]
)
(op38) juanlu@ephyra:~/Development/delivery-platform/mission/orbit-predictor$ pip list 2>/dev/null | grep sgp4 && python test_sgp4.py 
sgp4              2.3        
((1889.1426827871676, -3284.8744730940857, -5735.840931181199), (-4.611075211657226, 4.463677104958612, -4.08939018343622))
((1888.5188170546378, -3284.458541452594, -5736.306447887009), (-4.61115561156186, 4.464274162733721, -4.088615057046763))
((1888.5188165794123, -3284.458540992506, -5736.306448308384), (-4.611155611799821, 4.464274163147577, -4.088615056321933))

The difference is in the order of ~1 km in position and ~1 m/s in velocity. I find it interesting that there's a small numerical difference between the new API and the old one with the same gravity model, in the order of ~1mm.

[brandon-rhodes]

The difference is in the order of ~1 km in position and ~1 m/s in velocity.

Thanks for checking! Do we know which gravity model the various authorities use when integrating satellite positions? I am not inclined to complicate the new API at this point unless we know of a source that is definitely using one of the other gravity models.

I find it interesting that there's a small numerical difference between the new API and the old one with the same gravity model, in the order of ~1mm.

If you are curious about whether the difference is a tweak to Vallado's source code itself, you could go back and forth through its history looking at the tweaks that have gradually been made:

https://github.com/rirze/sgp4-cpp/tree/master

[astrojuanlu]

Do we know which gravity model the various authorities use when integrating satellite positions?

That's a very good question. After a bit of Googling I was not able to find what https://celestrak.com/ or https://www.space-track.org/ use.

Orekit uses WGS-72 for TLEs and only for TLEs, see the WGS-72 based constants and usage in the library.

If you are curious about whether the difference is a tweak to Vallado's source code itself, you could go back and forth through its history looking at the tweaks that have gradually been made:

Unfortunately, even when "no changes are reported", such changes tend to be massive. But you're right it's the correct place to look at.

I'll try to reach out to other people in the industry who might know better.

[zouhairm]

I know this is not a decision by consensus case, but as a 3rd party user of the API, I vote for allowing the option to change the gravity model.

I don't think it's a major complication to the API to make the gravity-constant a keyword argument to default to wgs72 to be consistent with the choice from the paper, and it allows people who have been using wgs84 to not worry about the implications of upgrading to the new API in order to take advantage of the speed-up and array functionality.

[zouhairm]

@brandon-rhodes - I saw that you assigned this. Were you intending on pushing this feature soon? Otherwise, I'm happy to take a crack at it and submit a PR in the next couple of days as this is currently blocking on one of my projects.

[brandon-rhodes]

No, not soon; I have too many other tasks. Assigning it to myself is merely an attempt to make sure it happens "someday" rather than "never".

This evening I'll do a little basic design work, though, to decide which way I want the approach to go — whether by making the init method more complicated, or giving folks a way of overwriting those values with whichever set of values they want. I'll comment back here, and at that point remove the assignment from myself.

[brandon-rhodes]

All right! I have re-read the source code.

1. This does look like it needs to be a value provided right at init time, rather than a setting that is overridden later through a "reset gravity constants" method, because both init methods (from a 2-line TLE, and from sgp4init itself) call a sgp4() run immediately as part of setting the struct up, which computes all kinds of derivative numbers that are affected by the gravity constant. Setting the gravity constant as a separate step would require a second sgp4() run, making initialization maybe about twice as expensive for alternative-gravity folks as for folks who use the standard settings.
2. So twoline2rv() will, alas, require a new argument, making it slower for everyone who uses the library, whether they need to adjust the gravity constant or not. The alternative would be a separate constructor that took the extra argument, but I suspect the bit of expense of an optional argument is better in the long run than having to maintain two separate versions of twoline2rv().
3. But to speed things up a bit, what if we had the extra argument be an integer rather than a string, so that a full string compare isn't necessary over and over again every time the constructor is called? Instead the integer could be the whichconst enumeration argument itself, and we could expose the various values whichconst can take as all-caps constants at the module level. That should make for a very quick check of the argument and hopefully keep the library close to lightning-fast.

@zouhairm — does that approach sound reasonable? I'll assign the issue to you now so that you can try implementing. Good luck!

[brandon-rhodes]

I have just released 2.5 that includes this improvement. Enjoy!