For #510, add and test hadec() position method

CHANGELOG.rst

@@ -5,6 +5,9 @@ Changelog
 v1.36 — ?
 ---------
 
+* Added :meth:`~skyfield.positionlib.ICRF.hadec()` position method that
+  returns hour angle and declination.
+
 * The default ``str()`` and ``repr()`` strings
   for geographic positions have been streamlined,
   and no longer raise ``ValueError`` when elevation is an array.

skyfield/documentation/api.rst

@@ -332,6 +332,7 @@ All positions support a basic set of methods:
    ICRF.distance
    ICRF.speed
    ICRF.radec
+   ICRF.hadec
    ICRF.separation_from
    ICRF.frame_xyz
    ICRF.frame_xyz_and_velocity

skyfield/positionlib.py

@@ -3,7 +3,7 @@
 
 from numpy import array, einsum, full, reshape, nan, nan_to_num
 from . import framelib
-from .constants import ANGVEL, AU_M, C, ERAD, DAY_S, RAD2DEG, tau
+from .constants import ANGVEL, AU_M, C, ERAD, DAY_S, RAD2DEG, pi, tau
 from .data.spice import inertial_frames
 from .descriptorlib import reify
 from .earthlib import compute_limb_angle
@@ -270,6 +270,34 @@ def radec(self, epoch=None):
                 Angle(radians=dec, signed=True),
                 Distance(r_au))
 
+    def hadec(self):
+        """Compute hour angle, declination, and distance.
+
+        Returns a tuple of two :class:`~skyfield.units.Angle` objects
+        plus the :class:`~skyfield.units.Distance` to the target.  The
+        angles are the hour angle (±12 hours) east or west of your
+        meridian along the ITRS celestial equator, and the declination
+        (±90 degees) above or below it.  This only works for positions
+        whose center is a geographic location; otherwise, there is no
+        local meridian from which to measure the hour angle.
+
+        Because this declination is measured from the plane of the
+        Earth’s physical geographic equator, it will be slightly
+        different than the declination returned by ``radec()`` if you
+        have loaded a :ref:`polar motion` file.
+
+        """
+        R = framelib.itrs.rotation_at(self.t)
+        r = mxv(R, self.position.au)
+        au, dec, sublongtiude = to_spherical(r)
+        ha = self.center.longitude.radians - sublongtiude
+        ha += pi
+        ha %= tau
+        ha -= pi
+        return (Angle(radians=ha, preference='hours', signed=True),
+                Angle(radians=dec, signed=True),
+                Distance(au))
+
     def separation_from(self, another_icrf):
         """Return the angle between this position and another.
 

skyfield/tests/test_positions.py

@@ -182,6 +182,29 @@ def test_velocity_in_ITRF_to_GCRS2():
     acceptable_error = 4e-12
     assert relative_error < acceptable_error
 
+def test_hadec():
+    # If the DE430 ephemeris excerpt is avaiable, this test can run
+    # locally against the HA number from first line of
+    # `moon_topo_4_6_2017_mkb_sf_v5_hadec.csv.txt` at:
+    # https://github.com/skyfielders/python-skyfield/issues/510
+    #planets = api.load('de430_1850-2150.bsp')
+    #expected_ha = -0.660078756021
+
+    # But in CI, we use DE421 for space and speed.
+    planets = api.load('de421.bsp')
+    expected_ha = -0.660078752
+
+    ts = api.load.timescale()
+    ts.polar_motion_table = [0.0], [0.009587], [0.384548]
+    t = ts.utc(2017, 4, 6)
+    topos = api.wgs84.latlon(-22.959748, -67.787260, elevation_m=5186.0)
+    earth = planets['Earth']
+    moon = planets['Moon']
+    a = (earth + topos).at(t).observe(moon).apparent()
+    ha, dec, distance = a.hadec()
+    difference_mas = (ha.hours - expected_ha) * 15 * 3600 * 1e3
+    assert abs(difference_mas) < 0.03
+
 # Test that the CIRS coordinate of the TIO is consistent with the Earth Rotation Angle
 # This is mostly an internal consistency check
 def test_cirs_era():