skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

Add True Equator, Mean Equinox (TEME) to framelib #931

Closed MichaelBonnet closed 9 months ago

MichaelBonnet commented 10 months ago

Right now, we have what seems to be TETE in framelib.py. It would be nice to have TEME as a reference frame as well, since a lot of (at least preliminary) positional astrodynamics happens in TEME as SGP4 is hilariously fast. I could theoretically do this with astropy, but astropy can be as bad as 70% slower than Skyfield.

For example, I have this function:

from datetime import datetime, timezone
from contextlib import closing
from skyfield import framelib
from skyfield.api import load

def get_earth_sun_vector_teme_at_epoch(epoch: datetime) -> np.ndarray:
    """Get the Earth-Sun position vector in TEME at given epoch."""

    ts = load.timescale()
    epoch = epoch.replace(tzinfo=timezone.utc)
    t = ts.from_datetime(epoch)
    bodies = load("de440s.bsp")
    with closing(bodies):
        sun = bodies["sun"]
        earth = bodies["earth"]
        apparent = earth.at(t).observe(sun).apparent()
        vector = apparent.frame_xyz(framelib.true_equator_and_equinox_of_date).km
        return np.array(vector)

but since framelib.true_equator_and_equinox_of_date seems to be TETE instead of TEME, I get some accuracy issues when checking against industry tool-generated validation data.

For what it's worth, when attempting to get a GCRS Earth-Sun vector, I do essentially the same as above without the frame transformation, which is somewhat more accurate:

earth_sun_vector_gcrs_at_epoch = earth.at(t).observe(sun).apparent().position.km

Having a Skyfield-native TEME frame transformation would be great.

brandon-rhodes commented 10 months ago

Happily, the TEME reference frame can be imported directly from sgp4lib— try:

from skyfield.sgp4lib import TEME

It is not in framelib because then anyone importing framelib would also have to incur the expense of importing sgp4lib, where the oddball definition of theta_GMST1982() lives, which is only ever needed with SGP4 because only it seems to use that reference frame.

The TEME class is mentioned in the API docs, but it's probably hard to find there. If I add it to the main Coordinates page, and mention "TEME" in the new section's title, would that have made it easier to find?

https://rhodesmill.org/skyfield/coordinates.html

MichaelBonnet commented 10 months ago

That would definitely make it easier to find. An example showing how using it for frame transformations does (or doesn't) work like with other reference frames using frame_xyz would also be helpful.

MichaelBonnet commented 10 months ago

Oh, and adding something in the docstring clearly stating that the TETE frame is not TEME would probably be good, since a developer new to astrodynamics may not catch the difference in naming immediately.

brandon-rhodes commented 9 months ago

There ­— I have committed an update to the docs that clarifies (hopefully!) both reference frames, and that links from the TETE docstring so that folks also find that TEME is there. The API docs now also show TEME in the list of reference frames. Hopefully that prevents future confusion; let me know if you see any typos.

I'll plan to push this update to the online documentation within the next few days, at which point the change (which currently you can see in diff form by clicking on the commit right above this comment) should appear on the website.