Fix GCRS to J2000 frame bias in tleHandling

docs/source/Support/bskKnownIssues.rst

@@ -7,6 +7,10 @@ Basilisk Known Issues
 Version |release|
 -----------------
 
+- GitHub issue 1621: ``tleHandling.frame_bias()`` applied the GCRS to mean J2000 frame bias rotations in the
+  reverse order and with the wrong sign of :math:`\xi_0`. The TEME to J2000/ICRF conversion used by
+  ``tleHandling`` and ``ommHandling`` therefore shifted TLE and OMM initial states by about 1.3 m in LEO.
+
 - The error message of ``atmosphericDrag()`` in ``orbitalMotion``, printed for a position with a non-positive
   altitude, showed the position components with an index shifted by one and read past the end of the position
   vector.

docs/source/Support/bskReleaseNotesSnippets/1621-tle-frame-bias.rst

@@ -0,0 +1 @@
+- Fixed the GCRS to J2000 frame bias in ``tleHandling.frame_bias()``.

src/tests/test_readTLE.py

@@ -266,6 +266,34 @@ def test_write_tle_6digit_norad_raises():
     with pytest.raises(ValueError, match="100000"):
         tleHandling.generateTle(tleData)
 
+
+# Reference values from Orekit 13.1 (orekit_jpype 13.1.8.0): GCRF to EME2000 rotation at J2000_EPOCH
+OREKIT_GCRF2EME2000 = np.array([
+    [9.9999999999999423e-01, -7.0782797441991980e-08, 8.0562171469761338e-08],
+    [7.0782794778573375e-08, 9.9999999999999711e-01, 3.3060414542221364e-08],
+    [-8.0562173809869717e-08, -3.3060408839805523e-08, 9.9999999999999623e-01],
+])  # [-]
+
+# TEME state and its Orekit 13.1 FramesFactory.getTEME().getTransformTo(FramesFactory.getGCRF(), date) result
+TEME_EPOCH = datetime(2004, 4, 6, 7, 51, 28, 386009)  # UTC
+TEME_R = np.array([5094.18016210, 6127.64465950, 6380.34453270]) * 1e3  # [m]
+TEME_V = np.array([-4.746131487, 0.785818041, 5.531931288]) * 1e3  # [m/s]
+OREKIT_GCRF_R = np.array([5102509.519317023, 6123010.947714332, 6378136.913898652])  # [m]
+OREKIT_GCRF_V = np.array([-4743.219940214507, 790.5366985409497, 5533.7558470487875])  # [m/s]
+
+
+def test_frame_bias_matches_orekit():
+    """frame_bias() must match the Orekit GCRF to EME2000 rotation (IERS frame bias)."""
+    np.testing.assert_allclose(tleHandling.frame_bias(), OREKIT_GCRF2EME2000, rtol=0.0, atol=1e-11)  # [-]
+
+
+def test_teme2j2000_matches_orekit():
+    """_teme2j2000() must match the Orekit TEME to GCRF transform."""
+    r_m, v_m_s = tleHandling._teme2j2000(TEME_R, TEME_V, TEME_EPOCH)
+    np.testing.assert_allclose(r_m, OREKIT_GCRF_R, rtol=0.0, atol=0.15)  # [m]
+    np.testing.assert_allclose(v_m_s, OREKIT_GCRF_V, rtol=0.0, atol=5e-4)  # [m/s]
+
+
 if __name__ == "__main__":
     # Reading TLE
     test_read_tle(DATA_DIR / "hypso1.tle", EXPECTED_OE_HYPSO)
@@ -277,3 +305,7 @@ def test_write_tle_6digit_norad_raises():
 
     # Take TLE, generate orbital elements, then generate TLE again and compare to original
     test_read_write_tle(DATA_DIR / "hypso1.tle")
+
+    # Frame bias and TEME to GCRF conversion against Orekit
+    test_frame_bias_matches_orekit()
+    test_teme2j2000_matches_orekit()

src/utilities/tleHandling.py

@@ -72,14 +72,26 @@ def utc_to_tt(jd_utc):
     return jd_utc + delta_t
 
 def frame_bias():
-    d_alpha = -0.0146 * ARCSEC2RAD
-    xi      = -0.016617 * ARCSEC2RAD
-    eta     = -0.0068192 * ARCSEC2RAD
+    """
+    Frame bias matrix from the GCRS to the mean equator and equinox of J2000.
+
+    Uses the IERS Conventions (2010) form :math:`B = R_1(-\\eta_0) R_2(\\xi_0) R_3(d\\alpha_0)`.
+    Here :math:`R_i(\\theta)` is the axes (frame) rotation returned by :func:`rot`, which rotates the
+    coordinate frame by :math:`\\theta` about axis :math:`i` (for example :math:`R_3(\\theta)` has
+    :math:`+\\sin\\theta` in row 1, column 2); under the vector-rotation convention the same product
+    reads as the transpose.
+
+    Returns:
+        numpy.ndarray: 3x3 rotation matrix that maps GCRS vector components to mean J2000 components.
+    """
+    d_alpha = -0.0146 * ARCSEC2RAD  # [rad] ICRS right ascension offset of the J2000 mean equinox
+    xi      = -0.016617 * ARCSEC2RAD  # [rad] celestial pole offset in x
+    eta     = -0.0068192 * ARCSEC2RAD  # [rad] celestial pole offset in y
 
     return (
-        rot(3, d_alpha) @
-        rot(2, -xi) @
-        rot(1, -eta)
+        rot(1, -eta) @
+        rot(2, xi) @
+        rot(3, d_alpha)
     )
 
 def nutation_iau1980(T):
@@ -297,6 +309,18 @@ def rot(axis, angle):
 
 # Helper functions: Coordinate conversion; TEME -> J2000/ICRS
 def _teme2j2000(r, v, epoch):
+    """
+    Rotate a TEME state into the Basilisk inertial frame.
+
+    The rotation chain is equation of equinoxes, IAU 1980 nutation, IAU 1976 precession, and the inverse of
+    :func:`frame_bias`, so the output is expressed in the GCRF. These are the ICRF-aligned axes that SPICE
+    calls ``J2000`` and that Basilisk uses as its inertial frame.
+
+    :param r: TEME position vector [m]
+    :param v: TEME velocity vector [m/s]
+    :param epoch: UTC epoch of the state as a ``datetime``
+    :return: tuple ``(r, v)`` with the GCRF position [m] and velocity [m/s]
+    """
 
     jd_utc = datetime_to_jd(epoch)
     jd_tt  = utc_to_tt(jd_utc)