Ephemeris output from Skyfield

Is it possible to use Skyfield to propagate a TLE forwards and then output its ephemeris, rather than geometric x,y,z at some point in the future?

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:

• propagate a TLE to the new time
• output position.km
• output velocity.km_per_s
• assembling these together with the new epoch to produce a new (SGP4 derived) state vector at that time
• derive orbital elements from the pos-vel

This seems rather a kludge if the capability exists already.

• Position and velocity at some epoch time are a form of ephemeris. In my mind, position/velocity at an epoch time is a better form of ephemeris than is a TLE -- so long as the position and velocity at that epoch time are not derived from a TLE. SGP4 is not that good of a propagator (and that's being nice). The goal was to make something that could propagate orbits using technology from the previous century (previous millennium!), and do so using single precision floating point technology. Mar 1, 2021 at 0:17
• An ephemeris is a table of positions (and optionally velocities) at various times. Orbital elements are not an ephemeris -- they are instructions for how to compute an ephemeris, which is sometimes more useful than the ephemeris itself, and sometimes not. A TLE is instructions for how to configure just one particular piece of software, SGP4, to compute an ephemeris, and that is the only thing for which SGP4 should ever be used: have it turn TLEs into tables of pos & vel, and then process them further with something else, having better physics and less confusing definitions. Mar 1, 2021 at 4:55

2 Answers

Skyfield, like most other tools that claim the ability to process TLE data, depends on a bunch of people's educated but aged guesses about what SGP4 might have been doing some time ago, rather than giving access to what SGP4 actually does right now. This made sense for a long time because there was no other way for most people to get hold of SGP4, but happily that is no longer the case. Anyone who registers for a free account can now download the official US government implementation of SGP4 and associated utilities, some of which were significantly updated in November 2020, from https://www.space-track.org/documentation#/sgp4 .

To see the code I'm talking about below, you have to agree to the SGP4 open source license restrictions (which will pop up when you click the download link, and are later available for your reference as \Sgp4Prop_small\SampleCode\*\SGP4_Open_License.txt), download that file, unzip it, and look at the contents of Sgp4Prop_small\SampleCode\Python\wrappers\ and Sgp4Prop_small\SampleCode\Python\DriverExamples\Sgp4Prop\src\ . Wrapper code is also provided in C, C#, Fortran, Java, Matlab, and Visual Basic, (the actual libraries are provided as DLLs for Windows and .so's for Linux) but you asked about Python so I'll describe it all from the Python point of view.

As part of that update, the provided Python wrapper has finally changed from Python 2 to Python 3 (Sgp4Prop.py as included with library versions up through 7.9, unchanged since January 2013 but still being distributed in October 2020, crashed under Python 3 because it wasn't using parentheses in print statements, among other basic flaws). Personally, I stick in a few tweaks of my own that I think make the libraries easier to use, like adding an empty __init__.py to the directory and prepending dots to convert the library wrappers' calls to each other from absolute to relative imports, but that's entirely optional and not part of the official distribution, which is all I'm going to detail at the moment.

The one distinctly non-Pythonic thing to contend with is that the interface provided by the default wrappers follows the old C and Fortran style exactly, so all the functions only return an integer status code. The data comes back as written into structures you must pre-declare using ctypes and provide as input arguments, which isn't hard but is very much not the way Python is typically used. It's quite possible to hide all of that by writing your own wrappers, which I have done for myself but not yet made available to anyone outside work.

The task of converting a TLE into a table of osculating Keplerian elements in the best available way, using as correctly as possible everything SGP4 models, is one of the basic things shown as an example in the demonstration script, Sgp4Prop_small\SampleCode\Python\DriverExamples\Sgp4Prop\src\Sgp4Prop.py . That script makes the unfortunate choice (line 129 calls TimeFunc.TConLoadFile) of wanting you to supply the start, stop, and step times for the file in the painfully archaic and needlessly confusing "6P card" format (e.g. , as described on pages 7–9 of Sgp4Prop_small\Documents\librarydocuments\TimeFunc.doc), but the TimeFunc library provides all the necessary conversion functions to generate the time format the propagator wants to receive (decimal days since 1950, UTC) from whatever other time data you may have.

...so that I can see what is happening between published TLEs.

This is a great reason to want such a feature; watching the evolution of orbital elements over time can provide some intuitive insight as to what's going on. (examples: 1, 2). It does not immediately tell you where the spacecraft is at any given point nor how fast it's going (that's what state vector customers usually want) but it would be a great option!

I don't know if it's in Skyfield yet or not (but I'm pretty sure not), but the Github site is quite active and lively and there is even a discussions page in addition to the issues page.

As I see it there are several ways to go:

1. Calculate an osculating Keplerian element converter from a single or array of state vectors to a single or array of osculating orbit elements, which is what Horizons has. (but see problems with Keplerian elements below)
2. Calculate an SGP-like approximate mean orbital element (would it contain SGP's drag terms?)
3. Calculate something better or more standard than SGP, perhaps Brouwer-Lyddane mean orbital elements.
4. All of the above!

Space SE certainly has all the equations to convert a state vector to a Keplerian orbit based on Earth's central force $$GM$$, but those are guaranteed to be wrong by order one part per thousand because they don't even include the Earth's oblateness as expressed by $$J2$$ nor any other effect, so I don't think that this is in the running for what you want.

This is easy to see; obtain the osculating elements for some objects in Horizons for say ten times distributed over one period and see how much each orbital element actually varies within a single orbit.

Further reading here:

Further reading offsite: