Discrepancy between Skyfield and manually calculated Orbital Parameters

[karatemir]

Hi there,

I am calculating semi-major axis, perigee and apogee altitude values for certain LEO satellites and came across a discrepancy between my theoretically calculated values and Skyfield provided values. To reproduce I am providing the following script:

# https://orbit.ing-now.com/satellite/39444/2013-066ae/funcube-1/
    # Orbit : #704394.00 @ 577 x 658 km x 97.6 deg
    # Spacecraft Age : 8 years, 10 months
    # Distance Travelled : ≈2.1 billion km
from skyfield.api import EarthSatellite, load, utc, wgs84
import numpy as np

mu = 3.986004418 * 1e14 # gravitational parameter of Earth (m^3s^−2)

fun = """\
1 39444U 13066AE  22257.61551127  .00002978  00000+0  35585-3 0  9990
2 39444  97.6337 227.5603 0057973  99.5629 261.2143 14.83718775474743
"""

# satrec = Satrec.twoline2rv(*fun.splitlines(), WGS84)
# sat    = EarthSatellite.from_satrec(satrec, ts)
sat = EarthSatellite(*fun.splitlines())
r_e = sat.model.radiusearthkm
print(f'Skyfield mu value: {sat.model.mu} vs theory: {mu/1e9}')

print(f'Skyfield results -- Axis: {sat.model.a * r_e}, Perigee: {sat.model.altp * r_e}, Apogee: {sat.model.alta * r_e}')

# theoretical computation

T   = 2 * np.pi / sat.model.no_kozai * 60
semi_major = (T**2 * mu / (4 * np.pi**2))**(1/3) / 1e3
print(f'''
        Theoretical results --
        Axis: {semi_major},
        Perigee: {semi_major*(1+sat.model.ecco) - sat.model.radiusearthkm},
        Apogee: {semi_major*(1-sat.model.ecco) - sat.model.radiusearthkm}''')

The results of the above script are:

Skyfield mu value: 398600.8 vs theory: 398600.4418
Skyfield results -- Axis: 6992.769095396103, Perigee: 574.0949151193616,  
Apogee: 655.1732756728427

        Theoretical results --
        Axis: 6995.747293563754,
        Perigee: 658.1687393487309,
        Apogee: 577.0558477787772

Please note that the values (577 x 658 km) I have manually computed seems to be in agreement with the values in the website.

I have tried changing the ellipsoid from WGS72 to WGS84 using from_satrec but that didn't really change much. There is also some difference in the gravitational constant but that doesn't seem to really explain the discrepancy either. I have also tested a lower inclination starlink satellite, and in that case the difference seems negligible.

Can you please advise? Best regards

[brandon-rhodes]

Three thoughts:

1. You might be interested in tweaking your script to use the pure sgp4 library that Skyfield uses, so that you can rule out whether Skyfield makes any difference here—your script would then be talking to the official SGP4 engine itself from Vallado et al.
2. What are those callto: strings inside of fun?
3. My guess is that you are seeing the difference between Kepler and SGP4, which are different ways of modeling an orbit, though they do have some things in common of course. To answer your question in detail my guess would be that you'll want to work through the math in this routine, to catch the differences between it and Kepler's simpler theory:

https://github.com/brandon-rhodes/python-sgp4/blob/fca7e90de2467308e7cb62ac2d019f383011dbc8/extension/SGP4.cpp#L1746

Let me know if those thoughts make sense!

[karatemir]

@brandon-rhodes sorry replaced the fun. I copied it from my mail (on a different machine now) and I guess outlook just added those.

I will go through the c++ routine to see the differences as you recommended. I believe the Kepler product description of an orbit in terms of Keplerian apogee x perigee altitude is commonly used, so might be useful to have it implemented in EarthSatellite class.

[karatemir]

@brandon-rhodes The reason for the discrepancy is: My formula when simplified corresponds to (sat.no_kozai * sat.tumin) ** (-2.0/3.0) whereas sgp4 and therefore skyfield returns satrec.a = pow( satrec.no_unkozai*satrec.tumin , (-2.0/3.0) ) as here.

no_unkozai is defined here. I don't know enough orbital mechanics to comment on this matter but from what I see most websites use the computation I proposed (sat.no_kozai * sat.tumin) ** (-2.0/3.0) to compute the semi-major axis and the related apogee/perigee values. The numbers don't match when using skyfield/sgp4 returned sat.a, sat.alta, sat.altp values.

Might be worth looking into.

[brandon-rhodes]

Thank you for tracking that down!

Alas, I'm not free to alter that line of code—that's an official line of C++ code from the SGP4 implementation:

https://github.com/brandon-rhodes/python-sgp4/blob/fca7e90de2467308e7cb62ac2d019f383011dbc8/extension/SGP4.cpp#L1476

(Why does it appear twice, you might ask? The code appears in both propagation.py (which I do control) and in the C++ code because propagation.py is a slow Python implementation of SGP4 that's used as a fallback on platforms that can't compile the C++.)

Here's the source of the official code that the sgp4 module bundles:

https://celestrak.org/publications/AIAA/2006-6753/

If you have found libraries that use alternative math for the semi-major axis, then it may be those libraries that need to be updated—they might not have updated their code to match the latest official code from Vallado, Crawford, Hujsak, and Kelso?

In any case, I can't, alas, alter that code without breaking the promise that I made in the docs that the sgp4 library uses a verbatim copy of the official SGP4 algorithm. You maybe could contact one of the authors to ask about why they use the formula they do to compute the semi-major axis? If you do, feel free to comment here on what you learn, since I'm also curious about their choice of formula.

[karatemir]

@brandon-rhodes

funnily enough celestrak itself seems to use my computation in their website:

To verify:

from skyfield.api import EarthSatellite
tle = """\
1 39444U 13066AE  22259.90845277  .00002789  00000+0  33351-3 0  9995
2 39444  97.6341 229.7533 0058030  92.2597 268.5264 14.83731891475089"""
sat = EarthSatellite(*tle.splitlines())

semi_major = (sat.model.no_kozai * sat.model.tumin) ** (-2.0/3.0)
alta = semi_major * (1+sat.model.ecco) - 1.0
altp = semi_major * (1-sat.model.ecco) - 1.0
print(f'Theory: {semi_major * sat.model.radiusearthkm}, {altp * sat.model.radiusearthkm}, {alta * sat.model.radiusearthkm}')
print(f'SGP4: {sat.model.a * sat.model.radiusearthkm}, {sat.model.altp * sat.model.radiusearthkm}, {sat.model.alta * sat.model.radiusearthkm}')

Which prints:

Theory: 6995.708161379013, 576.9770669185305, 658.169255839496
SGP4: 6992.727867151441, 574.014067338361, 655.1716669645218

Definitely looks weird. Someone who knows SGP4 people should probably ask them about this at some point. Anyhow it was fun looking at it :)

[brandon-rhodes]

It's possible that the SGP4-computed properties like alta are for some special purpose that we don't know about, and that a straight re-computation of the apogee and perigee per Kepler's formulae are better for display to users? Let me know if you run into further information clarifying things!

I suppose if one of us had time, we could run a loop with that satellite to find out its actual in-practice apogee and perigee, and see if those numbers are closer to your Kepler-computed one, or the one that SGP4 computes. But it could be that it varies from one satellite to another, whether their real apogee/perigee are close to the SGP4 or Kepler numbers, so one satellite's result wouldn't necessarily generalize to all.

Again, good job tracking down the exact difference in formulae!