eteq · Sep 29, 2015 · Sep 28, 2015 · May 2, 2015 · May 5, 2015 · Sep 28, 2015
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -39,6 +39,9 @@ New Features
   - Added ``Supergalactic`` frame to support de Vaucouleurs supergalactic
     coordinates. [#3892]
 
+  - Added ecliptic coordinates, including ``GeocentricTrueEcliptic``, 
+    ``BarycentricTrueEcliptic``, and ``HeliocentricTrueEcliptic``. [#3749]
+
 - ``astropy.cosmology``
 
   - Add Planck 2015 cosmology [#3476]

diff --git a/astropy/coordinates/__init__.py b/astropy/coordinates/__init__.py
@@ -21,4 +21,15 @@
 from .sky_coordinate import *
 from .funcs import *
 
-__doc__ += builtin_frames._transform_graph_docs
+__doc__ += builtin_frames._transform_graph_docs + """
+
+.. note::
+
+    The ecliptic coordinate systems (added in Astropy v1.1) have not been
+    extensively tested for accuracy or consistency with other implementations of
+    ecliptic coordinates.  We welcome contributions to add such testing, but in
+    the meantime, users who depend on consistency with other implementations may
+    wish to check test inputs against good datasets before using Astropy's
+    ecliptic coordinates.
+
+"""
diff --git a/astropy/coordinates/builtin_frames/__init__.py b/astropy/coordinates/builtin_frames/__init__.py
@@ -34,6 +34,7 @@
 from .gcrs import GCRS, PrecessedGeocentric
 from .cirs import CIRS
 from .itrs import ITRS
+from .ecliptic import GeocentricTrueEcliptic, BarycentricTrueEcliptic, HeliocentricTrueEcliptic
 
 #need to import transformations so that they get registered in the graph
 from . import icrs_fk5_transforms
@@ -43,11 +44,13 @@
 from . import icrs_cirs_transforms
 from . import cirs_observed_transforms
 from . import intermediate_rotation_transforms
+from . import ecliptic_transforms
 
 # we define an __all__ because otherwise the transformation modules get included
 __all__ = ['ICRS', 'FK5', 'FK4', 'FK4NoETerms', 'Galactic', 'Galactocentric',
            'Supergalactic', 'AltAz', 'GCRS', 'CIRS', 'ITRS',
-           'PrecessedGeocentric']
+           'PrecessedGeocentric', 'GeocentricTrueEcliptic',
+           'BarycentricTrueEcliptic', 'HeliocentricTrueEcliptic']
 
 def _make_transform_graph_docs():
     """

diff --git a/astropy/coordinates/builtin_frames/ecliptic.py b/astropy/coordinates/builtin_frames/ecliptic.py
@@ -0,0 +1,159 @@
+# -*- coding: utf-8 -*-
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import (absolute_import, unicode_literals, division,
+                        print_function)
+
+from ..representation import SphericalRepresentation
+from ..baseframe import BaseCoordinateFrame, RepresentationMapping, TimeFrameAttribute
+from .utils import EQUINOX_J2000
+
+
+__all__ = ['GeocentricTrueEcliptic', 'BarycentricTrueEcliptic', 'HeliocentricTrueEcliptic']
+
+
+class GeocentricTrueEcliptic(BaseCoordinateFrame):
+    """
+    Geocentric ecliptic coordinates.  These origin of the coordinates are the
+    geocenter (Earth), with the x axis pointing to the *true* (not mean) equinox
+    at the time specified by the ``equinox`` attribute, and the xy-plane in the
+    plane of the ecliptic for that date.
+
+    Be aware that the definition of "geocentric" here means that this frame
+    *includes* light deflection from the sun, abberation, etc when transfoming
+    to/from e.g. ICRS.
+
+    This frame has one frame attribute:
+
+    * ``equinox``
+        The date to assume for this frame.  Determines the location of the
+        x-axis and the location of the Earth (necessary for transformation to
+        non-geocentric systems).
+
+    .. warning::
+        In the current version of astropy, the ecliptic frames do not yet have
+        stringent accuracy tests.  We recommend you test to "known-good" cases
+        to ensure this frames are what you are looking for. (and then ideally
+        you would contribute these tests to Astropy!)
+
+    Parameters
+    ----------
+    representation : `BaseRepresentation` or None
+        A representation object or None to have no data (or use the other keywords)
+    lambda : `Angle`, optional, must be keyword
+        The ecliptic longitude for this object (``beta`` must also be given and
+        ``representation`` must be None).
+    beta : `Angle`, optional, must be keyword
+        The ecliptic latitde for this object (``lambda`` must also be given and
+        ``representation`` must be None).
+    delta : `~astropy.units.Quantity`, optional, must be keyword
+        The Distance for this object from the geocenter.
+        (``representation`` must be None).
+    """
+
+    frame_specific_representation_info = {
+        'unitspherical': [RepresentationMapping('lon', 'lambda'),
+                          RepresentationMapping('lat', 'beta')]
+    }
+    frame_specific_representation_info['spherical'] = \
+        frame_specific_representation_info['unitspherical'] + \
+        [RepresentationMapping('distance', 'delta')]
+
+    default_representation = SphericalRepresentation
+
+    equinox = TimeFrameAttribute(default=EQUINOX_J2000)
+
+
+class BarycentricTrueEcliptic(BaseCoordinateFrame):
+    """
+    Barycentric ecliptic coordinates.  These origin of the coordinates are the
+    barycenter of the solar system, with the x axis pointing in the direction of
+    the *true* (not mean) equinox as at the time specified by the ``equinox``
+    attribute (as seen from Earth), and the xy-plane in the plane of the
+    ecliptic for that date.
+
+    This frame has one frame attribute:
+
+    * ``equinox``
+        The date to assume for this frame.  Determines the location of the
+        x-axis and the location of the Earth and Sun.
+
+    .. warning::
+        In the current version of astropy, the ecliptic frames do not yet have
+        stringent accuracy tests.  We recommend you test to "known-good" cases
+        to ensure this frames are what you are looking for. (and then ideally
+        you would contribute these tests to Astropy!)
+
+    Parameters
+    ----------
+    representation : `BaseRepresentation` or None
+        A representation object or None to have no data (or use the other keywords)
+    l : `Angle`, optional, must be keyword
+        The ecliptic longitude for this object (``b`` must also be given and
+        ``representation`` must be None).
+    b : `Angle`, optional, must be keyword
+        The ecliptic latitde for this object (``l`` must also be given and
+        ``representation`` must be None).
+    r : `~astropy.units.Quantity`, optional, must be keyword
+        The Distance for this object from the sun's center.
+        (``representation`` must be None).
+    """
+
+    frame_specific_representation_info = {
+        'unitspherical': [RepresentationMapping('lon', 'l'),
+                          RepresentationMapping('lat', 'b')]
+    }
+    frame_specific_representation_info['spherical'] = \
+        frame_specific_representation_info['unitspherical'] + \
+        [RepresentationMapping('distance', 'r')]
+
+    default_representation = SphericalRepresentation
+
+    equinox = TimeFrameAttribute(default=EQUINOX_J2000)
+
+
+class HeliocentricTrueEcliptic(BaseCoordinateFrame):
+    """
+    Heliocentric ecliptic coordinates.  These origin of the coordinates are the
+    center of the sun, with the x axis pointing in the direction of
+    the *true* (not mean) equinox as at the time specified by the ``equinox``
+    attribute (as seen from Earth), and the xy-plane in the plane of the
+    ecliptic for that date.
+
+    This frame has one frame attribute:
+
+    * ``equinox``
+        The date to assume for this frame.  Determines the location of the
+        x-axis and the location of the Earth and Sun.
+
+    .. warning::
+        In the current version of astropy, the ecliptic frames do not yet have
+        stringent accuracy tests.  We recommend you test to "known-good" cases
+        to ensure this frames are what you are looking for. (and then ideally
+        you would contribute these tests to Astropy!)
+
+    Parameters
+    ----------
+    representation : `BaseRepresentation` or None
+        A representation object or None to have no data (or use the other keywords)
+    l : `Angle`, optional, must be keyword
+        The ecliptic longitude for this object (``b`` must also be given and
+        ``representation`` must be None).
+    b : `Angle`, optional, must be keyword
+        The ecliptic latitde for this object (``l`` must also be given and
+        ``representation`` must be None).
+    r : `~astropy.units.Quantity`, optional, must be keyword
+        The Distance for this object from the sun's center.
+        (``representation`` must be None).
+    """
+
+    frame_specific_representation_info = {
+        'unitspherical': [RepresentationMapping('lon', 'l'),
+                          RepresentationMapping('lat', 'b')]
+    }
+    frame_specific_representation_info['spherical'] = \
+        frame_specific_representation_info['unitspherical'] + \
+        [RepresentationMapping('distance', 'r')]
+
+    default_representation = SphericalRepresentation
+
+    equinox = TimeFrameAttribute(default=EQUINOX_J2000)
diff --git a/astropy/coordinates/builtin_frames/ecliptic_transforms.py b/astropy/coordinates/builtin_frames/ecliptic_transforms.py
@@ -0,0 +1,97 @@
+# -*- coding: utf-8 -*-
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+Contains the transformation functions for getting to/from ecliptic systems.
+"""
+from __future__ import (absolute_import, unicode_literals, division,
+                        print_function)
+
+import numpy as np
+
+from ... import units as u
+from ..baseframe import frame_transform_graph
+from ..transformations import FunctionTransform, DynamicMatrixTransform
+from ..angles import rotation_matrix
+from ..representation import CartesianRepresentation
+from ... import _erfa as erfa
+
+from .icrs import ICRS
+from .gcrs import GCRS
+from .ecliptic import GeocentricTrueEcliptic, BarycentricTrueEcliptic, HeliocentricTrueEcliptic
+from .utils import cartrepr_from_matmul
+
+
+def _ecliptic_rotation_matrix(equinox):
+    rnpb = erfa.pnm06a(equinox.jd1, equinox.jd2)
+    obl = erfa.obl06(equinox.jd1, equinox.jd2)*u.radian
+    return np.asarray(np.dot(rotation_matrix(obl, 'x'), rnpb))
+
+
+@frame_transform_graph.transform(FunctionTransform, GCRS, GeocentricTrueEcliptic)
+def gcrs_to_geoecliptic(from_coo, to_frame):
+    if np.any(from_coo.obstime != to_frame.equinox):
+        # if they GCRS obstime and ecliptic equinox are not the same, we first
+        # have to move to a GCRS where they are.
+        frameattrs = from_coo.get_frame_attr_names()
+        frameattrs['obstime'] = to_frame.equinox
+        from_coo = from_coo.transform_to(GCRS(**frameattrs))
+
+    rmat = _ecliptic_rotation_matrix(to_frame.equinox)
+    newrepr = cartrepr_from_matmul(rmat, from_coo)
+    return to_frame.realize_frame(newrepr)
+
+
+@frame_transform_graph.transform(FunctionTransform, GeocentricTrueEcliptic, GCRS)
+def geoecliptic_to_gcrs(from_coo, to_frame):
+    rmat = _ecliptic_rotation_matrix(to_frame.equinox)
+    newrepr = cartrepr_from_matmul(rmat, from_coo, transpose=True)
+
+    if np.all(from_coo.equinox == to_frame.obstime):
+        return to_frame.realize_frame(newrepr)
+    else:
+        # if the GCRS obstime and ecliptic equinox don't match, need to move
+        # to one where they do
+        frameattrs = to_frame.get_frame_attr_names()
+        frameattrs['obstime'] = from_coo.equinox
+        return GCRS(newrepr, **frameattrs).transform_to(to_frame)
+
+
+@frame_transform_graph.transform(DynamicMatrixTransform, ICRS, BarycentricTrueEcliptic)
+def icrs_to_baryecliptic(from_coo, to_frame):
+    return _ecliptic_rotation_matrix(to_frame.equinox)
+
+
+@frame_transform_graph.transform(DynamicMatrixTransform, BarycentricTrueEcliptic, ICRS)
+def baryecliptic_to_icrs(from_coo, to_frame):
+    return icrs_to_baryecliptic(to_frame, from_coo).T
+
+
+@frame_transform_graph.transform(FunctionTransform, ICRS, HeliocentricTrueEcliptic)
+def icrs_to_helioecliptic(from_coo, to_frame):
+    pvh, pvb = erfa.epv00(to_frame.equinox.jd1, to_frame.equinox.jd2)
+    delta_bary_to_helio = pvh[..., 0, :] - pvb[..., 0, :]
+
+    #first offset to heliocentric
+    heliocart = from_coo.cartesian.xyz + delta_bary_to_helio * u.au
+
+    # now compute the matrix to precess to the right orientation
+    rmat = _ecliptic_rotation_matrix(to_frame.equinox)
+
+    # it's not really ICRS because of the offset, but this is digestible by cartrepr_from_matmul
+    newrepr = cartrepr_from_matmul(rmat, ICRS(CartesianRepresentation(heliocart)))
+    return to_frame.realize_frame(newrepr)
+
+
+@frame_transform_graph.transform(FunctionTransform, HeliocentricTrueEcliptic, ICRS)
+def helioecliptic_to_icrs(from_coo, to_frame):
+
+    # first un-precess from ecliptic to ICRS orientation
+    rmat = _ecliptic_rotation_matrix(from_coo.equinox)
+    intermed_repr = cartrepr_from_matmul(rmat, from_coo, transpose=True)
+
+    # now offset back to barycentric, which is the correct center for ICRS
+    pvh, pvb = erfa.epv00(from_coo.equinox.jd1, from_coo.equinox.jd2)
+    delta_bary_to_helio = pvh[..., 0, :] - pvb[..., 0, :]
+    newrepr = CartesianRepresentation(intermed_repr.cartesian.xyz - delta_bary_to_helio*u.au)
+
+    return to_frame.realize_frame(newrepr)
diff --git a/astropy/coordinates/builtin_frames/intermediate_rotation_transforms.py b/astropy/coordinates/builtin_frames/intermediate_rotation_transforms.py
@@ -12,7 +12,7 @@
 
 from ..baseframe import frame_transform_graph
 from ..transformations import FunctionTransform
-from ..representation import CartesianRepresentation
+from .utils import cartrepr_from_matmul
 from ... import _erfa as erfa
 
 from .gcrs import GCRS, PrecessedGeocentric
@@ -54,28 +54,6 @@ def gcrs_precession_mat(equinox):
     return erfa.fw2m(gamb, phib, psib, epsa)
 
 
-def cartrepr_from_matmul(pmat, coo, transpose=False):
-    if pmat.shape[-2:] != (3, 3):
-        raise ValueError("tried to do matrix multiplication with an array that "
-                         "doesn't end in 3x3")
-    if coo.isscalar:
-        # a simpler path for scalar coordinates
-        if transpose:
-            pmat = pmat.T
-        newxyz = np.sum(pmat * coo.cartesian.xyz, axis=-1)
-    else:
-        xyz = coo.cartesian.xyz.T
-        # these expression are the same as iterating over the first dimension of
-        # pmat and xyz and doing matrix multiplication on each in turn.  resulting
-        # dimension is <coo shape> x 3
-        pmat = pmat.reshape(pmat.size//9, 3, 3)
-        if transpose:
-            pmat = pmat.transpose(0, 2, 1)
-        newxyz = np.sum(pmat * xyz.reshape(xyz.size//3, 1, 3), axis=-1).T
-
-    return CartesianRepresentation(newxyz)
-
-
 # now the actual transforms
 
 # the priority for the GCRS<->ITRS trasnsforms are higher (=less traveled) to

diff --git a/astropy/coordinates/builtin_frames/utils.py b/astropy/coordinates/builtin_frames/utils.py
@@ -15,6 +15,7 @@
 from ...time import Time
 from ...utils import iers
 from ...utils.exceptions import AstropyWarning
+from ..representation import CartesianRepresentation
 
 # The UTC time scale is not properly defined prior to 1960, so Time('B1950',
 # scale='utc') will emit a warning. Instead, we use Time('B1950', scale='tai')
@@ -41,6 +42,31 @@
 """
 
 
+def cartrepr_from_matmul(pmat, coo, transpose=False):
+    """
+    Note that pmat should be an ndarray, *not* a matrix.
+    """
+    if pmat.shape[-2:] != (3, 3):
+        raise ValueError("tried to do matrix multiplication with an array that "
+                         "doesn't end in 3x3")
+    if coo.isscalar:
+        # a simpler path for scalar coordinates
+        if transpose:
+            pmat = pmat.T
+        newxyz = np.sum(pmat * coo.cartesian.xyz, axis=-1)
+    else:
+        xyz = coo.cartesian.xyz.T
+        # these expression are the same as iterating over the first dimension of
+        # pmat and xyz and doing matrix multiplication on each in turn.  resulting
+        # dimension is <coo shape> x 3
+        pmat = pmat.reshape(pmat.size//9, 3, 3)
+        if transpose:
+            pmat = pmat.transpose(0, 2, 1)
+        newxyz = np.sum(pmat * xyz.reshape(xyz.size//3, 1, 3), axis=-1).T
+
+    return CartesianRepresentation(newxyz)
+
+
 def get_polar_motion(time):
     """
     gets the two polar motion components in radians for use with apio13