Discover how to make ITRF->GCRS very high accuracy

But don’t turn it on by default because, wow, it’s expensive!

skyfield/positionlib.py

@@ -747,17 +747,29 @@ def ITRF_to_GCRS(t, rITRF):
     position = mxv(spin, array(rITRF))
     return mxv(t.MT, position)
 
-def ITRF_to_GCRS2(t, rITRF, vITRF):
+def ITRF_to_GCRS2(t, rITRF, vITRF, _high_accuracy=False):
     # TODO: wobble
 
     spin = rot_z(t.gast / 24.0 * tau)
     position = mxv(spin, array(rITRF))
     velocity = mxv(spin, array(vITRF))
 
-    # TODO: Would it increase accuracy to use the actual rate of spin
-    # for this date, instead of the average ANGVEL?
-    velocity[0] += DAY_S * ANGVEL * - position[1]
-    velocity[1] += DAY_S * ANGVEL * position[0]
+    # TODO: This is expensive, and should be extensively trimmed to only
+    # include the most important terms underlying GAST.  But it improves
+    # the precision by something like 1e5 times when compared to using
+    # the round number skyfield.constants.ANGVEL!
+    #
+    # See the test `test_velocity_in_ITRF_to_GCRS2()`.
+    #
+    if _high_accuracy:
+        _one_second = 1.0 / DAY_S
+        t_later = t.ts.tt_jd(t.whole, t.tt_fraction + _one_second)
+        angvel = (t_later.gast - t.gast) / 24.0 * tau
+    else:
+        angvel = ANGVEL
+
+    velocity[0] += DAY_S * angvel * - position[1]
+    velocity[1] += DAY_S * angvel * position[0]
 
     position = mxv(t.MT, position)
     velocity = mxv(t.MT, velocity)

skyfield/tests/test_positions.py

@@ -1,8 +1,8 @@
 import numpy as np
 from skyfield import api
-from skyfield.constants import DAY_S
+from skyfield.constants import DAY_S, tau
 from skyfield.earthlib import earth_rotation_angle
-from skyfield.functions import length_of
+from skyfield.functions import length_of, mxv, rot_z
 from skyfield.positionlib import ICRF, ITRF_to_GCRS2, _GIGAPARSEC_AU
 from skyfield.starlib import Star
 from .fixes import low_precision_ERA
@@ -92,16 +92,41 @@ def test_position_from_radec():
     assert length_of(p.position.au - [0, 1, 0]) < 1e-16
 
 def test_velocity_in_ITRF_to_GCRS2():
+    # TODO: Get test working with these vectors too, showing it works
+    # with a non-zero velocity vector, but in that case the test will
+    # have to be fancier in how it corrects.
+    # r = np.array([(1, 0, 0), (1, 1 / DAY_S, 0)]).T
+    # v = np.array([(0, 1, 0), (0, 1, 0)]).T
+
     ts = api.load.timescale(builtin=True)
     t = ts.utc(2020, 7, 17, 8, 51, [0, 1])
-    r = np.array([(1, 0, 0), (1, 1.0 / DAY_S, 0)]).T
-    v = np.array([(0, 1, 0), (0, 1, 0)]).T
-    r, v = ITRF_to_GCRS2(t, r, v)
+    r = np.array([(1, 0, 0), (1, 0, 0)]).T
+    v = np.array([(0, 0, 0), (0, 0, 0)]).T
+
+    r, v = ITRF_to_GCRS2(t, r, v, True)
+
+    # Rotate back to equinox-of-date before applying correction.
+    r = mxv(t.M, r)
+    v = mxv(t.M, v)
+
     r0, r1 = r.T
     v0 = v[:,0]
-    actual_velocity = r1 - r0
-    stated_velocity = v0 / DAY_S
-    assert length_of(actual_velocity - stated_velocity) < 4e-9
+
+    # Apply a correction: the instantaneous velocity does not in fact
+    # carry the position in a straight line, but in an arc around the
+    # origin; so use trigonometry to move the destination point to where
+    # linear motion would have carried it.
+    angvel = (t.gast[1] - t.gast[0]) / 24.0 * tau
+    r1 = mxv(rot_z(np.arctan(angvel) - angvel), r1)
+    r1 *= np.sqrt(1 + angvel*angvel)
+
+    actual_motion = r1 - r0
+    predicted_motion = v0 / DAY_S
+
+    relative_error = (length_of(actual_motion - predicted_motion)
+                      / length_of(actual_motion))
+
+    assert relative_error < 2e-12
 
 # Test that the CIRS coordinate of the TIO is consistent with the Earth Rotation Angle
 # This is mostly an internal consistency check

skyfield/timelib.py

@@ -711,6 +711,7 @@ class Time(object):
         # TODO: move this into an eqeq function?
         c_terms = equation_of_the_equinoxes_complimentary_terms(tt)
         eq_eq = d_psi * cos(self._mean_obliquity_radians) + c_terms
+        # TODO: constrain to 24 hours?
         return self.gmst + eq_eq / tau * 24.0
 
     # Low-precision floats generated from internal float pairs.