eteq · May 1, 2020 · Apr 15, 2020 · Apr 15, 2020 · Apr 15, 2020 · Apr 15, 2020
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -25,6 +25,8 @@ astropy.coordinates
   This provides limited support for updating coordinate data values from another
   coordinate object of the same class and equivalent frame attributes. [#9857]
 
+- Added the True Equator Mean Equinox (TEME) frame. [#10149]
+
 astropy.cosmology
 ^^^^^^^^^^^^^^^^^
 
@@ -41,7 +43,7 @@ astropy.io.ascii
 astropy.io.misc
 ^^^^^^^^^^^^^^^
 - Added serialization of parameter constraints fixed and bounds.  [#10082]
-- Added 'functional_models.py' and 'physical_models.py' to asdf/tags/transform, 
+- Added 'functional_models.py' and 'physical_models.py' to asdf/tags/transform,
   with to allow serialization of all functional and physical models. [#10028]
 
 - Fix ASDF serialization of circular model inverses, and remove explicit calls

diff --git a/astropy/coordinates/builtin_frames/__init__.py b/astropy/coordinates/builtin_frames/__init__.py
@@ -37,6 +37,7 @@
 from .cirs import CIRS
 from .itrs import ITRS
 from .hcrs import HCRS
+from .teme import TEME
 from .ecliptic import *  # there are a lot of these so we don't list them all explicitly
 from .skyoffset import SkyOffsetFrame
 # need to import transformations so that they get registered in the graph
@@ -55,7 +56,7 @@
 # get included
 __all__ = ['ICRS', 'FK5', 'FK4', 'FK4NoETerms', 'Galactic', 'Galactocentric',
            'galactocentric_frame_defaults',
-           'Supergalactic', 'AltAz', 'GCRS', 'CIRS', 'ITRS', 'HCRS',
+           'Supergalactic', 'AltAz', 'GCRS', 'CIRS', 'ITRS', 'HCRS', 'TEME',
            'PrecessedGeocentric', 'GeocentricMeanEcliptic',
            'BarycentricMeanEcliptic', 'HeliocentricMeanEcliptic',
            'GeocentricTrueEcliptic', 'BarycentricTrueEcliptic',

diff --git a/astropy/coordinates/builtin_frames/intermediate_rotation_transforms.py b/astropy/coordinates/builtin_frames/intermediate_rotation_transforms.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """
-Contains the transformation functions for getting to/from ITRS, GCRS, and CIRS.
+Contains the transformation functions for getting to/from ITRS, TEME, GCRS, and CIRS.
 These are distinct from the ICRS and AltAz functions because they are just
 rotations without aberration corrections or offsets.
 """
@@ -16,11 +16,31 @@
 from .gcrs import GCRS, PrecessedGeocentric
 from .cirs import CIRS
 from .itrs import ITRS
+from .teme import TEME
 from .utils import get_polar_motion, get_jd12
 
 # # first define helper functions
 
 
+def teme_to_itrs_mat(time):
+    # Sidereal time, rotates from ITRS to mean equinox
+    # Use 1982 model for consistency with Vallado et al (2006)
+    # http://www.celestrak.com/publications/aiaa/2006-6753/AIAA-2006-6753.pdf
+    gst = erfa.gmst82(*get_jd12(time, 'ut1'))
+
+    # Polar Motion
+    # Do not include TIO locator s' because it is not used in Vallado 2006
+    xp, yp = get_polar_motion(time)
+    pmmat = erfa.pom00(xp, yp, 0)
+
+    # rotation matrix
+    # c2tcio expects a GCRS->CIRS matrix as it's first argument.
+    # Here, we just set that to an I-matrix, because we're already
+    # in TEME and the difference between TEME and CIRS is just the
+    # rotation by the sidereal time rather than the Earth Rotation Angle
+    return erfa.c2tcio(np.eye(3), gst, pmmat)
+
+
 def gcrs_to_cirs_mat(time):
     # celestial-to-intermediate matrix
     return erfa.c2i06a(*get_jd12(time, 'tt'))
@@ -129,3 +149,26 @@ def precessedgeo_to_gcrs(from_coo, to_frame):
 
     # then move to the GCRS that's actually desired
     return gcrs_coo.transform_to(to_frame)
+
+
+@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, TEME, ITRS)
+def teme_to_itrs(teme_coo, itrs_frame):
+    # first get us to TEME at the target obstime
+    # TODO: self transform?
+    teme_coo2 = teme_coo.transform_to(TEME(obstime=itrs_frame.obstime))
+
+    # now get the pmatrix
+    pmat = teme_to_itrs_mat(itrs_frame.obstime)
+    crepr = teme_coo2.cartesian.transform(pmat)
+    return itrs_frame.realize_frame(crepr)
+
+
+@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ITRS, TEME)
+def itrs_to_teme(itrs_coo, teme_frame):
+    # compute the pmatrix, and then multiply by its transpose
+    pmat = teme_to_itrs_mat(itrs_coo.obstime)
+    newrepr = itrs_coo.cartesian.transform(matrix_transpose(pmat))
+    teme = TEME(newrepr, obstime=itrs_coo.obstime)
+
+    # now do any needed offsets (no-op if same obstime)
+    return teme.transform_to(teme_frame)
diff --git a/astropy/coordinates/builtin_frames/teme.py b/astropy/coordinates/builtin_frames/teme.py
@@ -0,0 +1,52 @@
+# -*- coding: utf-8 -*-
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+A coordinate or frame in the TEME (True Equator, Mean Equinox) system.
+
+TEME is a True equator, Mean Equinox coordinate frame used in NORAD TLE
+satellite files.
+"""
+# -*- coding: utf-8 -*-
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+from astropy.utils.decorators import format_doc
+from astropy.coordinates.representation import CartesianRepresentation, CartesianDifferential
+from astropy.coordinates.baseframe import BaseCoordinateFrame, base_doc
+from astropy.coordinates.attributes import TimeAttribute
+from .utils import DEFAULT_OBSTIME
+
+__all__ = ['TEME']
+
+doc_footer_teme = """
+    Other parameters
+    ----------------
+    obstime : `~astropy.time.Time`
+        The time at which the frame is defined.  Used for determining the
+        position of the Earth.
+"""
+
+
+@format_doc(base_doc, components="", footer=doc_footer_teme)
+class TEME(BaseCoordinateFrame):
+    """
+    A coordinate or frame in the True Equator Mean Equinox frame (TEME).
+
+    This frame is a geocentric system similar to CIRS or geocentric apparent place,
+    except that the mean sidereal time is used to rotate from TIRS. TEME coordinates
+    are most often used in combination with orbital data for satellites in the
+    two-line-ephemeris format.
+
+    Different implementations of the TEME frame exist. For clarity, this frame follows the
+    conventions and relations to other frames that are set out in Vallado et al (2006).
+
+    For more background on TEME, see the references provided in the
+    :ref:`astropy-coordinates-seealso` section of the documentation.
+    """
+
+    default_representation = CartesianRepresentation
+    default_differential = CartesianDifferential
+
+    obstime = TimeAttribute()
+
+# Transformation functions for getting to/from TEME and ITRS are in
+# intermediate rotation transforms.py
diff --git a/astropy/coordinates/tests/test_intermediate_transformations.py b/astropy/coordinates/tests/test_intermediate_transformations.py
@@ -12,8 +12,8 @@
 from astropy.time import Time
 from astropy.coordinates import (EarthLocation, get_sun, ICRS, GCRS, CIRS, ITRS, AltAz,
                 PrecessedGeocentric, CartesianRepresentation, SkyCoord,
-                SphericalRepresentation, UnitSphericalRepresentation,
-                HCRS, HeliocentricMeanEcliptic)
+                CartesianDifferential, SphericalRepresentation, UnitSphericalRepresentation,
+                HCRS, HeliocentricMeanEcliptic, TEME)
 from astropy.utils import iers
 
 
@@ -492,6 +492,46 @@ def test_gcrs_self_transform_closeby():
     assert_allclose(delta, 0.0*u.m, atol=1*u.m)
 
 
+@pytest.mark.remote_data
+def test_teme_itrf():
+    """
+    Test case transform from TEME to ITRF.
+
+    Test case derives from example on appendix C of Vallado, Crawford, Hujsak & Kelso (2006).
+    See https://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf
+    """
+    v_itrf = CartesianDifferential(-3.225636520, -2.872451450, 5.531924446,
+                                   unit=u.km/u.s)
+    p_itrf = CartesianRepresentation(-1033.479383, 7901.2952740, 6380.35659580,
+                                     unit=u.km, differentials={'s': v_itrf})
+    t = Time("2004-04-06T07:51:28.386")
+
+    teme = ITRS(p_itrf, obstime=t).transform_to(TEME(obstime=t))
+    v_teme = CartesianDifferential(-4.746131487, 0.785818041, 5.531931288,
+                                   unit=u.km/u.s)
+    p_teme = CartesianRepresentation(5094.18016210, 6127.64465050, 6380.34453270,
+                                     unit=u.km, differentials={'s': v_teme})
+
+    assert_allclose(teme.cartesian.without_differentials().xyz,
+                    p_teme.without_differentials().xyz, atol=30*u.cm)
+
+    assert_allclose(teme.cartesian.differentials['s'].d_xyz,
+                    p_teme.differentials['s'].d_xyz, atol=1.0*u.cm/u.s)
+
+    # test round trip
+    itrf = teme.transform_to(ITRS(obstime=t))
+    assert_allclose(
+        itrf.cartesian.without_differentials().xyz,
+        p_itrf.without_differentials().xyz,
+        atol=100*u.cm
+    )
+    assert_allclose(
+        itrf.cartesian.differentials['s'].d_xyz,
+        p_itrf.differentials['s'].d_xyz,
+        atol=1*u.cm/u.s
+    )
+
+
 @pytest.mark.remote_data
 def test_earth_orientation_table():
     """Check that we can set the IERS table used as Earth Reference.

diff --git a/docs/coordinates/index.rst b/docs/coordinates/index.rst
@@ -448,6 +448,7 @@ listed below.
    skycoord
    transforming
    solarsystem
+   satellites
    formatting
    matchsep
    representations
@@ -496,6 +497,9 @@ coordinate systems implemented here include:
 * Meeus, J. "Astronomical Algorithms"
     A valuable text describing details of a wide range of coordinate-related
     problems and concepts.
+* `Revisiting Spacetrack Report #3 <https://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf>`_
+    A discussion of the simplified general pertubation (SGP) for satellite orbits, with a description of
+    the True Equator Mean Equinox (TEME) coordinate frame.
 
 
 Built-in Frame Classes

diff --git a/docs/coordinates/satellites.rst b/docs/coordinates/satellites.rst
@@ -0,0 +1,106 @@
+.. include:: references.txt
+
+.. _astropy-coordinates-satellites:
+
+Working with Earth Satellites Using Astropy Coordinates
+*******************************************************
+
+Satellite data is normally provided in the Two-Line Element (TLE) format
+(see `here <https://www.celestrak.com/NORAD/documentation/tle-fmt.php>`_
+for a definition). These datasets are designed to be used in combination
+with a theory for orbital propagation model to predict the positions
+of satellites.
+
+The history of such models is discussed in detail in
+`Vallado et al (2006) <https://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf>`_
+who also provide a reference implementation of the SGP4 orbital propagation
+code, designed to be compatible with the TLE sets provided by the United
+States Department of Defense, which are available from a source like
+`Celestrak <http://celestrak.com/>`_.
+
+The output coordinate frame of the SGP4 model is the True Equator, Mean Equinox
+frame (TEME), which is one of the frames built-in to `astropy.coordinates`.
+TEME is an Earth-centered inertial frame (i.e., it does not rotate with respect
+to the stars). Several definitions exist; ``astropy`` uses the implementation described
+in `Vallado et al (2006) <https://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf>`_.
+
+Finding TEME Coordinates from TLE Data
+==========================================
+
+There is currently no support in `astropy.coordinates` for computing satellite orbits
+from TLE orbital element sets. Full support for handling TLE files is available in
+the `Skyfield <http://rhodesmill.org/skyfield/>`_ library, but some advice for dealing
+with satellite data in ``astropy`` is below.
+
+.. EXAMPLE START Using sgp4 to get a TEME coordinate
+
+You will need some external library to compute the position and velocity of the satellite from the
+TLE orbital elements. The `SGP4 <https://pypi.org/project/sgp4/>`_ library can do this. An example
+of using this library to find the  `~astropy.coordinates.TEME` coordinates of a satellite is::
+
+.. doctest-requires:: sgp4
+
+        >>> from sgp4.api import Satrec
+        >>> from sgp4.api import SGP4_ERRORS
+        >>> s = '1 25544U 98067A   19343.69339541  .00001764  00000-0  38792-4 0  9991'
+        >>> t = '2 25544  51.6439 211.2001 0007417  17.6667  85.6398 15.50103472202482'
+        >>> satellite = Satrec.twoline2rv(s, t)
+
+The ``satellite`` object has a method, ``satellite.sgp4``, that will try to compute the TEME position
+and velocity at a given time::
+
+.. doctest-requires:: sgp4
+
+        >>> from astropy.time import Time
+        >>> t = Time(2458827.362605, format='jd')
+        >>> error_code, teme_p, teme_v = satellite.sgp4(t.jd1, t.jd2)  # in km and km/s
+        >>> if error_code != 0:
+        ...     raise RuntimeError(SGP4_ERRORS[error_code])
+
+Now that we have the position and velocity in kilometers and kilometers per second, we can create a
+position in the `~astropy.coordinates.TEME` reference frame::
+
+.. doctest-requires:: sgp4
+
+        >>> from astropy.coordinates import TEME, CartesianDifferential, CartesianRepresentation
+        >>> from astropy import units as u
+        >>> teme_p = CartesianRepresentation(teme_p*u.km)
+        >>> teme_v = CartesianDifferential(teme_v*u.km/u.s)
+        >>> teme = TEME(teme_p.with_differentials(teme_v), obstime=t)
+
+.. EXAMPLE END
+
+Note how we are careful to set the observed time of the `~astropy.coordinates.TEME` frame to
+the time at which we calculated satellite position.
+
+Transforming TEME to Other Coordinate Systems
+==================================================
+
+Once you have satellite positions in `~astropy.coordinates.TEME` coordinates they can be transformed
+into any `astropy.coordinates` frame.
+
+For example, to find the overhead latitude, longitude, and height of the satellite::
+
+.. EXAMPLE START Transforming TEME
+
+.. doctest-requires:: sgp4
+
+        >>> from astropy.coordinates import ITRS
+        >>> itrs = teme.transform_to(ITRS(obstime=t))  # doctest: +REMOTE_DATA
+        >>> location = itrs.earth_location  # doctest: +REMOTE_DATA
+        >>> location.geodetic  # doctest: +REMOTE_DATA +FLOAT_CMP
+        GeodeticLocation(lon=<Longitude 160.34199789 deg>, lat=<Latitude -24.6609379 deg>, height=<Quantity 420.17927591 km>)
+
+Or, if you want to find the altitude and azimuth of the satellite from a particular location::
+
+.. doctest-requires:: sgp4
+
+        >>> from astropy.coordinates import EarthLocation, AltAz
+        >>> siding_spring = EarthLocation.of_site('aao')  # doctest: +REMOTE_DATA
+        >>> aa = teme.transform_to(AltAz(obstime=t, location=siding_spring))  # doctest: +REMOTE_DATA
+        >>> aa.alt  # doctest: +REMOTE_DATA +FLOAT_CMP
+        <Latitude 10.94798427 deg>
+        >>> aa.az  # doctest: +REMOTE_DATA +FLOAT_CMP
+        <Longitude 59.28807348 deg>
+
+.. EXAMPLE END
diff --git a/docs/install.rst b/docs/install.rst
@@ -103,6 +103,8 @@ The following packages can optionally be used when testing:
 - `skyfield <https://rhodesmill.org/skyfield/>`_: Used for testing Solar System
   coordinates.
 
+- `spgp4 <https://pypi.org/project/sgp4/>`_: Used for testing satellite positions.
+
 - `tox <https://tox.readthedocs.io/en/latest/>`_: Used to automate testing
   and documentation builds.
 

diff --git a/docs/whatsnew/4.1.rst b/docs/whatsnew/4.1.rst
@@ -19,9 +19,9 @@ In particular, this release includes:
 
 2. Support for units on otherwise unitless models via the ``Model.coerce_units`` method.
 
-3. Ability to create a new ``CompoundModel`` by modifying an existing one
-   with the ``replace_submodel`` method.
+3. Added the True Equator Mean Equinox frame.
 
+4. Added the True Equator Mean Equinox frame.
 
 Support for writing Dask arrays to FITS files
 =============================================
@@ -38,3 +38,11 @@ written, which will avoid using excessive memory:
     >>> from astropy.io import fits
     >>> hdu = fits.PrimaryHDU(data=array)
     >>> hdu.writeto('test_dask.fits')
+
+Added True Equator Mean Equinox (TEME) frame for satellite two-line ephemeris data
+==================================================================================
+
+The True Equator Mean Equinox (`~astropy.coordinates.TEME`) frame has been added to
+the built-in frames within :ref:`astropy.coordinates <astropy-coordinates>`.
+
+For more details, see :ref:`astropy-coordinates-satellites`.
diff --git a/setup.cfg b/setup.cfg
@@ -52,6 +52,7 @@ test =
     ipython
     coverage
     skyfield>=1.20
+    sgp4
 all =
     scipy>=0.18
     dask[array]