Difference between the TLE and the model at the epoch of the TLE file

[matvidal]

When I propagate the orbit using the twoline2rv and propagate functions to the epoch day, I don't obtain the same orbital parameters that appear on the TLE file. Am I doing something wrong or is it a limitation of the model itself?

I made the following example using the current TLE of the ISS:

from numpy import pi, sqrt
from sgp4.earth_gravity import wgs72, wgs84
from sgp4.ext import rv2coe
from sgp4.io import twoline2rv

line1 = "1 25544U 98067A   21041.71818700  .00002037  00000-0  45184-4 0  9998"
line2 = "2 25544  51.6445 248.2549 0002750 357.9058 122.1739 15.48954337268990"
rad2deg = 180/pi
G = 6.67408*10**(-20) # Gravitational constant
M = 5.9722*10**24     # Earth's mass
mu = G*M              # Standard gravitational parameter

#====================== TLE orbital parameters and epoch =====================#
TLE_i = 51.6445       # Inclination
TLE_RAAN = 248.2549   # Right Ascention of the Ascending Node
TLE_e = 0.0002750     # Eccentricity
TLE_w = 357.9058      # Argument of the Perigee
TLE_MA = 122.1739     # Mean Anomaly
TLE_MM = 15.48954337  # Mean Motion (revolutions per day)
year = 2021
month = 2
epoch_day = 41.71818700
day = epoch_day - 31
hour = (day - int(day))*24
minute = (hour - int(hour))*60
second = (minute - int(minute))*60

#=================================  SGP4 model ===============================#
sat_model = twoline2rv(line1, line2, wgs84)#wgs72)
pos, vel = sat_model.propagate(year, month, int(day),
                               int(hour), int(minute), second)
p, a, e, i, RAAN, w, theta, MA, argl, tlon, lonp = rv2coe(pos, vel, mu)

#================== Difference between the TLE and the model =================#
i_diff = abs(TLE_i - i*rad2deg)
RAAN_diff = abs(TLE_RAAN - RAAN*rad2deg)
e_diff = abs(TLE_e - e)
w_diff = abs(TLE_w - w*rad2deg)
MA_diff = abs(TLE_MA - MA*rad2deg)
a = a
n = sqrt(mu/(a**3))
period = 2*pi/n
MM = 86400/period
MM_diff = abs(TLE_MM - MM)
def print_differences():
    print("Incl. (TLE): {}°, incl. (model): {}°, diff.: {}°".format(TLE_i, i*rad2deg, i_diff))
    print("RAAN (TLE): {}°, RAAN (model): {}°, diff.: {}°".format(TLE_RAAN, RAAN*rad2deg, RAAN_diff))
    print("Eccentricity (TLE): {}, Eccentricity (model): {}, diff.: {}".format(TLE_e, e, e_diff))
    print("w (TLE): {}°, w (model): {}°, diff.: {}°".format(TLE_w, w*rad2deg, w_diff))
    print("MA (TLE): {}°, MA (model): {}°, diff.: {}°".format(TLE_MA, MA*rad2deg, MA_diff))
    print("Mean motion (TLE): {}, mean motion (model): {}, diff.: {}".format(TLE_MM, MM, MM_diff))

print_differences()

[brandon-rhodes]

Good question, and, thanks for a complete working example!

First, I would run a few dozen positions over the subsequent few weeks after the epoch, using both the original parameters and the computed ones, to see whether the two element sets are really very distinct in what they predict. Especially when an orbit is very nearly round, with an eccentricity close to zero, it can become almost arbitrary which pair of angles to put in for w and MA, so long as they put the satellite at the right place at the moment of epoch.

Second, you could also see if you could download a few historic ISS orbit sets separated by a few weeks. The orbits should be nearly the same (as long as there's no course correction in between), and you could see whether even for the published elements the w tends to be stable over time, or whether it wildly varies because of how close to circular the orbit is.

[matvidal]

Thanks for your quick response. Actually, the position of the satellite is my real problem. What I'm really trying to do is the following:

1. Using a given TLE file propagate the orbit from time t0 (original epoch) to time t1.
2. At t1 generate a new TLE file using the orbit parameters at time t1 and the epoch in t1.
3. Propagate the orbit again in t1 to see the new orbital parameters and to check if the position is the same or not.

The thing is, with the FLOCK 4P-1 satellite, I have a distance error of around 14 km between the TLE 1 at t1 and the TLE 2 at t1. The positions should be the same because they are from the same moment. The difference between t0 and t1 is less than a day, and I can't find the problem.

I created a new example code but simplified it to eliminate variables such as the difference between t1 and to. In this code, I use a TLE 1 at t0, propagate the orbit to t0, create a TLE 2 at t0, and then propagate again at t0. I used the MOLNIYA 2-9 satellite to show that the eccentricity is not the issue. I still believe that fixing the first example should fix the next one.

from numpy import array, pi, sqrt
from sgp4.earth_gravity import wgs72, wgs84
from sgp4.ext import rv2coe
from sgp4.io import twoline2rv

def checksum(line):
    check = 0
    for char in line[:-1]:
        if char.isdigit():
            check += int(char)
        if char == "-":
            check += 1
    return check % 10

#====================================== Original TLE file ===================================#

line1 = "1 07276U 74026A   21040.94023498  .00000042  00000-0  00000+0 0  9992"
line2 = "2 07276  64.1795 301.0115 6721457 286.2818  13.8575  2.45094233237311"
rad2deg = 180/pi
G = 6.67408*10**(-20) # Gravitational constant
M = 5.9722*10**24     # Earth's mass
mu = G*M              # Standard gravitational parameter
print(line1)
print(line2)

#======================================  First propagation ==================================#
year = 2000 + int(line1[18:20])
month = 2
epoch_day = float(line1[20:32])
day = epoch_day - 31
hour = (day - int(day))*24
minute = (hour - int(hour))*60
second = (minute - int(minute))*60
sat_model = twoline2rv(line1, line2, wgs84)#wgs72)
pos0, vel0 = sat_model.propagate(year, month, int(day),
                                 int(hour), int(minute), second)
print("Position after first propagation from the original TLE file:\n{}".format(pos0))
p, a, e, i, RAAN, w, theta, MA, argl, tlon, lonp = rv2coe(pos0, vel0, mu)

#========== New TLE file generated by the new orbital parameters at the same epoch ==========#

a = a
n = sqrt(mu/(a**3))
period = 2*pi/n
MM = 86400/period
aux="{:.8f}".format(e)
formatted_e = "{}".format(aux[2:-1])
new_line2 = "{} {:8.4f} {:8.4f} {} {:8.4f} {:8.4f} {:11.8f}{}7".format(line2[0:7],
                                                                       i*rad2deg,
                                                                       RAAN*rad2deg,
                                                                       formatted_e,
                                                                       w*rad2deg,
                                                                       MA*rad2deg,
                                                                       MM,
                                                                       line2[63:68])
checksum2 = checksum(new_line2)
new_line2 = "{}{}".format(new_line2[0:-1], checksum2)
print(line1)
print(new_line2)

#====================================  Second propagation ===================================#
new_sat_model = twoline2rv(line1, new_line2, wgs84)#wgs72)
pos1, vel1 = new_sat_model.propagate(year, month, int(day),
                                     int(hour), int(minute), second)
print("Position after second propagation from the second TLE file:\n{}".format(pos1))
p, a, e, i, RAAN, w, theta, MA, argl, tlon, lonp = rv2coe(pos1, vel1, mu)
print("Difference between positions: {}".format(array(pos0) - array(pos1)))
distance = sqrt((pos0[0] - pos1[0])**2 + (pos0[1] - pos1[1])**2 + (pos0[2] - pos1[2])**2)
print("Distance: {} km".format(distance))

[brandon-rhodes]

Thanks, that example is simpler! So the core problem is that rv2coe() → propagate() is not performing a round trip: the orbit that comes back doesn't propagate to the requested vectors.

I know that TLEs are produced by curve fitting (least squares or somesuch) and not from a single vector pair, but I've never studied how the rv2coe() parameters behave if fed back into the SGP4 propagator instead of a standard simpler Kepler routine (which presumably would produce the same vectors exactly?). The routine is only mentioned in the Report —

https://www.celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev1.pdf

— as being present for testing, and indeed I see in tests.py that it's used to write the extended fields in the longer lines in the https://github.com/brandon-rhodes/python-sgp4/blob/master/sgp4/tcppver.out#L3 file, where not only the time and x,y,z,xdot,ydot,zdot are printed, but orbital elements showing how they evolve over the hours since the epoch.

I'm not familiar enough with the SGP4 internals (this Python library merely wraps the official C++ code) to know which internal corrections are responsible for SGP4 not behaving like a simple Kepler propagation.

If my goal were to create orbit elements that reproduced exactly a vector position and velocity, I would try starting with the rv2coe() parameters as a guess, and then use something like https://docs.scipy.org/doc/scipy/reference/optimize.html to narrow down the parameters until the output vectors at the epoch matched exactly. You might try that approach and see what happens. Studying the SGP4 source and the publications surrounding the algorithm might suggest which corrections inside of it prevent it from simply returning raw Kepler positions at a given time.

In the meantime, I'll mark this issue as an open question; maybe someone more expert on SGP4 than I am will see the question and know more details!

[dcajacob]

Brandon's intuition is correct. rv2coe returns the osculating elements, not the SGP4 elements. They're different. He's also correct about the generation of TLEs from state vectors. It's not a simple process. You need to take many state vectors over time and "fit": an estimated TLE that when propagated, minimizes the error between the resultant state vectors and the source state vectors.

On Thu, Feb 11, 2021 at 8:30 AM Brandon Rhodes notifications@github.com wrote:

Thanks, that example is simpler! So the core problem is that rv2coe() → propagate() is not performing a round trip: the orbit that comes back doesn't propagate to the requested vectors.

I know that TLEs are produced by curve fitting (least squares or somesuch) and not from a single vector pair, but I've never studied how the rv2coe() parameters behave if fed back into the SGP4 propagator instead of a standard simpler Kepler routine (which presumably would produce the same vectors exactly?). The routine is only mentioned in the Report —

https://www.celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev1.pdf

— as being present for testing, and indeed I see in tests.py that it's used to write the extended fields in the longer lines in the https://github.com/brandon-rhodes/python-sgp4/blob/master/sgp4/tcppver.out#L3 file, where not only the time and x,y,z,xdot,ydot,zdot are printed, but orbital elements showing how they evolve over the hours since the epoch.

I'm not familiar enough with the SGP4 internals (this Python library merely wraps the official C++ code) to know which internal corrections are responsible for SGP4 not behaving like a simple Kepler propagation.

If my goal were to create orbit elements that reproduced exactly a vector position and velocity, I would try starting with the rv2coe() parameters as a guess, and then use something like https://docs.scipy.org/doc/scipy/reference/optimize.html to narrow down the parameters until the output vectors at the epoch matched exactly. You might try that approach and see what happens. Studying the SGP4 source and the publications surrounding the algorithm might suggest which corrections inside of it prevent it from simply returning raw Kepler positions at a given time.

In the meantime, I'll mark this issue as an open question; maybe someone more expert on SGP4 than I am will see the question and know more details!

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/brandon-rhodes/python-sgp4/issues/86#issuecomment-777457796, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAYP3CRCTXLVLGCJTGL2WC3S6PLVTANCNFSM4XNUZZFA .

-- Very Respectfully,

Dan CaJacob

[brandon-rhodes]

Thanks very much for weighing in with your greater experience of the process, @dcajacob!

I’m going to go ahead and close this issue, @matvidal, since I think we now have an answer to your question. But please feel free to add more comments here if you want to follow up further. I would be interested in someday learning more about the art of producing real-world TLE elements, as more questions on that topic are arriving each year as more groups around the world have satellites they are managing — but, alas, my impression is that the best approaches are all more or less proprietary at this point.

[matvidal]

Thank you both @brandon-rhodes and @dcajacob for your answers. I haven't replied because I have been doing research on this topic. I found this document from JPL that I'm studying right now: https://trs.jpl.nasa.gov/bitstream/handle/2014/44785/10-0441_A1b.pdf?sequence=1&isAllowed=y

If I implement a solution in python I will post it here. Thanks again and take care.

[brandon-rhodes]

@matvidal — Interesting, thanks for that link! I’ll be interested to hear whether you’re able to make progress, and in any case will plan to read that paper someday.