From 5a85654e40fcb63136e17cce6b5f829cd4a831d7 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Thu, 7 Mar 2024 23:07:54 -0500
Subject: [PATCH 01/10] start supporting lunar ellipsoid models

---
 lunarsky/moon.py                  | 83 ++++++++++++++++++++++++++-----
 lunarsky/tests/test_transforms.py |  9 +++-
 2 files changed, 78 insertions(+), 14 deletions(-)

diff --git a/lunarsky/moon.py b/lunarsky/moon.py
index d5ee239..c8ed0f0 100644
--- a/lunarsky/moon.py
+++ b/lunarsky/moon.py
@@ -4,6 +4,7 @@
 from astropy.units.quantity import QuantityInfoBase
 from astropy.coordinates.angles import Longitude, Latitude
 from astropy.coordinates.earth import GeodeticLocation
+from astropy.coordinates.representation.geodetic import BaseGeodeticRepresentation
 from astropy.coordinates.representation import (
     CartesianRepresentation,
     SphericalRepresentation,
@@ -14,6 +15,53 @@
 
 __all__ = ["MoonLocation", "MoonLocationAttribute"]
 
+class SPHERESelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid as a sphere
+
+        Radius defined by lunarsky.spice_utils.LUNAR_RADIUS
+    """
+    _equatorial_radius = LUNAR_RADIUS * u.m 
+    _flattening = 0.0
+
+class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid defined by gravimetry of GRAIL data.
+
+    https://doi.org/10.1007/s40328-023-00415-w
+    """
+    _equatorial_radius = 1737576.6 * u.m
+    _flattening = 0.000305
+
+class CE1LAM10SelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid from Chang'e 1 laser altimetry.
+
+    Rotation ellipsoid = CE-1-LAM-GEO
+
+    https://doi.org/10.1007/s11430-010-4060-6
+    """
+    _equatorial_radius = 1737.632 * u.km
+    _flattening = 1/973.463
+
+
+# Define reference ellipsoids
+LUNAR_ELLIPSOIDS = {
+    "SPHERE" : SPHERESelenodeticRepresentation,
+    "GRAIL23" : GRAIL23SelenodeticRepresentation,
+    "CE-1-LAM-GEO" : CE1LAM10SelenodeticRepresentation,
+}
+
+
+class SelenodeticLocation(GeodeticLocation):
+    """Rename GeodeticLocation class for clarity"""
+
+
+def _check_ellipsoid(ellipsoid=None, default="SPHERE"):
+    """Defaulting lunar ellipsoid"""
+    if ellipsoid is None:
+        ellipsoid = default
+    if ellipsoid not in LUNAR_ELLIPSOIDS:
+        raise ValueError(f"Ellipsoid {ellipsoid} not among known ones ({ELLIPSOIDS})")
+    return ellipsoid
+
 
 class MoonLocationInfo(QuantityInfoBase):
     """
@@ -22,10 +70,12 @@ class MoonLocationInfo(QuantityInfoBase):
     be used as a general way to store meta information.
     """
 
-    _represent_as_dict_attrs = ("x", "y", "z")
+    _represent_as_dict_attrs = ("x", "y", "z", "ellipsoid")
 
     def _construct_from_dict(self, map):
+        ellipsoid = map.pop("ellipsoid")
         out = self._parent_cls(**map)
+        out.ellipsoid = ellipsoid
         return out
 
     def new_like(self, cols, length, metadata_conflicts="warn", name=None):
@@ -53,7 +103,7 @@ def new_like(self, cols, length, metadata_conflicts="warn", name=None):
             Empty instance of this class consistent with ``cols``
         """
         # Very similar to QuantityInfo.new_like, but the creation of the
-        # map is different enough that this needs its own rouinte.
+        # map is different enough that this needs its own routine.
         # Get merged info attributes shape, dtype, format, description.
         attrs = self.merge_cols_attributes(
             cols, metadata_conflicts, name, ("meta", "format", "description")
@@ -110,6 +160,7 @@ class MoonLocation(u.Quantity):
     property.
     """
 
+    _ellipsoid = "SPHERE"
     _location_dtype = np.dtype({"names": ["x", "y", "z"], "formats": [np.float64] * 3})
     _array_dtype = np.dtype((np.float64, (3,)))
 
@@ -248,7 +299,7 @@ def from_selenocentric(cls, x, y, z, unit=None):
         return inst
 
     @classmethod
-    def from_selenodetic(cls, lon, lat, height=0.0):
+    def from_selenodetic(cls, lon, lat, height=0.0, ellipsoid=None):
         """
         Location on the Moon, from latitude and longitude.
 
@@ -264,6 +315,10 @@ def from_selenodetic(cls, lon, lat, height=0.0):
             Height above reference sphere (if float, in meters; default: 0).
             The reference sphere is a sphere of radius 1737.1 kilometers,
             from the center of mass of the Moon.
+        ellipsoid : str, optional
+            Name of the reference ellipsoid to use (default: 'SPHERE').
+            Available ellipsoids are:  'SPHERE', 'GRAIL23', 'CE-1-LAM-GEO'.
+            See docstrings for classes in ELLIPSOIDS dictionary for references.
 
         Raises
         ------
@@ -281,7 +336,11 @@ def from_selenodetic(cls, lon, lat, height=0.0):
         relative to a prime meridian, which is itself given by
         the mean position of the "sub-Earth" point on the lunar surface.
 
+        For the conversion to selenocentric coordinates, the ERFA routine
+        ``gd2gce`` is used.  See https://github.com/liberfa/erfa
+
         """
+        ellipsoid = _check_ellipsoid(ellipsoid, default=cls._ellipsoid)
         lon = Longitude(lon, u.degree, copy=False).wrap_at(180 * u.degree)
         lat = Latitude(lat, u.degree, copy=False)
         # don't convert to m by default, so we can use the height unit below.
@@ -351,11 +410,10 @@ def selenodetic(self):
         """Convert to selenodetic coordinates."""
         return self.to_selenodetic()
 
-    def to_selenodetic(self):
+    def to_selenodetic(self, ellipsoid=None):
         """Convert to selenodetic coordinates (lat, lon, height).
 
-        Height is in reference to a sphere with radius `_lunar_radius`,
-        centered at the center of mass.
+        Height is in reference to the ellipsoid.
 
         Returns
         -------
@@ -364,15 +422,16 @@ def to_selenodetic(self):
             `~astropy.coordinates.Latitude`, and `~astropy.units.Quantity`
 
         """
+        ellipsoid = _check_ellipsoid(ellipsoid)
         xyz = self.view(self._array_dtype, u.Quantity)
-        lld = CartesianRepresentation(xyz, xyz_axis=-1, copy=False).represent_as(
-            SphericalRepresentation
+        llh = CartesianRepresentation(xyz, xyz_axis=-1, copy=False).represent_as(
+            LUNAR_ELLIPSOIDS[ellipsoid],
         )
-        return GeodeticLocation(
-            Longitude(lld.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
-            Latitude(lld.lat, u.degree, copy=False),
+        return SelenodeticLocation(
+            Longitude(llh.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
+            Latitude(llh.lat, u.degree, copy=False),
             u.Quantity(
-                lld.distance - (self._lunar_radius * u.meter), self.unit, copy=False
+                llh.height, self.unit, copy=False
             ),
         )
 
diff --git a/lunarsky/tests/test_transforms.py b/lunarsky/tests/test_transforms.py
index 1ead075..39befe7 100644
--- a/lunarsky/tests/test_transforms.py
+++ b/lunarsky/tests/test_transforms.py
@@ -8,6 +8,11 @@
 import lunarsky
 import pytest
 
+try:
+    from astropy.coordinates.angles.utils import angular_separation
+except ImportError:
+    from astropy.coordinates.angle_utilities import angular_separation
+
 # Lunar station positions
 Nangs = 3
 latitudes = np.linspace(-90, 90, Nangs)
@@ -97,7 +102,7 @@ def test_mcmf_to_mcmf():
     )
     sph0 = src.spherical
     sph1 = orig_pos.spherical
-    res = ac.angle_utilities.angular_separation(
+    res = angular_separation(
         sph0.lon, sph0.lat, sph1.lon, sph1.lat
     ).to("deg")
     assert_quantity_allclose(res, 177 * un.deg, atol=5 * un.deg)
@@ -238,7 +243,7 @@ def test_incompatible_shape_error(Nt, Nl, success):
     if success:
         lunarsky.LunarTopo(location=locs, obstime=times)
     else:
-        with pytest.raises(ValueError, match="non-scalar attributes"):
+        with pytest.raises(ValueError, match="non-scalar data and/or attributes"):
             lunarsky.LunarTopo(location=locs, obstime=times)
 
 

From 48f52ce49005e36dc400caccec0dcb6f93e0b54c Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Fri, 8 Mar 2024 15:56:34 -0500
Subject: [PATCH 02/10] add Matt's ellipsoid

---
 lunarsky/moon.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/lunarsky/moon.py b/lunarsky/moon.py
index c8ed0f0..f9924eb 100644
--- a/lunarsky/moon.py
+++ b/lunarsky/moon.py
@@ -23,6 +23,14 @@ class SPHERESelenodeticRepresentation(BaseGeodeticRepresentation):
     _equatorial_radius = LUNAR_RADIUS * u.m 
     _flattening = 0.0
 
+class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid from NASA/GSFC "Planetary Fact Sheet"
+
+       https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
+    """
+    _equatorial_radius = 1738.1 * u.m
+    _flattening = 0.0012
+
 class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation):
     """Lunar ellipsoid defined by gravimetry of GRAIL data.
 
@@ -45,6 +53,7 @@ class CE1LAM10SelenodeticRepresentation(BaseGeodeticRepresentation):
 # Define reference ellipsoids
 LUNAR_ELLIPSOIDS = {
     "SPHERE" : SPHERESelenodeticRepresentation,
+    "GSFC" : GSFCSelenodeticRepresentation,
     "GRAIL23" : GRAIL23SelenodeticRepresentation,
     "CE-1-LAM-GEO" : CE1LAM10SelenodeticRepresentation,
 }

From c0cb7823eee5dfc28da9bc59bf8f23c17ba88e19 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Tue, 12 Mar 2024 12:58:00 -0400
Subject: [PATCH 03/10] fix errors with maintaining ellipsoids through
 transformations; add some tests

---
 lunarsky/moon.py                  | 68 +++++++++++++++++++------------
 lunarsky/spice_utils.py           | 20 ++-------
 lunarsky/tests/test_transforms.py | 44 +++++++++++---------
 lunarsky/topo.py                  | 29 +++++--------
 4 files changed, 80 insertions(+), 81 deletions(-)

diff --git a/lunarsky/moon.py b/lunarsky/moon.py
index f9924eb..0a93833 100644
--- a/lunarsky/moon.py
+++ b/lunarsky/moon.py
@@ -11,7 +11,9 @@
 )
 from astropy.coordinates.attributes import Attribute
 
-from .spice_utils import remove_topo, LUNAR_RADIUS
+from .spice_utils import remove_topo
+
+LUNAR_RADIUS = 1737.1e3  # m
 
 __all__ = ["MoonLocation", "MoonLocationAttribute"]
 
@@ -28,7 +30,7 @@ class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation):
 
        https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
     """
-    _equatorial_radius = 1738.1 * u.m
+    _equatorial_radius = 1738.1e3 * u.m
     _flattening = 0.0012
 
 class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation):
@@ -51,7 +53,7 @@ class CE1LAM10SelenodeticRepresentation(BaseGeodeticRepresentation):
 
 
 # Define reference ellipsoids
-LUNAR_ELLIPSOIDS = {
+SELENOIDS = {
     "SPHERE" : SPHERESelenodeticRepresentation,
     "GSFC" : GSFCSelenodeticRepresentation,
     "GRAIL23" : GRAIL23SelenodeticRepresentation,
@@ -67,11 +69,21 @@ def _check_ellipsoid(ellipsoid=None, default="SPHERE"):
     """Defaulting lunar ellipsoid"""
     if ellipsoid is None:
         ellipsoid = default
-    if ellipsoid not in LUNAR_ELLIPSOIDS:
-        raise ValueError(f"Ellipsoid {ellipsoid} not among known ones ({ELLIPSOIDS})")
+    if ellipsoid not in SELENOIDS:
+        raise ValueError(f"Ellipsoid {ellipsoid} not among known ones ({SELENOIDS})")
     return ellipsoid
 
 
+def _add_rm_ellipsoid(ellipsoid, name, add=True):
+    """Add or remove an ellipsoid class from SELENOIDS.
+    Needed for testing."""
+    global SELENOIDS
+    if add:
+        SELENOIDS[name] = ellipsoid
+    else:
+        SELENOIDS.pop(name, None)
+
+
 class MoonLocationInfo(QuantityInfoBase):
     """
     Container for meta information like name, description, format.  This is
@@ -152,8 +164,10 @@ class MoonLocation(u.Quantity):
     This class uses the ME frame.
 
     Positions may be defined in Cartesian (x, y, z) coordinates with respect to the
-    center of mass of the Moon, or in ``selenodetic'' coordinates (longitude, latitude).
-    In selenodetic coordinates, positions are on the surface exactly.
+    center of mass of the Moon, or in ``selenodetic'' coordinates (longitude, latitude, height).
+
+    Selenodetic coordinates are defined with respect to a reference ellipsoid. The default is a
+    sphere of radius 1731.1e3 km, but other ellipsoids are available. See lunarsky.moon.SELENOIDS.
 
     See:
         "A Standardized Lunar Coordinate System for the Lunar Reconnaissance
@@ -173,8 +187,6 @@ class MoonLocation(u.Quantity):
     _location_dtype = np.dtype({"names": ["x", "y", "z"], "formats": [np.float64] * 3})
     _array_dtype = np.dtype((np.float64, (3,)))
 
-    _lunar_radius = LUNAR_RADIUS
-
     # Manage the set of defined ephemerides.
     # Class attributes only
     _inuse_stat_ids = []
@@ -230,7 +242,7 @@ def _set_site_id(cls, inst):
             raise ValueError("Too many unique MoonLocation objects open at once.")
 
         for llh in llh_arr:
-            lonlatheight = "_".join(["{:.4f}".format(ll) for ll in llh])
+            lonlatheight = "_".join(["{:.4f}".format(ll) for ll in llh] + [inst._ellipsoid])
             if lonlatheight not in cls._existing_locs:
                 new_stat_id = cls._avail_stat_ids.pop()
                 cls._existing_locs.append(lonlatheight)
@@ -362,23 +374,14 @@ def from_selenodetic(cls, lon, lat, height=0.0, ellipsoid=None):
                     str(lon.shape), str(lat.shape)
                 )
             )
-        # get selenocentric coordinates. Have to give one-dimensional array.
-
-        lunar_radius = u.Quantity(cls._lunar_radius, u.m, copy=False)
-
-        Npts = lon.size
-        xyz = np.zeros((Npts, 3))
-        xyz[:, 0] = (lunar_radius + height) * np.cos(lat) * np.cos(lon)
-        xyz[:, 1] = (lunar_radius + height) * np.cos(lat) * np.sin(lon)
-        xyz[:, 2] = (lunar_radius + height) * np.sin(lat)
 
-        xyz = np.squeeze(xyz)
-
-        self = xyz.ravel().view(cls._location_dtype, cls).reshape(xyz.shape[:-1])
-        self._unit = u.meter
-        inst = self.to(height.unit)
+        # get selenocentric coordinates. Have to give one-dimensional array.
 
-        return inst
+        selenodetic = SELENOIDS[ellipsoid](lon, lat, height, copy=False)
+        xyz = selenodetic.to_cartesian().get_xyz(xyz_axis=-1) << height.unit
+        self = xyz.view(cls._location_dtype, cls).reshape(selenodetic.shape)
+        self._ellipsoid = ellipsoid
+        return self
 
     def __str__(self):
         return self.__repr__()
@@ -419,6 +422,15 @@ def selenodetic(self):
         """Convert to selenodetic coordinates."""
         return self.to_selenodetic()
 
+    @property
+    def ellipsoid(self):
+        """The default ellipsoid used to convert to selenodetic coordinates."""
+        return self._ellipsoid
+
+    @ellipsoid.setter
+    def ellipsoid(self, ellipsoid):
+        self._ellipsoid = _check_ellipsoid(ellipsoid)
+
     def to_selenodetic(self, ellipsoid=None):
         """Convert to selenodetic coordinates (lat, lon, height).
 
@@ -431,10 +443,10 @@ def to_selenodetic(self, ellipsoid=None):
             `~astropy.coordinates.Latitude`, and `~astropy.units.Quantity`
 
         """
-        ellipsoid = _check_ellipsoid(ellipsoid)
+        ellipsoid = _check_ellipsoid(ellipsoid, default=self.ellipsoid)
         xyz = self.view(self._array_dtype, u.Quantity)
         llh = CartesianRepresentation(xyz, xyz_axis=-1, copy=False).represent_as(
-            LUNAR_ELLIPSOIDS[ellipsoid],
+            SELENOIDS[ellipsoid],
         )
         return SelenodeticLocation(
             Longitude(llh.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
@@ -544,6 +556,8 @@ def __getitem__(self, item):
 
     def __array_finalize__(self, obj):
         super().__array_finalize__(obj)
+        if hasattr(obj, "_ellipsoid"):
+            self._ellipsoid = obj._ellipsoid
 
     def __len__(self):
         if self.shape == ():
diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py
index 5563028..6c2826e 100644
--- a/lunarsky/spice_utils.py
+++ b/lunarsky/spice_utils.py
@@ -12,14 +12,9 @@
 from .data import DATA_PATH
 
 _J2000 = Time("J2000")
-LUNAR_RADIUS = 1737.1e3  # m
 
 TEMPORARY_KERNEL_DIR = tempfile.TemporaryDirectory()
 
-# TODO --> Look for updated kernels and coordinate definitions.
-#          Use https://naif.jpl.nasa.gov/pub/naif/MER/misc/pck/pck00008.tpc for lunar radii
-
-
 def check_is_loaded(search):
     """
     Search the kernel pool variable names for a given string.
@@ -73,7 +68,7 @@ def furnish_kernels():
     return kernel_paths
 
 
-def lunar_surface_ephem(latitude, longitude, station_num=98):
+def lunar_surface_ephem(pos_x, pos_y, pos_z, station_num=98):
     """
     Make an SPK for the point on the lunar surface
 
@@ -81,10 +76,8 @@ def lunar_surface_ephem(latitude, longitude, station_num=98):
 
     Parameters
     ----------
-    latitude: float
-        Mean-Earth frame selenodetic latitude in degrees.
-    longitude: float
-        Mean-Earth frame selenodetic longitude in degrees.
+    pos_x, pos_y, pos_z: float
+        MCMF frame cartesian position in km
     station_num: int
         Station number
 
@@ -95,15 +88,10 @@ def lunar_surface_ephem(latitude, longitude, station_num=98):
     """
     point_id = 301000 + station_num
 
-    lat = np.radians(latitude)
-    lon = np.radians(longitude)
     ets = np.array([spice.str2et("1950-01-01"), spice.str2et("2150-01-01")])
-    pos_mcmf = spice.latrec(
-        LUNAR_RADIUS / 1e3, lon, lat
-    )  # TODO Use MoonLocation instead?
 
     states = np.zeros((len(ets), 6))
-    states[:, :3] = np.repeat(pos_mcmf[None, :], len(ets), axis=0)
+    states[:, :3] = np.repeat([[pos_x], [pos_y], [pos_z]], len(ets), axis=1).T
 
     center = 301
     frame = "MOON_ME"
diff --git a/lunarsky/tests/test_transforms.py b/lunarsky/tests/test_transforms.py
index 39befe7..82e6403 100644
--- a/lunarsky/tests/test_transforms.py
+++ b/lunarsky/tests/test_transforms.py
@@ -6,6 +6,7 @@
 from astropy.utils import IncompatibleShapeError, exceptions
 from astropy.tests.helper import assert_quantity_allclose
 import lunarsky
+from lunarsky.moon import SELENOIDS
 import pytest
 
 try:
@@ -59,10 +60,11 @@ def test_icrs_to_topo_long_time(time, lat, lon, grcat):
         (jd_4mo[:2], [10.3, 11.2], [0.0, 1.4]),
     ],
 )
-def test_transform_loops(obj, path, time, lat, lon):
+@pytest.mark.parametrize("ell", SELENOIDS)
+def test_transform_loops(obj, path, time, lat, lon, ell):
     # Tests from ICRS -> path -> ICRS
     obj = lunarsky.SkyCoord(obj)  # Ensure we're working with lunarsky-compatible SkyCoord
-    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon)
+    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon, ellipsoid=ell)
     fdict = {
         "lunartopo": lunarsky.LunarTopo(location=loc, obstime=time),
         "mcmf": lunarsky.MCMF(obstime=time),
@@ -80,9 +82,10 @@ def test_transform_loops(obj, path, time, lat, lon):
     assert_quantity_allclose(obj.cartesian.xyz, orig_pos, atol=tol)
 
 
-def test_topo_to_topo():
+@pytest.mark.parametrize("ell", SELENOIDS)
+def test_topo_to_topo(ell):
     # Check that zenith source transforms properly
-    loc0 = lunarsky.MoonLocation.from_selenodetic(lat=0, lon=90)
+    loc0 = lunarsky.MoonLocation.from_selenodetic(lat=0, lon=90, ellipsoid=ell)
     loc1 = lunarsky.MoonLocation.from_selenodetic(lat=0, lon=0)
 
     src = lunarsky.SkyCoord(alt="90d", az="0d", frame="lunartopo", location=loc0)
@@ -109,7 +112,8 @@ def test_mcmf_to_mcmf():
     # Tolerance to allow for lunar precession / nutation.
 
 
-def test_earth_from_moon():
+@pytest.mark.parametrize("ell", SELENOIDS.keys())
+def test_earth_from_moon(ell):
     # Look at the position of the Earth in lunar reference frames.
 
     # The lunar apo/perigee vary over time. These are
@@ -138,7 +142,7 @@ def test_earth_from_moon():
 
     # Now test LunarTopo frame positions
     lat, lon = 0, 0  # deg
-    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon)
+    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon, ellipsoid=ell)
     topo = mcmf.transform_to(lunarsky.LunarTopo(location=loc, obstime=jd_4mo))
 
     # The Earth should stay within 10 deg of zenith of lat=lon=0 position
@@ -273,22 +277,24 @@ def test_incompatible_transform(fromframe):
         src.transform_to(ltop)
 
 
-def test_finite_vs_spherical():
-    # Transform MCMF coordinates with distance and without
+@pytest.mark.parametrize("ell", SELENOIDS)
+def test_finite_vs_spherical(ell):
+    # Transform MCMF coordinates with distance and without units
+    # Check consistency with ellipsoid equatorial radius
     # Assumes infinite distance if no unit given, as astropy does.
-    R0 = 2e3
-    N = 100
-    phi = np.linspace(0, 2 * np.pi, N)
-    xyz = R0 * un.km * np.array([np.cos(phi), np.sin(phi), np.zeros(N)])
-    with_units = lunarsky.SkyCoord(lunarsky.MCMF(*xyz))
-    xyz = R0 * np.array([np.cos(phi), np.sin(phi), np.zeros(N)])
-    sans_units = lunarsky.SkyCoord(lunarsky.MCMF(*xyz))
-
-    loc = lunarsky.MoonLocation.from_selenodetic(lon=180 * un.deg, lat=0, height=0)
+
+    R0 = 404789     # km
+    xyz = np.array([[R0, -R0], [0, 0], [0, 0]])
+    with_units = lunarsky.SkyCoord(lunarsky.MCMF(*(xyz * un.km)))
+    sans_units = lunarsky.SkyCoord(lunarsky.MCMF(*(xyz)))
+
+    loc = lunarsky.MoonLocation.from_selenodetic(
+        lon=180 * un.deg, lat=0, height=0, ellipsoid=ell
+    )
     altaz_with_units = with_units.transform_to(lunarsky.LunarTopo(location=loc))
     with pytest.warns(exceptions.AstropyUserWarning, match="Coordinates do not "):
         altaz_sans_units = sans_units.transform_to(lunarsky.LunarTopo(location=loc))
 
-    radius = lunarsky.spice_utils.LUNAR_RADIUS * un.m
-    assert np.all(altaz_with_units.distance < radius + R0 * un.km)
+    dists = R0 * un.km + np.array([1, -1]) * SELENOIDS[ell]._equatorial_radius
+    assert np.all(altaz_with_units.distance == dists)
     assert_quantity_allclose(altaz_sans_units.distance, 1.0)
diff --git a/lunarsky/topo.py b/lunarsky/topo.py
index d45def0..da790b3 100644
--- a/lunarsky/topo.py
+++ b/lunarsky/topo.py
@@ -87,31 +87,25 @@ def zen(self):
 # Helper functions
 # -----------------
 
-
-def _spice_setup(latitude, longitude, station_id):
-
-    latlonids = np.stack(
-        [np.atleast_1d(latitude), np.atleast_1d(longitude), station_id]
-    ).T
-    if latlonids.ndim == 1:
-        latlonids = latlonids[None, :]
-
-    for lat, lon, sid in latlonids:
-        sid = int(sid)  # Station IDs must be ints, but are converted to float above.
+def _spice_setup(locations, station_ids):
+    for li, loc in enumerate(np.atleast_1d(locations)):
+        sid = int(station_ids[li])
         frameloaded = check_is_loaded(f"FRAME_LUNAR-TOPO-{sid}")
         if not frameloaded:
             lunar_surface_ephem(
-                lat, lon, station_num=sid
+                loc.x.to_value('km'),
+                loc.y.to_value('km'),
+                loc.z.to_value('km'),
+                station_num=sid
             )  # Furnishes SPK for lunar surface point
             station_name, idnum, frame_specs, latlon = topo_frame_def(
-                lat, lon, moon=True, station_num=sid
+                loc.lat.deg, loc.lon.deg, moon=True, station_num=sid
             )
             frame_strs = ["{}={}".format(k, v) for (k, v) in frame_specs.items()]
             spice.lmpool(frame_strs)
 
 
 def make_transform(coo, toframe):
-
     ap_to_spice = {"icrs": ("J2000", 0), "mcmf": ("MOON_ME", 301)}
     # Get target frame and source frame names
     if isinstance(coo, LunarTopo):
@@ -138,7 +132,8 @@ def make_transform(coo, toframe):
     shape_out = check_broadcast(coo.shape, ets.shape, location.shape)
 
     # Set up SPICE ephemerides and frame details
-    _spice_setup(location.lat.deg, location.lon.deg, stat_ids)
+    ## Use MCMF cartesian position
+    _spice_setup(location, stat_ids)
 
     ets_ids = np.atleast_2d(np.stack(np.broadcast_arrays(ets, stat_ids)).T)
 
@@ -205,13 +200,11 @@ def make_transform(coo, toframe):
 # -----------------
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, LunarTopo)
 def icrs_to_lunartopo(icrs_coo, topo_frame):
-
     return make_transform(icrs_coo, topo_frame)
 
 
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, LunarTopo, ICRS)
 def lunartopo_to_icrs(topo_coo, icrs_frame):
-
     return make_transform(topo_coo, icrs_frame)
 
 
@@ -220,13 +213,11 @@ def mcmf_to_lunartopo(mcmf_coo, topo_frame):
     # TODO:
     #   > What if mcmf_coo and topo_frame have different obstimes?
     #   > What if location has obstime?
-
     return make_transform(mcmf_coo, topo_frame)
 
 
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, LunarTopo, MCMF)
 def lunartopo_to_mcmf(topo_coo, mcmf_frame):
-
     return make_transform(topo_coo, mcmf_frame)
 
 

From b30deb766ca0b540daaa2dd86609d371448e3eb7 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Tue, 19 Mar 2024 14:54:26 -0400
Subject: [PATCH 04/10] add test of ellipsoid geometry

---
 CHANGELOG.md                      |  5 ++++
 lunarsky/tests/test_transforms.py | 48 ++++++++++++++++++++++++++++++-
 2 files changed, 52 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c3bf968..45702df 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,10 @@
 # Changelog
 
+## Unreleased
+
+## Added
+- Support for ellipsoid models for selenodetic coordinates (non-spherical)
+
 ## [0.2.0] -- 2022-10-12
 
 ## Fixed
diff --git a/lunarsky/tests/test_transforms.py b/lunarsky/tests/test_transforms.py
index 82e6403..91a35c8 100644
--- a/lunarsky/tests/test_transforms.py
+++ b/lunarsky/tests/test_transforms.py
@@ -15,11 +15,16 @@
     from astropy.coordinates.angle_utilities import angular_separation
 
 # Lunar station positions
-Nangs = 3
+Nangs = 7
 latitudes = np.linspace(-90, 90, Nangs)
 longitudes = np.linspace(0, 360, Nangs)
 latlons_grid = [(lat, lon) for lon in longitudes for lat in latitudes]
 
+# Avoiding poles:
+latlons_grid_nopole = [
+    (lat, lon) for lon in np.linspace(0.1, 360, Nangs, endpoint=False)
+    for lat in np.linspace(-70.0, 70.0, Nangs, endpoint=False)
+]
 
 # Times
 t0 = Time("2020-10-28T15:30:00")
@@ -298,3 +303,44 @@ def test_finite_vs_spherical(ell):
     dists = R0 * un.km + np.array([1, -1]) * SELENOIDS[ell]._equatorial_radius
     assert np.all(altaz_with_units.distance == dists)
     assert_quantity_allclose(altaz_sans_units.distance, 1.0)
+
+@pytest.mark.parametrize("ell", SELENOIDS)
+@pytest.mark.parametrize("lat,lon", latlons_grid_nopole)
+def test_topo_zenith_shift(ell, lat, lon):
+    # Verify that a given source at zenith in one topo frame shifts
+    # predictably when viewed from a topo frame with the same lat/lon but
+    # different ellipsoid
+
+    # Checking that the ellipsoid is interpreted correctly
+
+    # This test is a little sketchy... will need to review this later.
+    #   Some discrepancies for GRAIL23 selenoid when the source distance is large. Chooseing 1000 km for now.
+    #   Also fails near poles.
+
+    # Comparing against the SPHERE ellipsoid. Test fails for this due to divide by zero
+    if ell == "SPHERE":
+        return
+
+    lat *= un.deg
+    lon *= un.deg
+
+    loc0 = lunarsky.MoonLocation.from_selenodetic(lon, lat, ellipsoid='SPHERE')     # For reference.
+    loc1 = lunarsky.MoonLocation.from_selenodetic(lon, lat, ellipsoid=ell)       # Same lat/lon = different place for different ellipsoid
+
+    # Zenith source at finite distance over loc0.
+    src0 = lunarsky.SkyCoord(alt='90d', az='0d', distance=1e3 * un.km, frame=lunarsky.LunarTopo(location=loc0, obstime=Time.now()))
+    src1 = src0.transform_to(lunarsky.LunarTopo(location=loc1))
+
+    # Law of sines:
+    #   Geometry among zenith angle, mcmf station vector,
+    #   and selenocentric vs. selenodetic latitudes
+
+    R1 = loc1.mcmf.cartesian.norm()
+    lat_cen = loc1.mcmf.spherical.lat
+    lat_det = loc1.lat
+    diff = (src1.distance / np.sin((np.abs(lat_det - lat_cen)).rad)) - (R1 / np.sin(src1.zen.rad))
+    assert_quantity_allclose(
+        src1.distance / np.sin((np.abs(lat_det - lat_cen)).rad),
+        R1 / np.sin(src1.zen.rad),
+        atol=1*un.km,
+    )

From 40ad27f0fa831c73575b172ac738cacaac4835e6 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Tue, 19 Mar 2024 14:55:27 -0400
Subject: [PATCH 05/10] linting and tests

---
 CHANGELOG.md                      |  2 +-
 lunarsky/moon.py                  | 33 +++++++++++++++++++------------
 lunarsky/spice_utils.py           |  1 +
 lunarsky/tests/test_moon.py       | 20 +++++++++++++++----
 lunarsky/tests/test_spice.py      |  6 +++---
 lunarsky/tests/test_transforms.py | 31 ++++++++++++++++++-----------
 lunarsky/topo.py                  | 10 +++++-----
 7 files changed, 66 insertions(+), 37 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 45702df..7d85795 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,6 @@
 # Changelog
 
-## Unreleased
+## [0.2.2] -- 2022-03-19
 
 ## Added
 - Support for ellipsoid models for selenodetic coordinates (non-spherical)
diff --git a/lunarsky/moon.py b/lunarsky/moon.py
index 0a93833..3590e11 100644
--- a/lunarsky/moon.py
+++ b/lunarsky/moon.py
@@ -7,7 +7,6 @@
 from astropy.coordinates.representation.geodetic import BaseGeodeticRepresentation
 from astropy.coordinates.representation import (
     CartesianRepresentation,
-    SphericalRepresentation,
 )
 from astropy.coordinates.attributes import Attribute
 
@@ -17,30 +16,37 @@
 
 __all__ = ["MoonLocation", "MoonLocationAttribute"]
 
+
 class SPHERESelenodeticRepresentation(BaseGeodeticRepresentation):
     """Lunar ellipsoid as a sphere
 
-        Radius defined by lunarsky.spice_utils.LUNAR_RADIUS
+    Radius defined by lunarsky.spice_utils.LUNAR_RADIUS
     """
-    _equatorial_radius = LUNAR_RADIUS * u.m 
+
+    _equatorial_radius = LUNAR_RADIUS * u.m
     _flattening = 0.0
 
+
 class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation):
     """Lunar ellipsoid from NASA/GSFC "Planetary Fact Sheet"
 
-       https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
+    https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
     """
+
     _equatorial_radius = 1738.1e3 * u.m
     _flattening = 0.0012
 
+
 class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation):
     """Lunar ellipsoid defined by gravimetry of GRAIL data.
 
     https://doi.org/10.1007/s40328-023-00415-w
     """
+
     _equatorial_radius = 1737576.6 * u.m
     _flattening = 0.000305
 
+
 class CE1LAM10SelenodeticRepresentation(BaseGeodeticRepresentation):
     """Lunar ellipsoid from Chang'e 1 laser altimetry.
 
@@ -48,16 +54,17 @@ class CE1LAM10SelenodeticRepresentation(BaseGeodeticRepresentation):
 
     https://doi.org/10.1007/s11430-010-4060-6
     """
+
     _equatorial_radius = 1737.632 * u.km
-    _flattening = 1/973.463
+    _flattening = 1 / 973.463
 
 
 # Define reference ellipsoids
 SELENOIDS = {
-    "SPHERE" : SPHERESelenodeticRepresentation,
-    "GSFC" : GSFCSelenodeticRepresentation,
-    "GRAIL23" : GRAIL23SelenodeticRepresentation,
-    "CE-1-LAM-GEO" : CE1LAM10SelenodeticRepresentation,
+    "SPHERE": SPHERESelenodeticRepresentation,
+    "GSFC": GSFCSelenodeticRepresentation,
+    "GRAIL23": GRAIL23SelenodeticRepresentation,
+    "CE-1-LAM-GEO": CE1LAM10SelenodeticRepresentation,
 }
 
 
@@ -242,7 +249,9 @@ def _set_site_id(cls, inst):
             raise ValueError("Too many unique MoonLocation objects open at once.")
 
         for llh in llh_arr:
-            lonlatheight = "_".join(["{:.4f}".format(ll) for ll in llh] + [inst._ellipsoid])
+            lonlatheight = "_".join(
+                ["{:.4f}".format(ll) for ll in llh] + [inst._ellipsoid]
+            )
             if lonlatheight not in cls._existing_locs:
                 new_stat_id = cls._avail_stat_ids.pop()
                 cls._existing_locs.append(lonlatheight)
@@ -451,9 +460,7 @@ def to_selenodetic(self, ellipsoid=None):
         return SelenodeticLocation(
             Longitude(llh.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
             Latitude(llh.lat, u.degree, copy=False),
-            u.Quantity(
-                llh.height, self.unit, copy=False
-            ),
+            u.Quantity(llh.height, self.unit, copy=False),
         )
 
     @property
diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py
index 6c2826e..06611ce 100644
--- a/lunarsky/spice_utils.py
+++ b/lunarsky/spice_utils.py
@@ -15,6 +15,7 @@
 
 TEMPORARY_KERNEL_DIR = tempfile.TemporaryDirectory()
 
+
 def check_is_loaded(search):
     """
     Search the kernel pool variable names for a given string.
diff --git a/lunarsky/tests/test_moon.py b/lunarsky/tests/test_moon.py
index 71bf1b9..bb5dc12 100644
--- a/lunarsky/tests/test_moon.py
+++ b/lunarsky/tests/test_moon.py
@@ -3,7 +3,9 @@
 import numpy as np
 import astropy.units as unit
 from astropy.coordinates import Longitude, Latitude
+from astropy.tests.helper import quantity_allclose
 from lunarsky import MoonLocation, MoonLocationAttribute, MCMF
+from lunarsky.moon import SELENOIDS
 
 
 class TestsWithObject:
@@ -27,7 +29,12 @@ def test_input(self):
         cartesian = MoonLocation(self.x.value, self.y.value, self.z.value, self.x.unit)
         assert np.all(cartesian == self.location)
         spherical = MoonLocation(self.lon.deg, self.lat.deg, self.h.to(unit.m))
-        assert np.all(spherical == self.location)
+        # TODO this was just np.all(... == ...). Seeing floating point diff since v0.2.2
+        assert quantity_allclose(
+            spherical.mcmf.cartesian.xyz,
+            self.location.mcmf.cartesian.xyz,
+            atol=1e-9 * unit.m,
+        )
 
     def test_invalid(self):
         # incomprehensible by either raises TypeError
@@ -130,7 +137,8 @@ def test_moonlocation_delete():
     assert check1 == before
 
 
-def test_station_ids():
+@pytest.mark.parametrize("ell", SELENOIDS)
+def test_station_ids(ell):
     # Check that when a MoonLocations are made, the appropriate station_ids are assigned.
 
     # Get whatever are already recorded in the class.
@@ -186,7 +194,11 @@ def test_station_ids():
 
     for gp in lonlatheights:
         lons, lats, heights = np.array(gp).T
-        locs.append(MoonLocation.from_selenodetic(lat=lats, lon=lons, height=heights))
+        locs.append(
+            MoonLocation.from_selenodetic(
+                lat=lats, lon=lons, height=heights, ellipsoid=ell
+            )
+        )
         locs[-1]._set_station_id()
 
     # Check that only unique positions got added
@@ -209,7 +221,7 @@ def test_station_ids():
         else:
             llh_arr = zip(inst.lon.deg, inst.lat.deg, inst.height.to_value("km"))
         for llh in llh_arr:
-            llhs.append("_".join(["{:.4f}".format(ll) for ll in llh]))
+            llhs.append("_".join(["{:.4f}".format(ll) for ll in llh] + [ell]))
         locstrs.append(llhs)
 
     # Check that saved location strings correspond with station IDs in each instance
diff --git a/lunarsky/tests/test_spice.py b/lunarsky/tests/test_spice.py
index 9cdf724..fcdfc51 100644
--- a/lunarsky/tests/test_spice.py
+++ b/lunarsky/tests/test_spice.py
@@ -108,12 +108,12 @@ def test_topo_kernel_setup():
 
     for filepath in lunarsky.spice_utils.KERNEL_PATHS:
         spice.furnsh(filepath)
-    lat, lon = 30, 20
+    loc = lunarsky.MoonLocation.from_selenodetic(lat="30d", lon="20d")
     station_name, idnum, frame_specs, latlon = lunarsky.spice_utils.topo_frame_def(
-        lat, lon, moon=True
+        loc.lat.deg, loc.lon.deg, moon=True
     )
     statnum = idnum - 1301000
-    lunarsky.topo._spice_setup(lat, lon, [statnum])
+    lunarsky.topo._spice_setup(loc, [statnum])
     try:
         assert lunarsky.spice_utils.check_is_loaded("*{}*".format(idnum))
     finally:
diff --git a/lunarsky/tests/test_transforms.py b/lunarsky/tests/test_transforms.py
index 91a35c8..9afa83d 100644
--- a/lunarsky/tests/test_transforms.py
+++ b/lunarsky/tests/test_transforms.py
@@ -22,7 +22,8 @@
 
 # Avoiding poles:
 latlons_grid_nopole = [
-    (lat, lon) for lon in np.linspace(0.1, 360, Nangs, endpoint=False)
+    (lat, lon)
+    for lon in np.linspace(0.1, 360, Nangs, endpoint=False)
     for lat in np.linspace(-70.0, 70.0, Nangs, endpoint=False)
 ]
 
@@ -110,9 +111,7 @@ def test_mcmf_to_mcmf():
     )
     sph0 = src.spherical
     sph1 = orig_pos.spherical
-    res = angular_separation(
-        sph0.lon, sph0.lat, sph1.lon, sph1.lat
-    ).to("deg")
+    res = angular_separation(sph0.lon, sph0.lat, sph1.lon, sph1.lat).to("deg")
     assert_quantity_allclose(res, 177 * un.deg, atol=5 * un.deg)
     # Tolerance to allow for lunar precession / nutation.
 
@@ -288,7 +287,7 @@ def test_finite_vs_spherical(ell):
     # Check consistency with ellipsoid equatorial radius
     # Assumes infinite distance if no unit given, as astropy does.
 
-    R0 = 404789     # km
+    R0 = 404789  # km
     xyz = np.array([[R0, -R0], [0, 0], [0, 0]])
     with_units = lunarsky.SkyCoord(lunarsky.MCMF(*(xyz * un.km)))
     sans_units = lunarsky.SkyCoord(lunarsky.MCMF(*(xyz)))
@@ -304,6 +303,7 @@ def test_finite_vs_spherical(ell):
     assert np.all(altaz_with_units.distance == dists)
     assert_quantity_allclose(altaz_sans_units.distance, 1.0)
 
+
 @pytest.mark.parametrize("ell", SELENOIDS)
 @pytest.mark.parametrize("lat,lon", latlons_grid_nopole)
 def test_topo_zenith_shift(ell, lat, lon):
@@ -314,7 +314,8 @@ def test_topo_zenith_shift(ell, lat, lon):
     # Checking that the ellipsoid is interpreted correctly
 
     # This test is a little sketchy... will need to review this later.
-    #   Some discrepancies for GRAIL23 selenoid when the source distance is large. Chooseing 1000 km for now.
+    #   Some discrepancies for GRAIL23 selenoid when the source distance is large.
+    #       Choosing 1000 km for now.
     #   Also fails near poles.
 
     # Comparing against the SPHERE ellipsoid. Test fails for this due to divide by zero
@@ -324,11 +325,20 @@ def test_topo_zenith_shift(ell, lat, lon):
     lat *= un.deg
     lon *= un.deg
 
-    loc0 = lunarsky.MoonLocation.from_selenodetic(lon, lat, ellipsoid='SPHERE')     # For reference.
-    loc1 = lunarsky.MoonLocation.from_selenodetic(lon, lat, ellipsoid=ell)       # Same lat/lon = different place for different ellipsoid
+    loc0 = lunarsky.MoonLocation.from_selenodetic(
+        lon, lat, ellipsoid="SPHERE"
+    )  # For reference.
+    loc1 = lunarsky.MoonLocation.from_selenodetic(
+        lon, lat, ellipsoid=ell
+    )  # Same lat/lon = different place for different ellipsoid
 
     # Zenith source at finite distance over loc0.
-    src0 = lunarsky.SkyCoord(alt='90d', az='0d', distance=1e3 * un.km, frame=lunarsky.LunarTopo(location=loc0, obstime=Time.now()))
+    src0 = lunarsky.SkyCoord(
+        alt="90d",
+        az="0d",
+        distance=1e3 * un.km,
+        frame=lunarsky.LunarTopo(location=loc0, obstime=Time.now()),
+    )
     src1 = src0.transform_to(lunarsky.LunarTopo(location=loc1))
 
     # Law of sines:
@@ -338,9 +348,8 @@ def test_topo_zenith_shift(ell, lat, lon):
     R1 = loc1.mcmf.cartesian.norm()
     lat_cen = loc1.mcmf.spherical.lat
     lat_det = loc1.lat
-    diff = (src1.distance / np.sin((np.abs(lat_det - lat_cen)).rad)) - (R1 / np.sin(src1.zen.rad))
     assert_quantity_allclose(
         src1.distance / np.sin((np.abs(lat_det - lat_cen)).rad),
         R1 / np.sin(src1.zen.rad),
-        atol=1*un.km,
+        atol=1 * un.km,
     )
diff --git a/lunarsky/topo.py b/lunarsky/topo.py
index da790b3..bc39ff4 100644
--- a/lunarsky/topo.py
+++ b/lunarsky/topo.py
@@ -87,16 +87,17 @@ def zen(self):
 # Helper functions
 # -----------------
 
+
 def _spice_setup(locations, station_ids):
     for li, loc in enumerate(np.atleast_1d(locations)):
         sid = int(station_ids[li])
         frameloaded = check_is_loaded(f"FRAME_LUNAR-TOPO-{sid}")
         if not frameloaded:
             lunar_surface_ephem(
-                loc.x.to_value('km'),
-                loc.y.to_value('km'),
-                loc.z.to_value('km'),
-                station_num=sid
+                loc.x.to_value("km"),
+                loc.y.to_value("km"),
+                loc.z.to_value("km"),
+                station_num=sid,
             )  # Furnishes SPK for lunar surface point
             station_name, idnum, frame_specs, latlon = topo_frame_def(
                 loc.lat.deg, loc.lon.deg, moon=True, station_num=sid
@@ -132,7 +133,6 @@ def make_transform(coo, toframe):
     shape_out = check_broadcast(coo.shape, ets.shape, location.shape)
 
     # Set up SPICE ephemerides and frame details
-    ## Use MCMF cartesian position
     _spice_setup(location, stat_ids)
 
     ets_ids = np.atleast_2d(np.stack(np.broadcast_arrays(ets, stat_ids)).T)

From 300c502d83fb2b2286bdf5d7f6a43fc1e6a60a46 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Tue, 19 Mar 2024 18:44:36 -0400
Subject: [PATCH 06/10] run CI tests on later python versions

---
 .github/workflows/testsuite.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testsuite.yaml b/.github/workflows/testsuite.yaml
index 1452365..af58e42 100644
--- a/.github/workflows/testsuite.yaml
+++ b/.github/workflows/testsuite.yaml
@@ -63,7 +63,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [3.7, 3.8, 3.9]
+        python-version: [3.7, 3.8, 3.9, 3.10, 3.11]
     steps:
       - uses: actions/checkout@v3
         with:

From e41afe46d4f404f2f20ff91ec7fe1fa93285c7f8 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Fri, 22 Mar 2024 14:20:53 -0400
Subject: [PATCH 07/10] deprecate py38

---
 .github/workflows/testsuite.yaml |  8 ++++----
 .pre-commit-config.yaml          |  2 +-
 CHANGELOG.md                     |  3 +++
 lunarsky/moon.py                 | 10 ----------
 setup.py                         |  2 +-
 5 files changed, 9 insertions(+), 16 deletions(-)

diff --git a/.github/workflows/testsuite.yaml b/.github/workflows/testsuite.yaml
index af58e42..a2a3654 100644
--- a/.github/workflows/testsuite.yaml
+++ b/.github/workflows/testsuite.yaml
@@ -15,10 +15,10 @@ jobs:
       - uses: actions/checkout@v3
         with:
          fetch-depth: 1
-      - name: Set up Python 3.8
+      - name: Set up Python 3.9
         uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: 3.9
       - name: Setup Environment
         run: |
           pip install pre-commit
@@ -36,7 +36,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: 3.9
       - name: Get current date
         id: date
         run: echo "::set-output name=date::$(date +'%Y-%m-%d')"
@@ -63,7 +63,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [3.7, 3.8, 3.9, 3.10, 3.11]
+        python-version: ["3.9", "3.10", "3.11"]
     steps:
       - uses: actions/checkout@v3
         with:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 8955147..2866e1b 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -26,6 +26,6 @@ repos:
       rev: 22.3.0
       hooks:
           - id: black
-            language_version: python3.8
+            language_version: python3.9
             args:
               - --line-length=90
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7d85795..75218fd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,9 @@
 ## Added
 - Support for ellipsoid models for selenodetic coordinates (non-spherical)
 
+## Deprecated
+- Removed support for Python < 3.9. Astropy >= 6.0 is required for ellipsoid support
+
 ## [0.2.0] -- 2022-10-12
 
 ## Fixed
diff --git a/lunarsky/moon.py b/lunarsky/moon.py
index 3590e11..8c9cdf0 100644
--- a/lunarsky/moon.py
+++ b/lunarsky/moon.py
@@ -81,16 +81,6 @@ def _check_ellipsoid(ellipsoid=None, default="SPHERE"):
     return ellipsoid
 
 
-def _add_rm_ellipsoid(ellipsoid, name, add=True):
-    """Add or remove an ellipsoid class from SELENOIDS.
-    Needed for testing."""
-    global SELENOIDS
-    if add:
-        SELENOIDS[name] = ellipsoid
-    else:
-        SELENOIDS.pop(name, None)
-
-
 class MoonLocationInfo(QuantityInfoBase):
     """
     Container for meta information like name, description, format.  This is
diff --git a/setup.py b/setup.py
index 4fc2b8f..60461b5 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,7 @@
     "test_suite": "pytest",
     "tests_require": ["pytest"],
     "setup_requires": ["pytest-runner", "setuptools_scm"],
-    "install_requires": ["numpy>=1.15", "astropy>3.0", "spiceypy", "jplephem"],
+    "install_requires": ["numpy>=1.15", "astropy>=6.0.0", "spiceypy", "jplephem"],
     "classifiers": [
         "Development Status :: 3 - Alpha",
         "Intended Audience :: Science/Research",

From 0ad61ac1974f4175c972eb733c376231a84c9e5a Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Fri, 22 Mar 2024 14:54:51 -0400
Subject: [PATCH 08/10] test coverage

---
 lunarsky/moon.py            |  2 +-
 lunarsky/tests/test_moon.py | 11 ++++++++---
 2 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/lunarsky/moon.py b/lunarsky/moon.py
index 8c9cdf0..f3df176 100644
--- a/lunarsky/moon.py
+++ b/lunarsky/moon.py
@@ -379,7 +379,7 @@ def from_selenodetic(cls, lon, lat, height=0.0, ellipsoid=None):
         selenodetic = SELENOIDS[ellipsoid](lon, lat, height, copy=False)
         xyz = selenodetic.to_cartesian().get_xyz(xyz_axis=-1) << height.unit
         self = xyz.view(cls._location_dtype, cls).reshape(selenodetic.shape)
-        self._ellipsoid = ellipsoid
+        self.ellipsoid = ellipsoid
         return self
 
     def __str__(self):
diff --git a/lunarsky/tests/test_moon.py b/lunarsky/tests/test_moon.py
index bb5dc12..dd267e1 100644
--- a/lunarsky/tests/test_moon.py
+++ b/lunarsky/tests/test_moon.py
@@ -2,10 +2,11 @@
 import copy
 import numpy as np
 import astropy.units as unit
-from astropy.coordinates import Longitude, Latitude
+from astropy.coordinates import Longitude, Latitude, SphericalRepresentation, SphericalDifferential
 from astropy.tests.helper import quantity_allclose
+from astropy.coordinates.tests.test_representation import representation_equal
 from lunarsky import MoonLocation, MoonLocationAttribute, MCMF
-from lunarsky.moon import SELENOIDS
+from lunarsky.moon import SELENOIDS, MoonLocationInfo
 
 
 class TestsWithObject:
@@ -21,7 +22,7 @@ def setup_method(self):
         self.lat = Latitude([+0.0, 30.0, 60.0, +90.0, -90.0, -60.0, -30.0, 0.0], unit.deg)
         self.h = unit.Quantity([0.1, 0.5, 1.0, -0.5, -1.0, +4.2, -11.0, -0.1], unit.m)
         self.location = MoonLocation.from_selenodetic(self.lon, self.lat, self.h)
-        self.x, self.y, self.z = self.location.to_selenocentric()
+        self.x, self.y, self.z = self.location.selenocentric
 
     def test_input(self):
         cartesian = MoonLocation(self.x, self.y, self.z)
@@ -68,6 +69,10 @@ def test_invalid(self):
         with pytest.raises(ValueError):
             MoonLocation.from_selenodetic(self.lon, self.lat, self.h[:5])
 
+        # invalid ellipsoid
+        with pytest.raises(ValueError):
+            loc = MoonLocation.from_selenodetic(self.lon, self.lat, self.h, ellipsoid="INVALID")
+
     def test_attributes(self):
         assert np.allclose(self.location.height, self.h)
         assert np.allclose(self.location.lon, self.lon)

From e05c1e4772f13bff61c0d87ff25de5bd3ca3ef73 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Fri, 22 Mar 2024 14:56:59 -0400
Subject: [PATCH 09/10] linting

---
 lunarsky/tests/test_moon.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/lunarsky/tests/test_moon.py b/lunarsky/tests/test_moon.py
index dd267e1..8decb46 100644
--- a/lunarsky/tests/test_moon.py
+++ b/lunarsky/tests/test_moon.py
@@ -2,11 +2,13 @@
 import copy
 import numpy as np
 import astropy.units as unit
-from astropy.coordinates import Longitude, Latitude, SphericalRepresentation, SphericalDifferential
+from astropy.coordinates import (
+    Longitude,
+    Latitude,
+)
 from astropy.tests.helper import quantity_allclose
-from astropy.coordinates.tests.test_representation import representation_equal
 from lunarsky import MoonLocation, MoonLocationAttribute, MCMF
-from lunarsky.moon import SELENOIDS, MoonLocationInfo
+from lunarsky.moon import SELENOIDS
 
 
 class TestsWithObject:
@@ -71,7 +73,7 @@ def test_invalid(self):
 
         # invalid ellipsoid
         with pytest.raises(ValueError):
-            loc = MoonLocation.from_selenodetic(self.lon, self.lat, self.h, ellipsoid="INVALID")
+            MoonLocation.from_selenodetic(self.lon, self.lat, self.h, ellipsoid="INVALID")
 
     def test_attributes(self):
         assert np.allclose(self.location.height, self.h)

From 172a77a6833ec0d4dc5fdfbb3b7613bfe374d125 Mon Sep 17 00:00:00 2001
From: aelanman <adamelanman@pm.me>
Date: Fri, 22 Mar 2024 15:46:21 -0400
Subject: [PATCH 10/10] test in python 3.12 too

---
 .github/workflows/testsuite.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testsuite.yaml b/.github/workflows/testsuite.yaml
index a2a3654..cc6ffb3 100644
--- a/.github/workflows/testsuite.yaml
+++ b/.github/workflows/testsuite.yaml
@@ -63,7 +63,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: ["3.9", "3.10", "3.11"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
     steps:
       - uses: actions/checkout@v3
         with: