Add a GeocentricSolarEcliptic Frame

src/poliastro/frames.py

@@ -5,24 +5,31 @@
 from typing import Dict
 
 import numpy as np
-from astropy import _erfa, units as u
+from astropy import _erfa as erfa, units as u
 from astropy.coordinates import (
     GCRS,
     ICRS,
     AffineTransform,
     BaseEclipticFrame,
     BaseRADecFrame,
     CartesianDifferential,
+    CartesianRepresentation,
     FunctionTransformWithFiniteDifference,
     TimeAttribute,
     UnitSphericalRepresentation,
     frame_transform_graph,
+    get_body,
     get_body_barycentric,
     get_body_barycentric_posvel,
 )
 from astropy.coordinates.baseframe import FrameMeta
 from astropy.coordinates.builtin_frames.utils import DEFAULT_OBSTIME, get_jd12
-from astropy.coordinates.matrix_utilities import matrix_transpose, rotation_matrix
+from astropy.coordinates.matrix_utilities import (
+    matrix_product,
+    matrix_transpose,
+    rotation_matrix,
+)
+from astropy.coordinates.transformations import DynamicMatrixTransform
 
 from poliastro.bodies import (
     Earth,
@@ -60,7 +67,7 @@ class HeliocentricEclipticJ2000(BaseEclipticFrame):
 
 def _ecliptic_rotation_matrix():
     jd1, jd2 = get_jd12(J2000, J2000.scale)
-    obl = _erfa.obl80(jd1, jd2) * u.radian
+    obl = erfa.obl80(jd1, jd2) * u.radian
     assert obl.to(u.arcsec).value == 84381.448
     return rotation_matrix(obl, "x")
 
@@ -220,6 +227,81 @@ class PlutoICRS(_PlanetaryICRS):
     body = Pluto
 
 
+def _make_rotation_matrix_from_reprs(start_representation, end_representation):
+    """
+    Return the matrix for the direct rotation from one representation to a second representation.
+    The representations need not be normalized first.
+    """
+    A = start_representation.to_cartesian()
+    B = end_representation.to_cartesian()
+    rotation_axis = A.cross(B)
+    rotation_angle = -np.arccos(
+        A.dot(B) / (A.norm() * B.norm())
+    )  # negation is required
+
+    # This line works around some input/output quirks of Astropy's rotation_matrix()
+    matrix = np.array(rotation_matrix(rotation_angle, rotation_axis.xyz.value.tolist()))
+    return matrix
+
+
+def _obliquity_rotation_value(equinox):
+    """
+    Function to calculate obliquity of the earth.
+    This uses obl06 of erfa.
+    """
+    jd1, jd2 = get_jd12(equinox, "tt")
+    obl = erfa.obl06(jd1, jd2) * u.radian
+    return obl.to(u.deg)
+
+
+class GeocentricSolarEcliptic(BaseEclipticFrame):
+    """
+    This system has its X axis towards the Sun and its Z axis perpendicular to
+    the plane of the Earth's orbit around the Sun (positive North). This system
+    is fixed with respect to the Earth-Sun line. It is convenient for specifying
+    magnetospheric boundaries. It has also been widely adopted as the system for
+    representing vector quantities in space physics databases.
+
+    """
+
+    obstime = TimeAttribute(default=DEFAULT_OBSTIME)
+
+
+@frame_transform_graph.transform(DynamicMatrixTransform, GCRS, GeocentricSolarEcliptic)
+def gcrs_to_geosolarecliptic(gcrs_coo, to_frame):
+
+    if not to_frame.obstime.isscalar:
+        raise ValueError(
+            "To perform this transformation the obstime Attribute must be a scalar."
+        )
+
+    _earth_orbit_perpen_point_gcrs = UnitSphericalRepresentation(
+        lon=0 * u.deg, lat=(90 * u.deg - _obliquity_rotation_value(to_frame.obstime))
+    )
+
+    _earth_detilt_matrix = _make_rotation_matrix_from_reprs(
+        _earth_orbit_perpen_point_gcrs, CartesianRepresentation(0, 0, 1)
+    )
+
+    sun_pos_gcrs = get_body("sun", to_frame.obstime).cartesian
+    earth_pos_gcrs = get_body("earth", to_frame.obstime).cartesian
+    sun_earth = sun_pos_gcrs - earth_pos_gcrs
+
+    sun_earth_detilt = sun_earth.transform(_earth_detilt_matrix)
+
+    # Earth-Sun Line in Geocentric Solar Ecliptic Frame
+    x_axis = CartesianRepresentation(1, 0, 0)
+
+    rot_matrix = _make_rotation_matrix_from_reprs(sun_earth_detilt, x_axis)
+
+    return matrix_product(rot_matrix, _earth_detilt_matrix)
+
+
+@frame_transform_graph.transform(DynamicMatrixTransform, GeocentricSolarEcliptic, GCRS)
+def geosolarecliptic_to_gcrs(from_coo, gcrs_frame):
+    return matrix_transpose(gcrs_to_geosolarecliptic(gcrs_frame, from_coo))
+
+
 _FRAME_MAPPING = {
     Sun: {Planes.EARTH_EQUATOR: HCRS, Planes.EARTH_ECLIPTIC: HeliocentricEclipticJ2000},
     Mercury: {Planes.EARTH_EQUATOR: MercuryICRS},

src/poliastro/tests/test_frames.py

@@ -6,6 +6,7 @@
     solar_system_ephemeris,
 )
 from astropy.tests.helper import assert_quantity_allclose
+from astropy.time import Time
 
 from poliastro.bodies import (
     Earth,
@@ -24,6 +25,7 @@
     GCRS,
     HCRS,
     ICRS,
+    GeocentricSolarEcliptic,
     JupiterICRS,
     MarsICRS,
     MercuryICRS,
@@ -109,3 +111,30 @@ def test_icrs_body_position_to_planetary_frame_yields_zeros(body, frame):
         )
 
     assert_quantity_allclose(vector_result.xyz, [0, 0, 0] * u.km, atol=1e-7 * u.km)
+
+
+def test_round_trip_from_GeocentricSolarEcliptic_gives_same_results():
+    gcrs = GCRS(ra="02h31m49.09s", dec="+89d15m50.8s", distance=200 * u.km)
+    gse = gcrs.transform_to(GeocentricSolarEcliptic(obstime=Time("J2000")))
+    gcrs_back = gse.transform_to(GCRS(obstime=Time("J2000")))
+    assert_quantity_allclose(gcrs_back.dec.value, gcrs.dec.value, atol=1e-7)
+    assert_quantity_allclose(gcrs_back.ra.value, gcrs.ra.value, atol=1e-7)
+
+
+def test_GeocentricSolarEcliptic_against_data():
+    gcrs = GCRS(ra="02h31m49.09s", dec="+89d15m50.8s", distance=200 * u.km)
+    gse = gcrs.transform_to(GeocentricSolarEcliptic(obstime=Time("J2000")))
+    lon = 233.11663895663975
+    lat = 48.64652559835358
+    assert_quantity_allclose(gse.lat.value, lat, atol=1e-7)
+    assert_quantity_allclose(gse.lon.value, lon, atol=1e-7)
+
+
+def test_GeocentricSolarEcliptic_raises_error_nonscalar_obstime():
+    with pytest.raises(ValueError) as excinfo:
+        gcrs = GCRS(ra="02h31m49.09s", dec="+89d15m50.8s", distance=200 * u.km)
+        gcrs.transform_to(GeocentricSolarEcliptic(obstime=Time(["J3200", "J2000"])))
+    assert (
+        "To perform this transformation the "
+        "obstime Attribute must be a scalar." in str(excinfo.value)
+    )