From a880292318e0db5c3cd42e331f153c1c1b5c9109 Mon Sep 17 00:00:00 2001 From: Dan Garden Date: Sun, 25 Oct 2026 22:32:35 +1000 Subject: [PATCH] Published test vectors and perihelion_advance (physica #3, item 14) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - relativity::perihelion_advance(gm, a, e) = 6πGM/(c²a(1 − e²)) per orbit, from Will (2014) Eq. 64 with γ = β = 1. With JPL's J2000 elements it gives Mercury's 42.98″ per century (Will's Eq. 65). - New tests, each value read from its paper: - fractional_frequency_shift: JILA's 1.09e-19/mm at g = 9.796 and Wisconsin's 10.9e-19/cm at g = 9.803 (both measured within 1σ). - The solar redshift: 636.31, 633.35 and 633.10 m/s (González Hernández 2020). - shapiro_delay: Cassini's round trip, 262.4 µs, against Will's Eq. 63. - geodetic_precession: the lunar de Sitter precession, 1.919″ per century ("about 2"). - perihelion_advance: S2 from GRAVITY 2020's Table E.1, 12.16′ against their 12.1′. - cargo fmt over the crate; lib.rs and the 0.5.x frame-dragging test had drifted from it, and CI checks fmt. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 8 ++++ README.md | 10 ++++- src/lib.rs | 1 - src/relativity.rs | 111 +++++++++++++++++++++++++++++++++++++++++++--- src/transfer.rs | 16 +++++++ 5 files changed, 138 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5e83ddc..96c56c8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,14 @@ All notable changes to `geo5d`. Plain semver; tags `vX.Y.Z` at the published commit. +## Unreleased + +- **New:** `relativity::perihelion_advance`, `6πGM/(c²a(1 − e²))` per orbit. +- **Tests:** published values for the existing formulas, each read from its paper: JILA's and + Wisconsin's clock gradients, the solar redshift seen from the Earth, Cassini's Shapiro delay, + the lunar de Sitter precession, Mercury's perihelion and S2's Schwarzschild precession. +- `cargo fmt` over the crate, which CI checks. + ## 0.5.1 - Gravity Probe B check uses its real orbit (a = 7027.4 km) and the published predictions diff --git a/README.md b/README.md index c35230e..18af20f 100644 --- a/README.md +++ b/README.md @@ -42,7 +42,7 @@ assert!((h - 1000.0).abs() < 0.1); | `NormalGravity` | Gravity of the reference ellipsoid (GRS80, WGS84) on and above it ([Moritz 2000](https://doi.org/10.1007/s001900050278)) | | `GeopotentialNumber` | Heights from the geopotential: dynamic, normal and orthometric | | `SphericalHarmonicModel` | Gravity, geoid, gravity anomalies, deflections and the gradient tensor from a model such as an ICGEM `.gfc` file, through [`geopotential`](https://crates.io/crates/geopotential) | -| `relativity` | Clocks as height meters (chronometric levelling), clock rates in TT and TCG ([IERS Conventions 2010](https://www.iers.org/TN36)), and frame dragging: the Lense–Thirring drag of satellite orbits and the gyroscope precessions Gravity Probe B measured ([Everitt et al. 2011](https://arxiv.org/abs/1105.3456)) | +| `relativity` | Clocks as height meters (chronometric levelling), clock rates in TT and TCG ([IERS Conventions 2010](https://www.iers.org/TN36)), frame dragging: the Lense–Thirring drag of satellite orbits and the gyroscope precessions Gravity Probe B measured ([Everitt et al. 2011](https://arxiv.org/abs/1105.3456)), and the relativistic advance of a periapsis ([Will 2014](https://arxiv.org/abs/1403.7377)) | | `transfer` | Relativistic corrections for GNSS and time transfer (Sagnac, Shapiro), and proper time along a trajectory | | `constants` | The defining constants, including both conventional values of W₀ | @@ -74,6 +74,14 @@ println!("{:.2} mas/yr", drag * MAS_PER_YEAR); // 30.67 - **GOCO06s:** on the degree-300 model, its synthesis matches `geopotential` to about 1e-15. - **Frame dragging:** LAGEOS's node drag (30.67 mas/yr) and Gravity Probe B's predicted drifts, 6,606 mas/yr geodetic and 39.2 mas/yr frame dragging projected on its guide star. +- **Clocks across a sample:** the redshift gradients JILA and Wisconsin measured, 1.09×10⁻¹⁹ per + mm ([Bothwell et al. 2022](https://arxiv.org/abs/2109.12238)) and 10.9×10⁻¹⁹ per cm + ([Zheng et al. 2023](https://arxiv.org/abs/2207.07145)); the Sun's redshift seen from the + Earth, 633.10 m/s ([González Hernández et al. 2020](https://arxiv.org/abs/2009.10558)). +- **Orbits and light:** Mercury's 42.98″ per century and the Moon's de Sitter precession of about + 2″ per century ([Will 2014](https://arxiv.org/abs/1403.7377)); S2's 12.1′ per orbit round + Sgr A* ([GRAVITY 2020](https://arxiv.org/abs/2004.07187)); Cassini's 262.4 µs round-trip + Shapiro delay. **Not included yet:** - a classical-geoid option from topographic heights; diff --git a/src/lib.rs b/src/lib.rs index c11c6ac..2c108c9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,5 +1,4 @@ #![doc = include_str!("../README.md")] - #![forbid(unsafe_code)] #![warn(missing_docs)] diff --git a/src/relativity.rs b/src/relativity.rs index fd4581f..4653850 100644 --- a/src/relativity.rs +++ b/src/relativity.rs @@ -160,6 +160,23 @@ pub fn gyroscope_frame_dragging(j: f64, a: f64) -> f64 { crate::constants::G_NEWTON * j / (2.0 * C2 * a.powi(3)) } +/// The relativistic advance of an orbit's periapsis, `6πGM/(c²a(1 − e²))` radians per orbit, for +/// semi-major axis `a` (m) and eccentricity `e` about a body of gravitational parameter `gm`: the +/// first term of [Will (2014)](https://arxiv.org/abs/1403.7377), Eq. 64, with `γ = β = 1`. +/// +/// ``` +/// use geo5d::{constants::GM_SUN, relativity::perihelion_advance}; +/// // Mercury (JPL's mean elements at J2000): 42.98″ per century (Will 2014, Eq. 65). +/// let (au, arcsec) = (149_597_870_700.0, core::f64::consts::PI / 648_000.0); +/// let (a, e) = (0.387_099_27 * au, 0.205_635_93); +/// let period = 2.0 * core::f64::consts::PI * (a * a * a / GM_SUN).sqrt(); +/// let per_century = perihelion_advance(GM_SUN, a, e) * 36_525.0 * 86_400.0 / period / arcsec; +/// assert!((per_century - 42.98).abs() < 0.005); +/// ``` +pub fn perihelion_advance(gm: f64, a: f64, e: f64) -> f64 { + 6.0 * core::f64::consts::PI * gm / (C2 * a * (1.0 - e * e)) +} + #[cfg(test)] mod tests { use super::*; @@ -214,20 +231,102 @@ mod tests { assert!(close(per_day, 9.4e-10, 1e-11), "{per_day}"); } + #[test] + fn clock_gradients_across_a_millimetre_and_a_centimetre() { + // A clock Δh higher sits at ΔW = −gΔh and runs faster by gΔh/c². Bothwell et al. (2022, + // arXiv:2109.12238) at g = 9.796 m/s² (USGS): 1.09×10⁻¹⁹ per mm, measured 0.98(23)×10⁻¹⁹. + let jila = fractional_frequency_shift(-9.796 * 1e-3); + assert!( + close(jila, 1.09e-19, 5e-22) && (jila - 0.98e-19).abs() < 2.3e-20, + "{jila}" + ); + // Zheng et al. (2023, arXiv:2207.07145) at g = 9.803 m/s²: 10.9×10⁻¹⁹ per cm, measured + // 12.4 ± 0.7 ± 2.5, 0.57σ away. + let wisconsin = fractional_frequency_shift(-9.803 * 1e-2); + assert!(close(wisconsin, 10.9e-19, 5e-21), "{wisconsin}"); + assert!(((12.4e-19 - wisconsin) / 0.7f64.hypot(2.5) / 1e-19 - 0.575).abs() < 0.001); + } + + #[test] + fn the_suns_redshift_seen_from_the_earth() { + // González Hernández et al. (2020, arXiv:2009.10558): GM☉/(cR☉) = 636.31 m/s at the + // photosphere, 633.35 at 1 AU, 633.10 at the Earth on 2010 November 25 (147.75 million km, + // less the Earth's own 0.21; the text calls the distance correction +4 cm/s, but it lowers + // the value); measured 638 ± 6 m/s. + use crate::constants::{GM_EARTH, GM_SUN, SPEED_OF_LIGHT as C}; + let (r_sun, r_earth, au) = (6.957e8, 6.3781e6, 149_597_870_700.0); // IAU 2015 nominal + let speed = |w_clock: f64, w_ref: f64| -C * fractional_frequency_shift(w_clock - w_ref); + assert!(close(speed(GM_SUN / r_sun, 0.0), 636.31, 0.005)); + assert!(close(speed(GM_SUN / r_sun, GM_SUN / au), 633.35, 0.005)); + let seen = speed(GM_SUN / r_sun, GM_SUN / 147.75e9 + GM_EARTH / r_earth); + assert!(close(seen, 633.10, 0.01) && (638.0 - seen).abs() < 6.0, "{seen}"); + } + + #[test] + fn precession_of_orbits_and_gyroscopes() { + use crate::constants::GM_SUN; + let (au, arcsec) = (149_597_870_700.0, core::f64::consts::PI / 648_000.0); + let century = 36_525.0 * 86_400.0; + // The Earth–Moon system as a gyroscope round the Sun: "about 2 arcseconds per century" + // (Will 2014, after Eq. 70). + let de_sitter = geodetic_precession(GM_SUN, au) * century / arcsec; + assert!( + close(de_sitter, 1.919, 0.001) && close(de_sitter, 2.0, 0.1), + "{de_sitter}" + ); + // S2 round Sgr A* (GRAVITY 2020, arXiv:2004.07187): M = 4.261×10⁶ M☉, R₀ = 8246.7 pc, + // a = 125.058 mas, e = 0.884649 (Table E.1) give 12.16′ per orbit, against their 12.1′ + // (Eq. 1); they measured f_SP = 1.10 ± 0.19 times it. + let a = 125.058e-3 * 8246.7 * au; // an angle in arcseconds at a distance in parsecs is AU + let s2 = perihelion_advance(4.261e6 * GM_SUN, a, 0.884_649) / arcsec / 60.0; + assert!(close(s2, 12.156, 0.001) && close(s2, 12.1, 0.06), "{s2}"); + // A circle advances 6πGM/(c²a), and the advance diverges as e → 1 at fixed a. + assert!(close( + perihelion_advance(GM_SUN, au, 0.0), + 6.0 * core::f64::consts::PI * GM_SUN / (C2 * au), + 1e-25 + )); + assert!(perihelion_advance(GM_SUN, au, 0.99) > 50.0 * perihelion_advance(GM_SUN, au, 0.0)); + } + #[test] fn frame_dragging_lageos_and_gravity_probe_b() { use crate::constants::{EARTH_ANGULAR_MOMENTUM as J, GM_EARTH}; // LAGEOS node: 30.67 mas/yr. Perigee drag is −3 cos i times the node rate. - assert!(close(lense_thirring_node_rate(J, 12_270e3, 0.0045) * MAS_PER_YEAR, 30.67, 0.01)); + assert!(close( + lense_thirring_node_rate(J, 12_270e3, 0.0045) * MAS_PER_YEAR, + 30.67, + 0.01 + )); let i = 109.8f64.to_radians(); - assert!(close(lense_thirring_perigee_rate(J, 12_270e3, 0.0045, i) / lense_thirring_node_rate(J, 12_270e3, 0.0045), -3.0 * i.cos(), 1e-12)); + assert!(close( + lense_thirring_perigee_rate(J, 12_270e3, 0.0045, i) / lense_thirring_node_rate(J, 12_270e3, 0.0045), + -3.0 * i.cos(), + 1e-12 + )); // Gravity Probe B, polar orbit at a = 7027.4 km: 6603.9 mas/yr geodetic and 40.81 mas/yr // frame dragging. The published predictions are 6606.1 and 39.2: the latter is the drift // projected on the guide star IM Pegasi, 40.81 × cos 16.84° = 39.06 (Everitt et al. 2011; // measured −6601.8 ± 18.3 and −37.2 ± 7.2). - assert!(close(geodetic_precession(GM_EARTH, 7_027.4e3) * MAS_PER_YEAR, 6603.9, 0.5)); - assert!(close(gyroscope_frame_dragging(J, 7_027.4e3) * MAS_PER_YEAR, 40.81, 0.02)); - assert!(close(geodetic_precession(GM_EARTH, 7_027.4e3) * MAS_PER_YEAR, 6606.1, 3.0)); - assert!(close(gyroscope_frame_dragging(J, 7_027.4e3) * MAS_PER_YEAR * 16.84f64.to_radians().cos(), 39.2, 0.2)); + assert!(close( + geodetic_precession(GM_EARTH, 7_027.4e3) * MAS_PER_YEAR, + 6603.9, + 0.5 + )); + assert!(close( + gyroscope_frame_dragging(J, 7_027.4e3) * MAS_PER_YEAR, + 40.81, + 0.02 + )); + assert!(close( + geodetic_precession(GM_EARTH, 7_027.4e3) * MAS_PER_YEAR, + 6606.1, + 3.0 + )); + assert!(close( + gyroscope_frame_dragging(J, 7_027.4e3) * MAS_PER_YEAR * 16.84f64.to_radians().cos(), + 39.2, + 0.2 + )); } } diff --git a/src/transfer.rs b/src/transfer.rs index 3507edc..94c2cd1 100644 --- a/src/transfer.rs +++ b/src/transfer.rs @@ -128,6 +128,22 @@ mod tests { assert!(shap > 1e-12 && shap < 1e-9, "{shap}"); } + #[test] + fn cassini_round_trip_shapiro_delay() { + // Cassini's 2002 conjunction (Will 2014, arXiv:1403.7377): the spacecraft at 8.43 AU, the + // ray 1.6 R☉ from the Sun. Twice the one-way delay is 262.4 µs; Will's Eq. 63 shorthand, + // ½(1 + γ)[240 − 20 ln(d²/r)] µs, gives 263.8. + let (au, b): (f64, f64) = (149_597_870_700.0, 1.6 * 6.957e8); + let (r1, r2) = (au, 8.43 * au); + let r12 = (r1 * r1 - b * b).sqrt() + (r2 * r2 - b * b).sqrt(); + let round_trip = 2.0 * shapiro_delay(crate::constants::GM_SUN, r1, r2, r12) * 1e6; + let will = 240.0 - 20.0 * (1.6f64.powi(2) / 8.43).ln(); + assert!( + close(round_trip, 262.42, 0.01) && close(round_trip / will, 1.0, 0.006), + "{round_trip} µs" + ); + } + #[test] fn one_way_transfer_is_dominated_by_light_time() { let sat = Ecef::new(26_000_000.0, 5_000_000.0, 3_000_000.0);