Output element set from propagated TLE

[samokus]

I think this is a clarification request relating to the API webpages. Is it possible to use Skyfield to propagate a TLE forwards and then output an element set, rather than the geometric x,y,z ? I'd like to be able to produce plots of eccentricity and inclination progression between published TLEs.

I've worked through the examples in the API and may have misread something but it appears that all the available outputs are for eclipse events, rise and set times, viewing angles etc. I'd like to see the six elements (SMA, eccentricity, inclination, RAAN, ArgP, mean anomaly) so that I can see what is happening between published TLEs.

I can see a way that might work. It seems right in principle though I'm a bit wary as I'm new to Skyfield:

1. propagate a TLE to the new time
2. output position.km
3. output velocity.km_per_s
4. assembling these together with the new epoch to produce a new (albeit SGP4 derived) state vector at that time
5. derive orbital elements from the pos-vel (using non-skyfield routines)

Are steps 1..4 within the scope of how Skyfield is capable of being used, I mean for short propagation intervals of a day or two? (and obviously I need to take responsibility myself for step 5, unless there is a way within Skyfield already to do this.

Many thanks for any clues.

NB this is actually an adaptation of a question I originally posted here https://space.stackexchange.com/questions/50441/ephemeris-output-from-skyfield?noredirect=1#comment164336_50441

[brandon-rhodes]

Steps 1–4 are indeed features of the existing API. Step 5 might have one of two different approaches:

1. There is contributed code in elementslib.py and keplerlib.py that can perform the conversion, but I've never tried to use it to model objects in Earth orbit. But it might produce the numbers you need. The API is not quite mature yet, though, which is why it's not yet as heavily documented as other parts of Skyfield.
2. There's an SGP4 orbit element calculator called rv2coe() in the sgp4 Python library that Skyfield uses. I mentioned it over in another Skyfield issue:

https://github.com/skyfielders/python-skyfield/issues/278

If you are able to get either approach working and would like to share the code you come up with, I could use your example as the basis for a new section on the Earth satellite documentation page, to help other folks with a similar need.

[samokus]

Thank you for the prompts. So far rv2coe() seems to run ok. My next step is to use Satrec to create a new TLE from the propagated position and then see how much of a mismatch there is between the original TLE and the newly created one as they are both propagated ahead a few days. I don't actually need to turn the rv2coe() into a TLE for my own purposes but it seems like a convenient way to check whether it is working well enough. I'm not a software developer and am completely new to github though once I have things sorted out I'd be happy to show you what I've done. I'm assuming at this point doing this as a background task so it will be a week or two at least.

[brandon-rhodes]

There's no rush! I'll be happy to look at what you've accomplished whenever you can get around to writing it up; and if you're simply too busy, feel free to close the issue instead and not feel obligated.

[samokus]

I've come up with some routines that do more or less what I want in conjunction with rv2coe(). That routine requires the user to provide mu. I have had a look around web pages for skyfield and sgp4 and eventually found three values for mu for Earth of = 398600.79964, 398600.8 and 398600.5 at the bottom of https://github.com/brandon-rhodes/python-sgp4/blob/master/sgp4/propagation.py

Here is one of the subroutines I wrote, without any propagation step for simplicity, i.e. the elements output should be identical to the TLE. If I use the subroutine below then the values output by rv2coe() aren't close enough to those from the TLE with any of the above values for mu which leads me to think I'm doing something wrong.

def getElementsfromTLE(lineString0,lineString1,lineString2):
    # The linestrings correspond to the name and data lines of a TLE
    #
    ts = load.timescale()
    TLE = EarthSatellite(lineString1, lineString2, lineString0, ts)
    statePV = TLE.at(TLE.epoch)
    elements = rv2coe(statePV.position.km, statePV.velocity.km_per_s, mu)
    return elements

If I use the original TLE of

    line0 = 'Angosat'
    line1 = '1 43087U 17086A   21096.26189822 -.00000168  00000-0  00000-0 0  9996'
    line2 = '2 43087   2.8112  87.4413 0013051 120.2632 256.8210  0.99384049 11878'

then I get reasonable agreement except for the Argument of Perigee and Mean anomaly, i.e. 120.2632 vs 117.8832 and 256.8210 vs 258.7510. However summing them shows their errors cancel out a little, still leaving a net error of ~0.4 deg. All the angle outputs below are converted to degrees.

a                    ecc                incl               omega             argp                 nu             m  
42416.51946162975 0.0013094867564974984 2.9175498889347944 87.63418610196786 117.88322000147812 258.603895254296 258.7510214779837 

I'm sorry about the dreadful formatting, the "insert code" option seems not to respect the default Spyder console output formatting so I used the "quote" option which is marginally better. Let me know if you need to see more code / different formatting.

[brandon-rhodes]

Thanks, I think I'll have a chance to answer in more detail this afternoon — for the moment I've edited your comment to help you with the formatting. Putting a line with triple-backticks on it before and after a piece of code is how to get it formatted like code on GitHub. (Or you can indent it by 4 spaces, but it's annoying to have to add extra indentation to every line when pasting, I find.)

[brandon-rhodes]

Okay!

So, your mention of mu made me pause, because I felt like I had seen mu somewhere in the satellite record where it might already be exposed for users to access, so you would not have to copy it from elsewhere in the code. And sure enough, in the commit https://github.com/brandon-rhodes/python-sgp4/commit/2124370255d2a7d160038626e9fe809a525ab403 exactly one year ago tomorrow (!) another project contributor @interplanetarychris, while adding support for another feature, happened to expose mu where users can now access it. For a given satellite you should be able to:

print(sat.model.mu)

— to learn its value.

To discover that a mu attribute was available, folks either had to have an IDE that would show them, or an automated tool like pydoc sgp4.api.Satrec; the attribute wasn’t mentioned in the hand-written docs that appear on PyPI. So, given how many recent questions have been not about specific positions but about the shape of the orbits themselves, I rolled up my sleeves and yesterday and this morning have managed to put together the entire list of Satrec attributes.

So, you can now find out that mu is available simply by reading the main module docstring, which will become the documentation on PyPI the next time I make a release.

So I'd better arrange for a release this week to get these new docs available.

I also discovered that the Satrec attributes at the C++ level include a set of “averaged” orbital elements that are re-computed each time you ask for a position! I had never understood their significance, but thanks to your question, I now see that they might be available for folks like you who want to understand how an orbit is evolving over time. The commit that added my new documentation to the sgp4 library is here:

https://github.com/brandon-rhodes/python-sgp4/commit/db1a9cd6235eebc5b1da6eac64d92d0c13abb5aa

Here's the except specifically about the averaged elements:

Partway through each propagation, the SGP4 routine saves a set of
“singly averaged mean elements” that describe the orbit’s shape at the
moment of propagation.  They are averaged with respect to the mean
anomaly and include the effects of secular gravity, atmospheric drag,
and — in Deep Space mode — of those pertubations from the Sun and Moon
that SGP4 averages over an entire revolution of each of those bodies.
They omit both the shorter-term and longer-term periodic pertubations
from the Sun and Moon that SGP4 applies right before computing each
position.

| ``am`` — Average semi-major axis (earth radii).
| ``em`` — Average eccentricity.
| ``im`` — Average inclination (radians).
| ``Om`` — Average right ascension of ascending node (radians).
| ``om`` — Average argument of perigee (radians).
| ``mm`` — Average mean anomaly (radians).
| ``nm`` — Average mean motion (radians/minute).

These should all be available on the .model attribute where a Skyfield satellite keeps the raw underlying SGP4 object. Could you take a look and see whether the values are useful for your project?

[samokus]

Interesting development. I will need to take a while to wade through this though I wonder if the solution path I took may not work. Let me show you getPropagationElements(), which is another of the routines I wrote, very similar to the first except that the epoch is determined differently so as to allow user input, which in turn depends upon set_updated_time(), see below.

def set_updated_time(InputTLE,timestep_list,ts):
    # InputTLE is a satellite object, such as one made by the Skyfield method
    # EarthSatellite(line1, line2, line0, ts)
    # timestep_list has six elements in the form year,month,day,hour,min,second
    # where each value is the amount that the epoch of InputTLE is desired to
    # change, e.g. [0,0,0,1,0,0] to increase by one hour
    # ts is a Skyfield TimeScale
    # 
    # The return parameter is the new time object resulting from the changed epoch.
    # The internal workings are kept consistently in .tt rather than .utc
    currentcal = list(InputTLE.epoch.tt_calendar())    #currentcal and newcal are also calender dates
    newcal = [x + y for x, y in zip(currentcal, timestep_list)] #element-wise summation of the two date lists 
    return ts.tt(newcal[0],newcal[1],newcal[2],newcal[3],newcal[4],newcal[5])

def getPropagationElements(lineString0,lineString1,lineString2,t):
    # The linestrings correspond to the name and data lines of a TLE
    # t is the amount of time to add to get to the propagation point in seconds
    ts = load.timescale()
    TLE = EarthSatellite(lineString1, lineString2, lineString0, ts)
    tnew = set_updated_time(TLE,[0,0,0,0,0,t],ts)      #Calc new time step in calender terms
    statePV = TLE.at(tnew)
    elements1 = rv2coe(statePV.position.km,statePV.velocity.km_per_s,mu)
    return elements1

Am I right that I would extract the elements with a statement such as

sat.model.em
# or, using my variable names:
TLE.model.em

for eccentricity? If so then I presume I won't be able to use the line in getPropagationElements() statePV = TLE.at(tnew) as the propagation step because TLE.model.em and all the other elements would still point back to the original TLE and not the propagated state.

I will have a closer look through and see what I can get to work, I'm rather on the edge of my comprehension at the moment. I presume there is still the option of creating a new propagated TLE from satrec.sgp4init()

As an aside, I have already discovered a similar looking but different feature for accessing the range of elements, namely sat.model.ecco for eccentricity etc though I presume this is used in the same way as sat.model.em

[brandon-rhodes]

If so then I presume I won't be able to use the line in getPropagationElements() statePV = TLE.at(tnew) as the propagation step because TLE.model.em and all the other elements would still point back to the original TLE and not the propagated state.

Happily, they are not the original TLE elements. The elements like TLE.model.em refer to the propagated state, not the original state at the epoch. That is what I was trying to communicate above, with the sentence:

Partway through each propagation, the SGP4 routine saves a set of
“singly averaged mean elements” that describe the orbit’s shape at the
moment of propagation. …

Could you suggest how I might make the description clearer, so that it would better answer the question of whether .em is the original element or the propagated one? I haven’t released these new docs yet and would love to make them clearer!

As an aside, I have already discovered a similar looking but different feature for accessing the range of elements, namely sat.model.ecco for eccentricity etc though I presume this is used in the same way as sat.model.em

The .ecco attribute should be the exact value from the TLE, so it’s the original SGP4 eccentricity at the moment of epoch. I’ll go back to those docs I just wrote, and emphasize that .ecco and its friends are constant and don’t change once they’ve been copied into the record from the TLE!

[samokus]

Hi

Thanks so much for looking at this in detail. I'd be very happy to help provide a description so that you can up date documentation. I presume you'd like to get this out soon though I'm not able to respond fully immediately, I am doing some long days lately with outputs on other projects. This is also why I didn't grasp what you had previously unearthed about the propagated elements - I've simply been doing everything by the seat of my pants and didn't read it properly, sorry about that.

Now that I've paused for a moment I think the phrase you already used "The elements like TLE.model.em refer to the propagated state, not the original state at the epoch" is quite clear already.

I do still have a reservation or too on the previous description though. When I look back at the paragraph:

Partway through each propagation, the SGP4 routine saves a set of “singly averaged mean elements” that describe the orbit’s shape at the moment of propagation. They are averaged with respect to the mean anomaly and include the effects of secular gravity, atmospheric drag, and — in Deep Space mode — of those pertubations from the Sun and Moon that SGP4 averages over an entire revolution of each of those bodies. They omit both the shorter-term and longer-term periodic pertubations from the Sun and Moon that SGP4 applies right before computing each position.

then the bits that I find ambiguous are in bold above. Its probably a good point to mention here that my background with propagators is a) using them rather than writing and b) with an applied force model approach rather than whatever SGP4 gets up to. I have seen the lagrange planetary equations turned into a secular form though I don't really feel confident in terms such as "singly averaged" and "averaged with respect to the mean anomaly". That said, let me plough on a bit further with the bold text that I find ambiguous:

1. at the moment of because to me propagation isn't a "moment", it has a start and a finish (i.e. original state and propagated state). To address this we just we need to be confident which of the latter two points in time are referred to.
2. They omit ... that SGP4 applies right before computing each position seems to imply that the propagation occurs in two processes and that the saved elements are only subject to the first set of processes (secular ones) rather than the second set ("shorter-term and longer-term periodic pertubations from the Sun and Moon"). This seems rather unsettling to me, I originally interpreted it as being a propagation subject to some effects of physics and not others. I think the way to address this depends on the user level. If they are like me then they just need reassurance that my interpretation this way is incorrect, though some might want more background details I guess.

Once I've had time to settle and catch up with this project I will look through the documentation links that you've provided and also simply try out these newly found element outputs. if I see anything else in the documentation links that jump out at me I'll let you know.

In the mean time. I'm actually not clear how I apply the .model.em outputs. In the context of the getPropagationElements() method of my previous post, it is statePV that is updated at tnew, not "TLE". I don't see a an equivalent function to .at(tnew) that updates TLE (i.e. the variable) rather than statePV.

Lastly, notwithstanding that I haven't had a chance to check out the new model.mu parameter, what is the status now of the rv2coe() routine? What is its intended use?

[brandon-rhodes]

Thanks for the comments on the documentation! Yesterday I spent a bit of time trying to learn more about the different phases of the propagation calculation myself; like you, I've always just been a user and have never implemented anything like SGP4 myself, but merely call it from Python. Once I understand better the effects of (a) the secular adjustments that precede and are part of creating the average elements, and (b) the periodic adjustments that are applied after the average elements have been computed, I'll see if I can improve the paragraph you pointed out.

(Though: by “at the moment of” I meant “at the exact time t for which you want a satellite position”; not the clock time when you run the program. I’ll have to make that clearer.)

Meanwhile, I've gone ahead and pushed a quick release so the new list of attributes can be viewed right on the Package Index:

https://pypi.org/project/sgp4/#attributes

In the mean time. I'm actually not clear how I apply the .model.em outputs. In the context of the getPropagationElements() method of my previous post, it is statePV that is updated at tnew, not "TLE". I don't see a an equivalent function to .at(tnew) that updates TLE (i.e. the variable) rather than statePV.

Running TLE.at(tnew) not only returns vectors, but also writes a number of new values into the attributes of the TLE.model object itself. So each time you call TLE.at(tnew), you’ll see new average elements have been written to fields like TLE.model.em.

Lastly, notwithstanding that I haven't had a chance to check out the new model.mu parameter, what is the status now of the rv2coe() routine? What is its intended use?

I have no idea. David Vallado included it in SGP4 with the description:

    *  this function finds the classical orbital elements given the geocentric
    *    equatorial position and velocity vectors.

Its little history suggests it's been there for quite a while:

    *  revisions
    *    vallado     - fix special cases                              5 sep 2002
    *    vallado     - delete extra check in inclination code        16 oct 2002
    *    vallado     - add constant file use                         29 jun 2003
    *    vallado     - add mu                                         2 apr 2007

But I am not sure why it's there, unless it's to help folks turn position/velocity vectors into osculating elements that can then be iterated on to create good SGP4 elements.

[samokus]

OK.

The TLE.at(tnew) feature sounds ideal and thanks for the link. I'm still rather in the thick of it but will return to wade through the new documentation link. At a quick glance it looks very helpful.

For rv2coe() I had wondered if it was a generic cartesian to elements transform that was really aimed at keplerian representations of orbits, independent of the type of propagator. I think I have always assumed that in a force model propagator the state-vector is the starting point and the keplerian elements are largely an abstraction though because there is a one to one mapping between any given state vector and an associated set of elements that essentially means we can treat the keplerian elements as if they were the same thing.

I didn't quite follow your last point

can then be iterated on to create good SGP4 elements.

Are posvels generated from a TLE supposed to be compatible in any sense with posvels from osculating orbits or does the TLE averaging process introduce biases? I think I had assumed the latter. Its a bit late here (I'm in the UT +1 timezone, UK) and I'm not quite thinking straight but this seems as if it might be relevant to your last point.

[brandon-rhodes]

I didn't quite follow your last point “can then be iterated on to create good SGP4 elements.” Are posvels generated from a TLE supposed to be compatible in any sense with posvels from osculating orbits or does the TLE averaging process introduce biases? I think I had assumed the latter.

I'm not sure what you mean by “compatible”, but the posvels should certainly be similar. The SGP4 routine essentially takes the average orbit, then does some stretching and adjustment to account for how the Moon's and Sun's gravity warp the orbit away from the perfect ellipse it would be if only the satellite and Earth existed. So posvels generated from the average orbit parameters will be similar to, but not the same as, the posvels generated by SGP4 that take the warped shape of the orbit into account.

[brandon-rhodes]

I'm going to be hopeful and assume that the discussion in this issue went quiet because your needs were indeed met by the intermediate elements available after propagating a satellite. I'm going to go ahead and close this issue, but feel free to ask any further questions even once it's closed.

[dcajacob]

I've just released a proof of concept for this: https://github.com/dcajacob/tle-tailor