From e4b815e2f490b96e48d659d230b6f48782b9ff5a Mon Sep 17 00:00:00 2001 From: Carlo Cena Date: Fri, 16 Oct 2026 17:49:47 +0200 Subject: [PATCH] [#1621] Fix GCRS to J2000 frame bias in tleHandling --- docs/source/Support/bskKnownIssues.rst | 4 +++ .../1621-tle-frame-bias.rst | 1 + src/tests/test_readTLE.py | 32 +++++++++++++++++ src/utilities/tleHandling.py | 36 +++++++++++++++---- 4 files changed, 67 insertions(+), 6 deletions(-) create mode 100644 docs/source/Support/bskReleaseNotesSnippets/1621-tle-frame-bias.rst diff --git a/docs/source/Support/bskKnownIssues.rst b/docs/source/Support/bskKnownIssues.rst index cfe832bf708..c6377f32604 100644 --- a/docs/source/Support/bskKnownIssues.rst +++ b/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. diff --git a/docs/source/Support/bskReleaseNotesSnippets/1621-tle-frame-bias.rst b/docs/source/Support/bskReleaseNotesSnippets/1621-tle-frame-bias.rst new file mode 100644 index 00000000000..32f20e8f592 --- /dev/null +++ b/docs/source/Support/bskReleaseNotesSnippets/1621-tle-frame-bias.rst @@ -0,0 +1 @@ +- Fixed the GCRS to J2000 frame bias in ``tleHandling.frame_bias()``. diff --git a/src/tests/test_readTLE.py b/src/tests/test_readTLE.py index 882f5a2fa33..9dea194ca92 100644 --- a/src/tests/test_readTLE.py +++ b/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() diff --git a/src/utilities/tleHandling.py b/src/utilities/tleHandling.py index dbe8e90dc5b..b9049353399 100644 --- a/src/utilities/tleHandling.py +++ b/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)