Add new mean_equator_and_equinox_of_date frame

Inspired by discussion #1034.

CHANGELOG.rst

@@ -12,6 +12,10 @@ Released versions
 v1.50 — Not yet released
 ------------------------
 
+* A new :class:`~skyfield.framelib.mean_equator_and_equinox_of_date`
+  coordinate frame lets users generate the same coordinates that an
+  almanac might give.
+
 * Skyfield now offers a Solar System Barycenter object, so users don’t
   have to construct the position themselves: ``SSB.at(t)`` returns a
   position whose coordinates and velocity are both zero in the ICRS.

documentation/accuracy-efficiency.rst

@@ -23,7 +23,7 @@ This is used in:
 * Geographic positions generated by an Earth ellipsoid
   like the :data:`~skyfield.toposlib.wgs84` model.
 
-* The :data:`~skyfield.framelib.true_equator_and_equinox_of_date`
+* The :class:`~skyfield.framelib.true_equator_and_equinox_of_date`
   reference frame that can be used to generate |xyz| coordinates
   as described in the :doc:`coordinates` chapter.
 

documentation/api-framelib.rst

@@ -9,7 +9,9 @@ whose use is described in `reference_frames`.
 
 .. currentmodule:: skyfield.framelib
 
-.. autodata:: true_equator_and_equinox_of_date
+.. autoclass:: mean_equator_and_equinox_of_date
+
+.. autoclass:: true_equator_and_equinox_of_date
 
 .. autodata:: itrs
 

documentation/examples.rst

@@ -549,7 +549,7 @@ as the target moves across the sky:
    Or, if you instead want to know how fast the target is moving
    against the background of stars,
    you can pass Skyfield’s built-in
-   :data:`~skyfield.framelib.true_equator_and_equinox_of_date` reference frame
+   :class:`~skyfield.framelib.true_equator_and_equinox_of_date` reference frame
    to compute rates of moment in right ascension and declination:
 
    .. testcode::

skyfield/framelib.py

@@ -49,11 +49,18 @@ class ICRS(object):
     def rotation_at(t):
         return _identity
 
-def build_ecliptic_matrix(t):
-    # Build the matrix to rotate an ICRF vector into ecliptic coordinates.
-    _, d_eps = t._nutation_angles_radians
-    true_obliquity = t._mean_obliquity_radians + d_eps
-    return mxm(rot_x(- true_obliquity), t.M)
+class mean_equator_and_equinox_of_date(object):
+    """The coordinate frame of Earth’s mean equator and equinox.
+
+    This frame is used for measuring right ascension and declination.
+    It tracks the Earth’s ‘mean’ equator and equinox which shift slowly
+    across the sky due to precession, but ignores the smaller effects of
+    nutation.
+
+    """
+    @staticmethod
+    def rotation_at(t):
+        return t.P
 
 class true_equator_and_equinox_of_date(object):
     """The dynamical frame of Earth’s true equator and true equinox of date.
@@ -80,8 +87,6 @@ class true_equator_and_equinox_of_date(object):
     def rotation_at(t):
         return t.M
 
-true_equator_and_equinox_of_date = true_equator_and_equinox_of_date()
-
 _itrs_angvel_matrix = array((
     (0.0, DAY_S * ANGVEL, 0.0),
     (-DAY_S * ANGVEL, 0.0, 0.0),
@@ -140,6 +145,12 @@ def _dRdt_times_RT_at(t):
 
 itrs = itrs()
 
+def build_ecliptic_matrix(t):
+    # Build the matrix to rotate an ICRF vector into ecliptic coordinates.
+    _, d_eps = t._nutation_angles_radians
+    true_obliquity = t._mean_obliquity_radians + d_eps
+    return mxm(rot_x(- true_obliquity), t.M)
+
 class ecliptic_frame(object):
     """Reference frame of the true ecliptic and equinox of date."""
     @staticmethod

skyfield/tests/test_frames.py

@@ -1,6 +1,6 @@
 from numpy import cos
 from skyfield import framelib
-from skyfield.api import Topos, load, wgs84
+from skyfield.api import SSB, Star, Topos, load, wgs84
 from skyfield.constants import AU_M, ERAD
 from skyfield.positionlib import Geocentric
 
@@ -86,6 +86,23 @@ def test_frame_without_spin():
     r, v = g.frame_xyz_and_velocity(f)
     Geocentric.from_time_and_frame_vectors(t, f, r, v)
 
+def test_true_equator_and_equinox_of_date():
+    ts = load.timescale()
+    t = ts.utc(2025, 7, 2.375)
+    alpha_andromeda = Star(
+        ra_hours=0.13976888866666667, dec_degrees=29.09082805,
+        ra_mas_per_year=135.68, dec_mas_per_year=-162.95, parallax_mas=33.6,
+        epoch=ts.J(1991.25),
+    )
+    dec, ra, _ = SSB.at(t).observe(alpha_andromeda).frame_latlon(
+        framelib.mean_equator_and_equinox_of_date
+    )
+    ra.preference = 'hours'
+
+    # Position from page H2 of the 2025 Astronomical Almanac:
+    assert ra.hstr(places=1) == '00h 09m 42.7s'
+    assert dec.dstr() == '+29deg 13\' 52.0"'
+
 def test_tirs_at_least_runs():
     # TODO: find an external source for a TIRS vector to test against.
     # For now, just make sure it doesn't raise an exception.