physics: add the TEME to Earth-fixed transformation

crates/apsis-rs/Cargo.toml

@@ -189,3 +189,7 @@ required-features = ["network"]
 [[example]]
 name = "dump_iau2000"
 required-features = ["physics"]
+
+[[example]]
+name = "dump_teme"
+required-features = ["physics"]

crates/apsis-rs/examples/dump_teme.rs

@@ -0,0 +1,36 @@
+//! Writes the TEME to Earth-fixed transformation computed by this workspace,
+//! for cross-checking against astropy.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features physics --example dump_teme \
+//!     > validation/data/expected/rust_teme.csv
+//! ```
+
+use apsis_frames::iau2000::{PolarMotion, teme_to_terrestrial_matrix};
+
+/// The samples, as a UTC Modified Julian Date and a TEME position in
+/// kilometers.
+const SAMPLES: &[(f64, [f64; 3])] = &[
+    (60_481.5, [7_000.0, 1_000.0, 2_000.0]),
+    (51_544.5, [7_000.0, -3_000.0, 1_000.0]),
+    (40_587.0, [42_164.0, 0.0, 0.0]),
+    (60_000.25, [6_878.137, 100.0, -50.0]),
+];
+
+fn main() {
+    println!("mjd,x_ecef,y_ecef,z_ecef");
+    for &(utc_mjd, position) in SAMPLES {
+        // The transformation needs the UT1 time, because it is the Earth's
+        // rotation that carries the frame. UT1 differs from UTC by under
+        // 0.9 seconds, which is 0.0037 degrees of rotation, so the UTC value
+        // is used directly here and the check tolerates that difference.
+        // The Julian date is split into a whole number of days and a
+        // fraction, which preserves more of the day than a single value.
+        let matrix = teme_to_terrestrial_matrix(2_400_000.5, utc_mjd, PolarMotion::NONE);
+        let vector = apsis_math::Vec3::new(position[0], position[1], position[2]);
+        let result = matrix * vector;
+        println!("{utc_mjd},{:.9},{:.9},{:.9}", result.x, result.y, result.z);
+    }
+}

crates/physics/apsis-frames/src/iau2000/mod.rs

@@ -37,6 +37,7 @@ mod fundamental;
 mod nutation;
 mod nutation_table;
 mod precession;
+mod teme;
 mod terrestrial;
 
 pub use fundamental::{
@@ -52,6 +53,7 @@ pub use precession::{
     bias_precession_matrix, frame_bias_matrix, frame_to_mean_of_date, mean_obliquity,
     precession_angles, rotation_z,
 };
+pub use teme::{greenwich_mean_sidereal_time_1982, teme_to_terrestrial_matrix};
 pub use terrestrial::{
     PolarMotion, bias_precession_nutation_matrix, celestial_to_terrestrial_matrix,
     earth_rotation_angle, equation_of_the_equinoxes, equation_of_the_equinoxes_complementary_terms,

crates/physics/apsis-frames/src/iau2000/teme.rs

@@ -0,0 +1,105 @@
+//! The True Equator Mean Equinox frame and its transformation to the
+//! Earth-fixed frame.
+//!
+//! The SGP4 propagator produces its state in the True Equator Mean Equinox
+//! (TEME) frame. TEME is the true equator of date, so its third axis is the
+//! true celestial pole and it carries the nutation, but its first axis
+//! points at the mean equinox rather than at the true equinox.
+//!
+//! Because the first axis is the mean equinox, the transformation to the
+//! Earth-fixed frame uses Greenwich mean sidereal time, not the apparent
+//! value. That is the detail that is easy to get wrong: subtracting the
+//! equation of the equinoxes as well, as one would for a true-equinox frame,
+//! rotates the result by about a thousandth of a degree, which is a
+//! kilometre of ground position for a low Earth satellite.
+//!
+//! The time is therefore the UT1 Julian date, from which the mean sidereal
+//! time follows, and the polar motion is applied last.
+//!
+//! Source: the relation between TEME and the Earth-fixed frame is given in
+//! Vallado, D. A., *Fundamentals of Astrodynamics and Applications*, 4th
+//! edition, microcosm Press, 2013, section 3.7. The sidereal time
+//! expression is the one used by the SGP4 reference implementation, as
+//! given in the same work and implemented by `eraGmst82`.
+
+use super::terrestrial::PolarMotion;
+use apsis_math::Mat3;
+use core::f64::consts::PI;
+
+/// The number of seconds in a day.
+const SECONDS_PER_DAY: f64 = 86_400.0;
+
+/// Converts seconds of time to radians.
+///
+/// A full turn is 86 400 seconds of time, so one second of time is
+/// `2 pi / 86400` radians. Sidereal time is quoted in seconds of time, not in
+/// arcseconds, and the two differ by a factor of fifteen.
+const SECONDS_OF_TIME_TO_RADIAN: f64 = 2.0 * PI / SECONDS_PER_DAY;
+
+/// Returns Greenwich mean sidereal time from the 1982 expression, in
+/// radians.
+///
+/// The expression is a polynomial in the number of Julian centuries from
+/// J2000.0 plus the fraction of the day in seconds of time. It is the
+/// sidereal time the SGP4 theory was fitted against, which is why it is used
+/// here rather than the modern expression: the element sets are consistent
+/// with it.
+///
+/// The Julian date is taken in two parts, a whole number of days and a
+/// fraction, because a single `f64` near 2 461 000 has a resolution of about
+/// 40 microseconds and the sidereal time is wanted to better than that. The
+/// caller that has only one value may pass the whole part and the fraction
+/// separately, or pass a value twice: `(jd, 0.0)`, `(2400000.5, mjd)`, and
+/// `(jd - 0.5, 0.5)` are all accepted, but the split with the smaller whole
+/// part preserves more of the fraction.
+///
+/// Source: as implemented by `eraGmst82`.
+pub fn greenwich_mean_sidereal_time_1982(day: f64, fraction: f64) -> f64 {
+    // The constant absorbs the half day between the start of the Julian
+    // date and the start of the day.
+    let a = 24_110.548_41 - SECONDS_PER_DAY / 2.0;
+    let b = 8_640_184.812_866;
+    let c = 0.093_104;
+    let d = -6.2e-6;
+
+    // The reference orders the two parts so that the smaller is the first,
+    // which is the arrangement that preserves the most precision.
+    let (first, second) = if day < fraction {
+        (day, fraction)
+    } else {
+        (fraction, day)
+    };
+
+    let t = (first + (second - 2_451_545.0)) / 36_525.0;
+    let seconds_of_time = SECONDS_PER_DAY * (first.rem_euclid(1.0) + second.rem_euclid(1.0));
+    let seconds = (a + (b + (c + d * t) * t) * t) + seconds_of_time;
+    (seconds * SECONDS_OF_TIME_TO_RADIAN).rem_euclid(2.0 * PI)
+}
+
+/// Returns the rotation from the True Equator Mean Equinox frame to the
+/// Earth-fixed frame.
+///
+/// The argument is the UT1 Julian date, because the rotation is the Earth's
+/// rotation, and the polar motion is supplied by the caller because it is
+/// measured rather than modelled.
+///
+/// This is the transformation to apply to the position and velocity that
+/// `apsis-sgp4` produces.
+pub fn teme_to_terrestrial_matrix(
+    ut1_day: f64,
+    ut1_fraction: f64,
+    polar_motion: PolarMotion,
+) -> Mat3 {
+    let sidereal = greenwich_mean_sidereal_time_1982(ut1_day, ut1_fraction);
+    // The composition applies the Earth rotation first and the polar motion
+    // last, with each new rotation on the left of the accumulator.
+    polar_motion.matrix() * rotation_z(sidereal)
+}
+
+/// Returns a rotation of `angle` about the third axis, in radians.
+///
+/// The convention is the one the reference helpers use: the rotation is
+/// applied to the frame, so composing `a * b` applies `b` first.
+fn rotation_z(angle: f64) -> Mat3 {
+    super::precession::rotation_z(angle)
+}

crates/physics/apsis-frames/tests/iau2000.rs

@@ -396,3 +396,114 @@ fn a_matrix_at_two_epochs_separated_by_a_century_differs_by_the_precession() {
         "the rotation angle is {degrees} degrees"
     );
 }
+
+#[test]
+fn the_1982_sidereal_time_matches_the_reference() {
+    // Greenwich mean sidereal time from the 1982 expression, which is the
+    // one the SGP4 theory was fitted against. The values are from ERFA.
+    let cases = [
+        (51_544.5, 4.894_961_212_823_059),
+        (60_476.999_999_810_5, 4.620_858_198_832_686),
+        (40_587.0, 1.749_337_177_344_97),
+    ];
+    for (mjd, expected) in cases {
+        let actual = apsis_frames::iau2000::greenwich_mean_sidereal_time_1982(2_400_000.5, mjd);
+        assert!(
+            (actual - expected).abs() < 1e-15,
+            "MJD {mjd}: {actual} against {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_1982_sidereal_time_is_a_full_turn_per_day() {
+    // Sidereal time advances by one turn plus the excess of the sidereal
+    // day over the solar day, which is 0.9856 degrees per day.
+    let first = apsis_frames::iau2000::greenwich_mean_sidereal_time_1982(2_400_000.5, 60_000.0);
+    let second = apsis_frames::iau2000::greenwich_mean_sidereal_time_1982(2_400_000.5, 60_001.0);
+    let advance = (second - first).rem_euclid(2.0 * PI);
+    assert!(
+        (advance.to_degrees() - 0.985_6).abs() < 0.001,
+        "advance {} degrees",
+        advance.to_degrees()
+    );
+}
+
+#[test]
+fn the_teme_transformation_to_earth_fixed_matches_the_reference() {
+    use apsis_frames::iau2000::{PolarMotion, teme_to_terrestrial_matrix};
+    use apsis_math::Vec3;
+
+    // The reference values are ERFA's `gmst82` applied as a third axis
+    // rotation with the identity polar motion.
+    let cases = [
+        (
+            60_481.5,
+            [7_000.0, 1_000.0, 2_000.0],
+            [1_098.716_242_792, -6_985.185_940_104, 2_000.0],
+        ),
+        (
+            51_544.5,
+            [7_000.0, -3_000.0, 1_000.0],
+            [4_221.057_369_867, 6_338.980_571_061, 1_000.0],
+        ),
+    ];
+    for (mjd, position, expected) in cases {
+        let matrix = teme_to_terrestrial_matrix(2_400_000.5, mjd, PolarMotion::NONE);
+        let result = matrix * Vec3::new(position[0], position[1], position[2]);
+        let expected = Vec3::new(expected[0], expected[1], expected[2]);
+        assert!(
+            (result - expected).norm() < 1e-6,
+            "MJD {mjd}: {result} against {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_teme_transformation_preserves_length() {
+    use apsis_frames::iau2000::{PolarMotion, teme_to_terrestrial_matrix};
+    use apsis_math::Vec3;
+
+    let matrix = teme_to_terrestrial_matrix(2_400_000.5, 60_000.0, PolarMotion::NONE);
+    assert!(matrix.is_rotation(1e-12));
+    let vector = Vec3::new(7_000.0, -1_000.0, 500.0);
+    assert!(
+        ((matrix * vector).norm() - vector.norm()).abs() < 1e-9,
+        "the length changed"
+    );
+}
+
+#[test]
+fn the_teme_frame_uses_the_mean_equinox_not_the_true_one() {
+    // The distinction matters. TEME's first axis is the mean equinox, so it
+    // pairs with a mean sidereal time. Substituting the apparent value, or
+    // the modern mean value, rotates the result by a fraction of a degree,
+    // which is many kilometres of ground position.
+    use apsis_frames::iau2000::{
+        greenwich_apparent_sidereal_time, greenwich_mean_sidereal_time,
+        greenwich_mean_sidereal_time_1982,
+    };
+
+    let mjd = 60_481.5;
+    let tt = mjd + 2_400_000.5 + 69.184 / 86_400.0;
+    let sidereal_1982 = greenwich_mean_sidereal_time_1982(2_400_000.5, mjd);
+    let apparent = greenwich_apparent_sidereal_time(tt, tt);
+    let modern_mean = greenwich_mean_sidereal_time(tt, tt);
+
+    // The apparent value includes the equation of the equinoxes, which is a
+    // few thousandths of a degree.
+    let equation = (apparent - modern_mean).abs().to_degrees();
+    assert!(
+        equation > 1e-4 && equation < 1e-2,
+        "the equation of the equinoxes is {equation} degrees"
+    );
+
+    // The 1982 expression, which TEME is fitted against, differs from the
+    // modern mean value by a fraction of a degree, because the precession
+    // model changed between them.
+    let secular = (sidereal_1982 - modern_mean).abs().to_degrees();
+    assert!(
+        secular > 0.1 && secular < 1.0,
+        "the difference is {secular} degrees"
+    );
+}

validation/README.md

@@ -86,6 +86,10 @@ cargo run -p apsis-rs --features rf --example dump_gas_exact \
 # IAU 2000/2006 地球定向（岁差、章动、框架偏移、GCRS→ITRF）
 cargo run -p apsis-rs --features physics --example dump_iau2000 \
     > validation/data/expected/rust_iau2000.csv
+
+# TEME → 地固系（SGP4 输出链路）
+cargo run -p apsis-rs --features physics --example dump_teme \
+    > validation/data/expected/rust_teme.csv
 ```
 
 Python 侧读取同一批输入，用参考实现重算并输出误差表：
@@ -102,6 +106,7 @@ pixi run python validation/python/link/check_s465.py
 pixi run python validation/python/atmosphere/check_p839.py
 pixi run python validation/python/atmosphere/check_gas_exact.py
 pixi run python validation/python/time/check_iau2000.py
+pixi run python validation/python/coord/check_teme.py
 ```
 
 两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。
@@ -134,6 +139,7 @@ Earth orientation 数据自动下载；本 crate 涉及的换算不需要该数
 | `atmosphere/check_p839.py` | 零度等温线高度，15 个全球点 | `itur` 的 ITU-R P.839-4 实现 | 3.3e-13 km |
 | `atmosphere/check_gas_exact.py` | 精确斜路径，11 个场景（1–100 GHz，1°–90° 仰角） | `itur` 的 ITU-R P.676-12 exact 方法 | 9.6e-9（相对） |
 | `time/check_iau2000.py` | 岁差、章动、框架偏移、GAST、GCRS→ITRF，8 个历元 × 28 个量 | ERFA 2.0.1（IAU SOFA 派生） | 章动 6.6e-15 rad；BPN 矩阵 6.0e-15；GAST 1.1e-11 rad |
+| `coord/check_teme.py` | TEME → 地固系，4 个历元 | ERFA 的 `gmst82` 与 `pom00`（同为 0 极移） | 4.9e-10 km |
 
 ### 关于射频层对拍的说明
 

validation/python/coord/check_teme.py

@@ -0,0 +1,119 @@
+"""Cross-check the ApsisRS TEME to Earth-fixed transformation.
+
+Reads the positions written by ``cargo run -p apsis-rs --features physics
+--example dump_teme`` and recomputes each one with ERFA, applying the same
+Greenwich mean sidereal time and the same zero polar motion, so that the
+check compares the implementation rather than the choice of inputs.
+
+The reference is the sidereal time of 1982 and the polar motion composition
+of the IAU SOFA library, as implemented by ERFA. Note that astropy agrees
+with this transformation only to a few hundred metres, because astropy
+applies the measured polar motion and a measured UT1 - UTC while this check
+uses zero for both. That agreement is reported for information but is not
+the acceptance criterion.
+
+The script needs no network access.
+
+Usage:
+
+    pixi run python validation/python/coord/check_teme.py [path]
+
+Exit status is zero when every position agrees within tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import sys
+from pathlib import Path
+
+import erfa
+
+# Tolerance in kilometers. The transformation is a product of two rotations
+# evaluated in the same order, so the agreement is at the level of the last
+# bits of a double.
+TOLERANCE_KM = 1e-9
+
+# The two-part Julian date used by the Rust example for the whole days.
+DAY_PART = 2400000.5
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_teme.csv")
+
+
+def read_rows(path: Path) -> list[dict[str, str]]:
+    """Reads the CSV rows."""
+    with path.open(newline="") as handle:
+        reader = csv.DictReader(handle)
+        return [row for row in reader if row.get("mjd")]
+
+
+def reference_position(mjd: float, position: list[float]) -> list[float]:
+    """Recomputes the Earth-fixed position with ERFA."""
+    # The sidereal time, from the two-part Julian date.
+    sidereal = erfa.gmst82(DAY_PART, mjd)
+    # The rotation of the frame about its third axis, composed as the
+    # reference `rz` does, and the identity polar motion.
+    rotation = erfa.rz(sidereal, erfa.ir())
+    polar = erfa.pom00(0.0, 0.0, 0.0)
+    matrix = erfa.rxr(polar, rotation)
+
+    # The rotation maps the frame, so the position is the matrix applied to
+    # the components.
+    return [
+        sum(matrix[i][j] * position[j] for j in range(3)) for i in range(3)
+    ]
+
+
+def main(argv: list[str]) -> int:
+    path = Path(argv[1]) if len(argv) > 1 else DEFAULT_INPUT
+    if not path.exists():
+        print(f"missing input file: {path}", file=sys.stderr)
+        print(
+            "generate it with: cargo run -p apsis-rs --features physics "
+            "--example dump_teme > " + str(path),
+            file=sys.stderr,
+        )
+        return 2
+
+    rows = read_rows(path)
+    if not rows:
+        print("input file has no data rows", file=sys.stderr)
+        return 2
+
+    # The TEME inputs used by the example, in the same order.
+    inputs = [
+        [7000.0, 1000.0, 2000.0],
+        [7000.0, -3000.0, 1000.0],
+        [42164.0, 0.0, 0.0],
+        [6878.137, 100.0, -50.0],
+    ]
+
+    worst = 0.0
+    worst_at = None
+    for row, position in zip(rows, inputs):
+        mjd = float(row["mjd"])
+        actual = [
+            float(row["x_ecef"]),
+            float(row["y_ecef"]),
+            float(row["z_ecef"]),
+        ]
+        expected = reference_position(mjd, position)
+        error = max(abs(a - b) for a, b in zip(actual, expected))
+        if error > worst:
+            worst = error
+            worst_at = mjd
+
+    print(f"positions compared: {len(rows)}")
+    print(f"worst difference:   {worst:.3e} km at MJD {worst_at}")
+    print(f"tolerance:          {TOLERANCE_KM:.3e} km")
+
+    if worst <= TOLERANCE_KM:
+        print("result:             pass")
+        return 0
+    print("result:             FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))