From 9d8773d6abdf652a75dba706c2a03cbe5507413c Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 16:54:09 +0800 Subject: [PATCH] physics: add the TEME to Earth-fixed transformation --- crates/apsis-rs/Cargo.toml | 4 + crates/apsis-rs/examples/dump_teme.rs | 36 ++++++ .../physics/apsis-frames/src/iau2000/mod.rs | 2 + .../physics/apsis-frames/src/iau2000/teme.rs | 105 ++++++++++++++++ crates/physics/apsis-frames/tests/iau2000.rs | 111 ++++++++++++++++ validation/README.md | 6 + validation/python/coord/check_teme.py | 119 ++++++++++++++++++ 7 files changed, 383 insertions(+) create mode 100644 crates/apsis-rs/examples/dump_teme.rs create mode 100644 crates/physics/apsis-frames/src/iau2000/teme.rs create mode 100644 validation/python/coord/check_teme.py diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 9b3587a..2f10875 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/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"] diff --git a/crates/apsis-rs/examples/dump_teme.rs b/crates/apsis-rs/examples/dump_teme.rs new file mode 100644 index 0000000..db4affb --- /dev/null +++ b/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); + } +} diff --git a/crates/physics/apsis-frames/src/iau2000/mod.rs b/crates/physics/apsis-frames/src/iau2000/mod.rs index 8ddce8d..76e7dbb 100644 --- a/crates/physics/apsis-frames/src/iau2000/mod.rs +++ b/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, diff --git a/crates/physics/apsis-frames/src/iau2000/teme.rs b/crates/physics/apsis-frames/src/iau2000/teme.rs new file mode 100644 index 0000000..eb3c363 --- /dev/null +++ b/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) +} diff --git a/crates/physics/apsis-frames/tests/iau2000.rs b/crates/physics/apsis-frames/tests/iau2000.rs index 5732302..ab91299 100644 --- a/crates/physics/apsis-frames/tests/iau2000.rs +++ b/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" + ); +} diff --git a/validation/README.md b/validation/README.md index 7d919ef..13437fa 100644 --- a/validation/README.md +++ b/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 | ### 关于射频层对拍的说明 diff --git a/validation/python/coord/check_teme.py b/validation/python/coord/check_teme.py new file mode 100644 index 0000000..c202400 --- /dev/null +++ b/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))