ApsisRS
/
apsis-rs
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

physics: add IAU 2000A nutation and IAU 2006 precession

This commit is contained in:
ApsisRS contributors 2026-09-26 16:48:57 +08:00
parent c0bd0400bf
commit 841b586e2d
15 changed files with 7546 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
crates/apsis-rs/Cargo.toml
View File

@ -185,3 +185,7 @@ required-features = ["rf"]
[[example]] [[example]]
name = "integrate_network" name = "integrate_network"
required-features = ["network"] required-features = ["network"]
[[example]]
name = "dump_iau2000"
required-features = ["physics"]

92
crates/apsis-rs/examples/dump_iau2000.rs Normal file
View File

@ -0,0 +1,92 @@
//! Writes the IAU 2000/2006 Earth orientation quantities computed by this
//! workspace, for cross-checking against ERFA.
//!
//! Usage:
//!
//! ```text
//! cargo run -p apsis-rs --features physics --example dump_iau2000 \
//! > validation/data/expected/rust_iau2000.csv
//! ```
use apsis_frames::iau2000::{
PolarMotion, bias_precession_nutation_matrix, celestial_to_terrestrial_matrix,
greenwich_apparent_sidereal_time, mean_obliquity, nutation_2000a,
nutation_2000a_with_2006_adjustment, precession_angles,
};
use apsis_math::Mat3;
/// The Modified Julian Dates sampled. They span from 1970 to 2050 so that
/// the precession is exercised over a range of a century.
const MJD: &[f64] = &[
40_587.0, 51_544.5, 53_000.0, 55_000.0, 57_000.0, 58_000.0, 60_000.0, 62_000.0,
];
/// The Terrestrial Time offset from TAI, in seconds.
const TT_MINUS_TAI: f64 = 32.184;
/// The TAI minus UTC offset used for this check, in seconds. It is the value
/// in force from 2017 onward.
const TAI_MINUS_UTC: f64 = 37.0;
/// Returns the TT Julian date for a Modified Julian Date on the UTC scale.
fn tt_julian_date(utc_mjd: f64) -> f64 {
// TAI = UTC + 37, TT = TAI + 32.184, so TT = UTC + 69.184 seconds.
utc_mjd + 2_400_000.5 + (TT_MINUS_TAI + TAI_MINUS_UTC) / 86_400.0
}
fn main() {
println!("mjd,quantity,index,value");
for &mjd in MJD {
let tt = tt_julian_date(mjd);
let t = (tt - 2_451_545.0) / 36_525.0;
// The obliquity, in radians.
println!("{mjd},obliquity,0,{:.17e}", mean_obliquity(t));
// The nutation, from the IAU 2000A model.
let nutation = nutation_2000a(t);
println!("{mjd},nutation_longitude_a,0,{:.17e}", nutation.longitude);
println!("{mjd},nutation_obliquity_a,0,{:.17e}", nutation.obliquity);
// The nutation with the IAU 2006 adjustment.
let adjusted = nutation_2000a_with_2006_adjustment(t);
println!("{mjd},nutation_longitude_06a,0,{:.17e}", adjusted.longitude);
println!("{mjd},nutation_obliquity_06a,0,{:.17e}", adjusted.obliquity);
// The four precession angles.
let angles = precession_angles(t);
println!("{mjd},precession_gamma_bar,0,{:.17e}", angles.gamma_bar);
println!("{mjd},precession_phi_bar,0,{:.17e}", angles.phi_bar);
println!("{mjd},precession_psi_bar,0,{:.17e}", angles.psi_bar);
println!("{mjd},precession_epsilon,0,{:.17e}", angles.epsilon);
// The bias-precession-nutation matrix, which carries the celestial
// reference frame to the true equator and equinox of date.
let bpn = bias_precession_nutation_matrix(t);
for row in 0..3 {
for column in 0..3 {
let value = bpn.get(row, column).expect("the indices are in range");
println!("{mjd},bpn,{},{:.17e}", row * 3 + column, value);
}
}
// Greenwich apparent sidereal time, in radians.
println!(
"{mjd},gast,0,{:.17e}",
greenwich_apparent_sidereal_time(tt, tt)
);
// The celestial-to-terrestrial matrix, with zero polar motion.
let matrix = celestial_to_terrestrial_matrix(tt, tt, PolarMotion::NONE);
for row in 0..3 {
for column in 0..3 {
let value = matrix.get(row, column).expect("the indices are in range");
println!("{mjd},c2t,{},{:.17e}", row * 3 + column, value);
}
}
}
}
/// Keeps the matrix type in scope for the documentation build.
#[allow(dead_code)]
fn _uses_matrix(_matrix: Mat3) {}

158
crates/physics/apsis-frames/src/iau2000/eect_table.rs Normal file
View File

@ -0,0 +1,158 @@
//! The IAU 2000 complementary terms of the equation of the equinoxes.
//!
//! The table is generated from the reference implementation by
//! `validation/tools/generate_eect_table.py` and is reproduced without
//! change.
//!
//! Source: the MHB2000 model, as distributed with the IAU SOFA library
//! and implemented by ERFA.
use super::ComplementaryTerm;
/// The terms that do not depend on time.
///
/// The eight multipliers apply to the arguments `l`, `l'`, `F`, `D`, `Om`,
/// the mean longitude of Venus, the mean longitude of the Earth, and the
/// general accumulated precession, in that order.
pub static CONSTANT_TERMS: [ComplementaryTerm; 33] = [
ComplementaryTerm {
argument: [0, 0, 0, 0, 1, 0, 0, 0],
coefficients: [2640.96e-6, -0.39e-6],
},
ComplementaryTerm {
argument: [0, 0, 0, 0, 2, 0, 0, 0],
coefficients: [63.52e-6, -0.02e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, -2, 3, 0, 0, 0],
coefficients: [11.75e-6, 0.01e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, -2, 1, 0, 0, 0],
coefficients: [11.21e-6, 0.01e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, -2, 2, 0, 0, 0],
coefficients: [-4.55e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, 0, 3, 0, 0, 0],
coefficients: [2.02e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, 0, 1, 0, 0, 0],
coefficients: [1.98e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 0, 0, 3, 0, 0, 0],
coefficients: [-1.72e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 1, 0, 0, 1, 0, 0, 0],
coefficients: [-1.41e-6, -0.01e-6],
},
ComplementaryTerm {
argument: [0, 1, 0, 0, -1, 0, 0, 0],
coefficients: [-1.26e-6, -0.01e-6],
},
ComplementaryTerm {
argument: [1, 0, 0, 0, -1, 0, 0, 0],
coefficients: [-0.63e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [1, 0, 0, 0, 1, 0, 0, 0],
coefficients: [-0.63e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 1, 2, -2, 3, 0, 0, 0],
coefficients: [0.46e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 1, 2, -2, 1, 0, 0, 0],
coefficients: [0.45e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 4, -4, 4, 0, 0, 0],
coefficients: [0.36e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 1, -1, 1, -8, 12, 0],
coefficients: [-0.24e-6, -0.12e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, 0, 0, 0, 0, 0],
coefficients: [0.32e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, 0, 2, 0, 0, 0],
coefficients: [0.28e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [1, 0, 2, 0, 3, 0, 0, 0],
coefficients: [0.27e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [1, 0, 2, 0, 1, 0, 0, 0],
coefficients: [0.26e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, -2, 0, 0, 0, 0],
coefficients: [-0.21e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 1, -2, 2, -3, 0, 0, 0],
coefficients: [0.19e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 1, -2, 2, -1, 0, 0, 0],
coefficients: [0.18e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 0, 0, 0, 8, -13, -1],
coefficients: [-0.10e-6, 0.05e-6],
},
ComplementaryTerm {
argument: [0, 0, 0, 2, 0, 0, 0, 0],
coefficients: [0.15e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [2, 0, -2, 0, -1, 0, 0, 0],
coefficients: [-0.14e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [1, 0, 0, -2, 1, 0, 0, 0],
coefficients: [0.14e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 1, 2, -2, 2, 0, 0, 0],
coefficients: [-0.14e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [1, 0, 0, -2, -1, 0, 0, 0],
coefficients: [0.14e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 4, -2, 4, 0, 0, 0],
coefficients: [0.13e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [0, 0, 2, -2, 4, 0, 0, 0],
coefficients: [-0.11e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [1, 0, -2, 0, -3, 0, 0, 0],
coefficients: [0.11e-6, 0.00e-6],
},
ComplementaryTerm {
argument: [1, 0, -2, 0, -1, 0, 0, 0],
coefficients: [0.11e-6, 0.00e-6],
},
];
/// The terms whose amplitude grows linearly with time.
///
/// The arguments are the same as for the constant terms.
pub static LINEAR_TERMS: [ComplementaryTerm; 1] = [ComplementaryTerm {
argument: [0, 0, 0, 0, 1, 0, 0, 0],
coefficients: [-0.87e-6, 0.00e-6],
}];

273
crates/physics/apsis-frames/src/iau2000/fundamental.rs Normal file
View File

@ -0,0 +1,273 @@
//! The fundamental arguments of the IAU 2000 nutation theory.
//!
//! The nutation series is a sum of periodic terms whose arguments are linear
//! functions of time. Each term is identified by integer multipliers applied
//! to these arguments, so the arguments must be evaluated in a fixed order.
//!
//! Two argument sets appear in the theory and they are not the same, which
//! is a detail that is easy to get wrong:
//!
//! - The **luni-solar** series uses the five Delaunay arguments in the order
//! `l`, `l'`, `F`, `D`, `Om`: the mean anomaly of the Moon, the mean
//! anomaly of the Sun, the mean argument of latitude of the Moon, the mean
//! elongation of the Moon from the Sun, and the mean longitude of the
//! Moon's ascending node.
//! - The **planetary** series uses `l`, `F`, `D`, `Om`, the mean longitudes
//! of Mercury through Uranus, the mean longitude of Neptune, and the
//! general accumulated precession in longitude.
//!
//! The reference evaluates `l'` and `D` with inline expressions while using
//! functions for the others, and the inline expressions differ in the last
//! digits from the functions of the same name in the older IERS 2003
//! conventions. The values here are the ones the series was fitted with, so
//! the inline forms are reproduced.
//!
//! All arguments are in radians and are reduced into `[0, 2 pi)`.
//!
//! Source: the argument definitions of the MHB2000 model as implemented in
//! `eraNut00a`, with the planetary longitudes from `eraFame03`,
//! `eraFave03`, `eraFae03`, `eraFama03`, `eraFaju03`, `eraFasa03` and
//! `eraFaur03`.
use core::f64::consts::PI;
/// Converts arcseconds to radians.
const ARCSECOND_TO_RADIAN: f64 = PI / (180.0 * 3_600.0);
/// The number of arcs in a full turn, in arcseconds.
const TURN_ARCSECONDS: f64 = 1_296_000.0;
/// Returns a value reduced into `[0, 2 pi)`.
fn anp(angle: f64) -> f64 {
angle.rem_euclid(2.0 * PI)
}
/// Returns the mean anomaly of the Moon, in radians.
///
/// This is the Delaunay argument `l`, from the IERS 2003 expressions.
pub fn mean_anomaly_moon(t: f64) -> f64 {
let arcseconds = 485_868.249_036
+ t * (1_717_915_923.217_8 + t * (31.879_2 + t * (0.051_635 + t * (-0.000_244_70))));
anp((arcseconds % TURN_ARCSECONDS) * ARCSECOND_TO_RADIAN)
}
/// Returns the mean anomaly of the Sun, in radians.
///
/// This is the Delaunay argument `l'`. The reference evaluates it inline
/// rather than through the general function, and the two differ slightly.
pub fn mean_anomaly_sun(t: f64) -> f64 {
let arcseconds = 1_287_104.793_05
+ t * (129_596_581.048_1 + t * (-0.553_2 + t * (0.000_136 + t * (-0.000_011_49))));
anp((arcseconds % TURN_ARCSECONDS) * ARCSECOND_TO_RADIAN)
}
/// Returns the mean anomaly of the Sun from the IERS 2003 expressions.
///
/// This differs from [`mean_anomaly_sun`] in the last digits of the constant
/// and linear terms. The nutation series uses the inline form; the
/// complementary terms of the equation of the equinoxes use this one.
///
/// Source: as implemented by `eraFalp03`.
pub fn mean_anomaly_sun_iers2003(t: f64) -> f64 {
let arcseconds = 1_287_104.793_048
+ t * (129_596_581.048_1 + t * (-0.553_2 + t * (0.000_136 + t * (-0.000_011_49))));
anp((arcseconds % TURN_ARCSECONDS) * ARCSECOND_TO_RADIAN)
}
/// Returns the mean elongation of the Moon from the Sun from the IERS 2003
/// expressions.
///
/// This differs from [`mean_elongation_moon_sun`] in the last digits of the
/// constant term. The nutation series uses the inline form; the
/// complementary terms of the equation of the equinoxes use this one.
///
/// Source: as implemented by `eraFad03`.
pub fn mean_elongation_moon_sun_iers2003(t: f64) -> f64 {
let arcseconds = 1_072_260.703_692
+ t * (1_602_961_601.209_0 + t * (-6.370_6 + t * (0.006_593 + t * (-0.000_031_69))));
anp((arcseconds % TURN_ARCSECONDS) * ARCSECOND_TO_RADIAN)
}
/// Returns the mean argument of latitude of the Moon, in radians.
///
/// This is the Delaunay argument `F`, the mean longitude of the Moon minus
/// the mean longitude of its ascending node.
pub fn mean_latitude_moon(t: f64) -> f64 {
let arcseconds = 335_779.526_232
+ t * (1_739_527_262.847_8 + t * (-12.751_2 + t * (-0.001_037 + t * (0.000_004_17))));
anp((arcseconds % TURN_ARCSECONDS) * ARCSECOND_TO_RADIAN)
}
/// Returns the mean elongation of the Moon from the Sun, in radians.
///
/// This is the Delaunay argument `D`. The reference evaluates it inline.
pub fn mean_elongation_moon_sun(t: f64) -> f64 {
let arcseconds = 1_072_260.703_69
+ t * (1_602_961_601.209_0 + t * (-6.370_6 + t * (0.006_593 + t * (-0.000_031_69))));
anp((arcseconds % TURN_ARCSECONDS) * ARCSECOND_TO_RADIAN)
}
/// Returns the mean longitude of the Moon's ascending node, in radians.
///
/// This is the Delaunay argument `Om`.
pub fn mean_longitude_node(t: f64) -> f64 {
let arcseconds = 450_160.398_036
+ t * (-6_962_890.543_1 + t * (7.472_2 + t * (0.007_702 + t * (-0.000_059_39))));
anp((arcseconds % TURN_ARCSECONDS) * ARCSECOND_TO_RADIAN)
}
/// Returns the mean longitude of Mercury, in radians.
pub fn mean_longitude_mercury(t: f64) -> f64 {
anp(4.402_608_842 + 2_608.790_314_157_4 * t)
}
/// Returns the mean longitude of Venus, in radians.
pub fn mean_longitude_venus(t: f64) -> f64 {
anp(3.176_146_697 + 1_021.328_554_621_1 * t)
}
/// Returns the mean longitude of the Earth, in radians.
pub fn mean_longitude_earth(t: f64) -> f64 {
anp(1.753_470_314 + 628.307_584_999_1 * t)
}
/// Returns the mean longitude of Mars, in radians.
pub fn mean_longitude_mars(t: f64) -> f64 {
anp(6.203_480_913 + 334.061_242_670_0 * t)
}
/// Returns the mean longitude of Jupiter, in radians.
pub fn mean_longitude_jupiter(t: f64) -> f64 {
anp(0.599_546_497 + 52.969_096_264_1 * t)
}
/// Returns the mean longitude of Saturn, in radians.
pub fn mean_longitude_saturn(t: f64) -> f64 {
anp(0.874_016_757 + 21.329_910_496_0 * t)
}
/// Returns the mean longitude of Uranus, in radians.
pub fn mean_longitude_uranus(t: f64) -> f64 {
anp(5.481_293_872 + 7.478_159_856_7 * t)
}
/// Returns the mean longitude of Neptune, in radians.
///
/// The reference evaluates this inline, from the MHB2000 model rather than
/// from the IERS 2003 expressions used for the other planets.
pub fn mean_longitude_neptune(t: f64) -> f64 {
anp(5.321_159_000 + 3.812_777_400_0 * t)
}
/// Returns the general accumulated precession in longitude, in radians.
pub fn general_precession(t: f64) -> f64 {
(0.024_381_750 + 0.000_005_386_91 * t) * t
}
/// The five Delaunay arguments, in the order the luni-solar series uses.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct DelaunayArguments {
/// The mean anomaly of the Moon.
pub l: f64,
/// The mean anomaly of the Sun.
pub l_prime: f64,
/// The mean argument of latitude of the Moon.
pub f: f64,
/// The mean elongation of the Moon from the Sun.
pub d: f64,
/// The mean longitude of the Moon's ascending node.
pub omega: f64,
}
impl DelaunayArguments {
/// Returns the five arguments for a time in Julian centuries.
pub fn at(t: f64) -> Self {
Self {
l: mean_anomaly_moon(t),
l_prime: mean_anomaly_sun(t),
f: mean_latitude_moon(t),
d: mean_elongation_moon_sun(t),
omega: mean_longitude_node(t),
}
}
/// Returns the arguments as an array in series order.
pub const fn as_array(self) -> [f64; 5] {
[self.l, self.l_prime, self.f, self.d, self.omega]
}
}
/// The thirteen arguments of the planetary series, in series order.
///
/// The order is `l`, `F`, `D`, `Om`, Mercury through Uranus, Neptune, and the
/// general accumulated precession. Note that `l'` does not appear; the
/// planetary series uses the Moon's mean anomaly, not the Sun's.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PlanetaryArguments {
/// The mean anomaly of the Moon.
pub l: f64,
/// The mean argument of latitude of the Moon.
pub f: f64,
/// The mean elongation of the Moon from the Sun.
pub d: f64,
/// The mean longitude of the Moon's ascending node.
pub omega: f64,
/// The mean longitude of Mercury.
pub mercury: f64,
/// The mean longitude of Venus.
pub venus: f64,
/// The mean longitude of the Earth.
pub earth: f64,
/// The mean longitude of Mars.
pub mars: f64,
/// The mean longitude of Jupiter.
pub jupiter: f64,
/// The mean longitude of Saturn.
pub saturn: f64,
/// The mean longitude of Uranus.
pub uranus: f64,
/// The mean longitude of Neptune.
pub neptune: f64,
/// The general accumulated precession in longitude.
pub precession: f64,
}
impl PlanetaryArguments {
/// Returns the arguments for a time in Julian centuries.
pub fn at(t: f64) -> Self {
Self {
l: mean_anomaly_moon(t),
f: mean_latitude_moon(t),
d: mean_elongation_moon_sun(t),
omega: mean_longitude_node(t),
mercury: mean_longitude_mercury(t),
venus: mean_longitude_venus(t),
earth: mean_longitude_earth(t),
mars: mean_longitude_mars(t),
jupiter: mean_longitude_jupiter(t),
saturn: mean_longitude_saturn(t),
uranus: mean_longitude_uranus(t),
neptune: mean_longitude_neptune(t),
precession: general_precession(t),
}
}
/// Returns the arguments as an array in series order.
pub const fn as_array(self) -> [f64; 13] {
[
self.l,
self.f,
self.d,
self.omega,
self.mercury,
self.venus,
self.earth,
self.mars,
self.jupiter,
self.saturn,
self.uranus,
self.neptune,
self.precession,
]
}
}

105
crates/physics/apsis-frames/src/iau2000/mod.rs Normal file
View File

@ -0,0 +1,105 @@
//! The IAU 2000/2006 Earth orientation models.
//!
//! The transformation from an inertial frame to an Earth-fixed frame is a
//! sequence of rotations, each from a different model:
//!
//! 1. The **frame bias**, a constant rotation from the International
//! Celestial Reference System to the dynamical mean equator and equinox
//! of J2000.0.
//! 2. The **precession**, the secular motion of the pole, from the P03
//! theory adopted in 2006.
//! 3. The **nutation**, the periodic motion of the pole, from the MHB2000
//! series adopted in 2000.
//! 4. The **Earth rotation angle**, which carries the frame from the
//! celestial to the terrestrial meridian.
//! 5. The **polar motion**, the small motion of the pole within the Earth,
//! which a caller supplies because it is published as a daily measured
//! quantity.
//!
//! The composition of the first three is the celestial-to-intermediate
//! matrix, and the composition with the rotation angle and the polar motion
//! is the celestial-to-terrestrial matrix. For satellite work the polar
//! motion is often omitted, because it is under a tenth of an arcsecond,
//! but the nutation and precession are not: the nutation reaches seventeen
//! arcseconds, which is half a kilometer of ground position for a low Earth
//! satellite.
//!
//! All times are Julian centuries from J2000.0 on the Terrestrial Time
//! scale, unless a function says otherwise. The angles a caller supplies to
//! the polar motion are in radians.
//!
//! Source: the models and the order of the rotations follow the IAU SOFA
//! library and its ERFA derivative. The specific functions are named in each
//! module.
mod eect_table;
mod fundamental;
mod nutation;
mod nutation_table;
mod precession;
mod terrestrial;
pub use fundamental::{
DelaunayArguments, PlanetaryArguments, general_precession, mean_anomaly_moon, mean_anomaly_sun,
mean_anomaly_sun_iers2003, mean_elongation_moon_sun, mean_elongation_moon_sun_iers2003,
mean_latitude_moon, mean_longitude_earth, mean_longitude_jupiter, mean_longitude_mars,
mean_longitude_mercury, mean_longitude_neptune, mean_longitude_node, mean_longitude_saturn,
mean_longitude_uranus, mean_longitude_venus,
};
pub use nutation::{Nutation, nutation_2000a, nutation_2000a_with_2006_adjustment};
pub use precession::{
FRAME_BIAS_DECLINATION, FRAME_BIAS_LONGITUDE, FRAME_BIAS_RIGHT_ASCENSION, PrecessionAngles,
bias_precession_matrix, frame_bias_matrix, frame_to_mean_of_date, mean_obliquity,
precession_angles, rotation_z,
};
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,
greenwich_apparent_sidereal_time, greenwich_mean_sidereal_time,
intermediate_to_terrestrial_matrix,
};
/// One term of the luni-solar nutation series.
///
/// The five multipliers are applied to the Delaunay arguments in the order
/// `l`, `l'`, `F`, `D`, `Om`. The six coefficients are the longitude sine,
/// its time derivative, the longitude cosine, the obliquity cosine, its time
/// derivative, and the obliquity sine, in units of 0.1 microarcsecond.
#[derive(Debug, Clone, Copy)]
pub struct SeriesTerm {
/// The integer multipliers of the five Delaunay arguments.
pub argument: [i32; 5],
/// The six coefficients of the term, in the order the reference lists
/// them.
pub coefficients: [f64; 6],
}
/// One term of the complementary terms of the equation of the equinoxes.
///
/// The eight multipliers apply to the arguments `l`, `l'`, `F`, `D`, `Om`,
/// the mean longitude of Venus, the mean longitude of the Earth, and the
/// general accumulated precession. The two coefficients are the sine and
/// cosine amplitudes, in arcseconds.
#[derive(Debug, Clone, Copy)]
pub struct ComplementaryTerm {
/// The integer multipliers of the eight arguments.
pub argument: [i32; 8],
/// The sine and cosine coefficients.
pub coefficients: [f64; 2],
}
/// One term of the planetary nutation series.
///
/// The thirteen multipliers are applied to the planetary arguments in the
/// order `l`, `F`, `D`, `Om`, Mercury through Uranus, Neptune, and the
/// general accumulated precession. The four coefficients are the longitude
/// sine, the longitude cosine, the obliquity sine, and the obliquity cosine,
/// in units of 0.1 microarcsecond.
#[derive(Debug, Clone, Copy)]
pub struct PlanetaryTerm {
/// The integer multipliers of the thirteen planetary arguments.
pub argument: [i32; 13],
/// The four coefficients of the term, in the order the reference lists
/// them.
pub coefficients: [f64; 4],
}

154
crates/physics/apsis-frames/src/iau2000/nutation.rs Normal file
View File

@ -0,0 +1,154 @@
//! The IAU 2000A nutation model.
//!
//! Nutation is the periodic part of the motion of the Earth's rotation axis
//! about its mean position. It is caused mainly by the torque of the Moon
//! and the Sun on the equatorial bulge, with smaller contributions from the
//! planets. The model gives the offset of the true pole from the mean pole
//! as two angles: the nutation in longitude and the nutation in obliquity.
//!
//! The series has two parts whose arguments and coefficient groupings
//! differ. The luni-solar part has 678 terms built from the Delaunay
//! arguments, and its obliquity coefficient pairs the cosine with the
//! constant and time terms while the sine stands alone. The planetary part
//! has 687 terms built from the planetary longitudes, and there the sine
//! pairs with the constant and the cosine stands alone. Getting those two
//! groupings the same way round produces an error of arcseconds, so they are
//! kept separate here.
//!
//! Source: the MHB2000 model of Mathews, Herring and Buffett, adopted by the
//! International Astronomical Union in 2000, as implemented in the IAU SOFA
//! library and its ERFA derivative, functions `eraNut00a` and `eraNut06a`.
use super::fundamental::{DelaunayArguments, PlanetaryArguments};
use super::nutation_table::{LUNI_SOLAR, PLANETARY};
use super::{PlanetaryTerm, SeriesTerm};
use core::f64::consts::PI;
/// Converts a coefficient in units of 0.1 microarcsecond to radians.
const UNIT_TO_RADIAN: f64 = 1.0e-7 * PI / (180.0 * 3_600.0);
/// Converts arcseconds to radians.
const ARCSECOND_TO_RADIAN: f64 = PI / (180.0 * 3_600.0);
/// The nutation in longitude and obliquity, in radians.
#[derive(Debug, Clone, Copy, PartialEq, Default)]
pub struct Nutation {
/// The nutation in longitude, in radians.
pub longitude: f64,
/// The nutation in obliquity, in radians.
pub obliquity: f64,
}
impl Nutation {
/// Creates a nutation from its two angles.
pub const fn new(longitude: f64, obliquity: f64) -> Self {
Self {
longitude,
obliquity,
}
}
/// Returns the nutation in longitude, in arcseconds.
pub fn longitude_arcseconds(self) -> f64 {
self.longitude / ARCSECOND_TO_RADIAN
}
/// Returns the nutation in obliquity, in arcseconds.
pub fn obliquity_arcseconds(self) -> f64 {
self.obliquity / ARCSECOND_TO_RADIAN
}
/// Returns the sum of two nutations.
pub const fn plus(self, other: Self) -> Self {
Self::new(
self.longitude + other.longitude,
self.obliquity + other.obliquity,
)
}
}
/// Evaluates one term of the luni-solar series.
///
/// The five multipliers are applied to the Delaunay arguments. The longitude
/// is the sum of a time dependent sine, a constant sine, and a cosine. The
/// obliquity is the sum of a time dependent cosine, a constant cosine, and a
/// sine.
fn lunisolar_term(term: &SeriesTerm, arguments: &DelaunayArguments, t: f64) -> Nutation {
let multipliers = &term.argument;
let c = &term.coefficients;
let argument = f64::from(multipliers[0]) * arguments.l
+ f64::from(multipliers[1]) * arguments.l_prime
+ f64::from(multipliers[2]) * arguments.f
+ f64::from(multipliers[3]) * arguments.d
+ f64::from(multipliers[4]) * arguments.omega;
let (sine, cosine) = argument.sin_cos();
// The order is the one the reference uses: sp, spt, cp, ce, cet, se.
Nutation::new(
(c[0] + c[1] * t) * sine + c[2] * cosine,
(c[3] + c[4] * t) * cosine + c[5] * sine,
)
}
/// Evaluates one term of the planetary series.
///
/// The thirteen multipliers are applied to the planetary arguments. The
/// longitude is a sine plus a cosine, and the obliquity is a sine plus a
/// cosine, with no time dependence in either.
fn planetary_term(term: &PlanetaryTerm, arguments: &PlanetaryArguments) -> Nutation {
let multipliers = &term.argument;
let c = &term.coefficients;
let values = arguments.as_array();
let mut argument = 0.0;
for index in 0..13 {
argument += f64::from(multipliers[index]) * values[index];
}
let (sine, cosine) = argument.sin_cos();
// The order is the one the reference uses: sp, cp, se, ce.
Nutation::new(c[0] * sine + c[1] * cosine, c[2] * sine + c[3] * cosine)
}
/// Returns the nutation in longitude and obliquity for a time.
///
/// The time is in Julian centuries from J2000.0 on the Terrestrial Time
/// scale. This is the IAU 2000A model, which includes the planetary terms.
///
/// Source: IAU 2000A, as implemented by `eraNut00a`.
pub fn nutation_2000a(t: f64) -> Nutation {
let delaunay = DelaunayArguments::at(t);
let planetary = PlanetaryArguments::at(t);
let mut sum = Nutation::default();
for term in LUNI_SOLAR.iter() {
sum = sum.plus(lunisolar_term(term, &delaunay, t));
}
for term in PLANETARY.iter() {
sum = sum.plus(planetary_term(term, &planetary));
}
Nutation::new(
sum.longitude * UNIT_TO_RADIAN,
sum.obliquity * UNIT_TO_RADIAN,
)
}
/// Returns the nutation with the IAU 2006 precession adjustments applied.
///
/// The IAU 2000A model was built on the IAU 1976 precession, and the IAU
/// 2006 precession changed the precession rates. The adjustment brings the
/// nutation into consistency with the newer precession, which is the
/// combination used for the modern celestial-to-terrestrial transformation.
///
/// Source: IAU 2006, as implemented by `eraNut06a`, which applies the
/// expressions of Wallace and Capitaine (2006), equations 5.
pub fn nutation_2000a_with_2006_adjustment(t: f64) -> Nutation {
let nutation = nutation_2000a(t);
// The adjustment scales the nutation itself rather than adding an
// independent term, which is why it is written as a factor.
let fj2 = -2.777_4e-6 * t;
Nutation::new(
nutation.longitude * (1.0 + 0.469_7e-6 + fj2),
nutation.obliquity * (1.0 + fj2),
)
}

5501
crates/physics/apsis-frames/src/iau2000/nutation_table.rs Normal file
View File

File diff suppressed because it is too large Load Diff

180
crates/physics/apsis-frames/src/iau2000/precession.rs Normal file
View File

@ -0,0 +1,180 @@
//! The IAU 2006 precession model and the frame bias.
//!
//! Precession is the slow, secular motion of the Earth's rotation axis about
//! the pole of the ecliptic, caused mainly by the torque of the Moon and the
//! Sun on the equatorial bulge. The IAU 2006 model expresses it as a
//! rotation from the celestial reference frame to the mean equator and
//! equinox of date, built from four angles.
//!
//! The frame bias accounts for the small offset between the orientation of
//! the International Celestial Reference System and the dynamical mean
//! equator and equinox of J2000.0. It is a constant rotation, and omitting
//! it places a position about twenty milliarcseconds away, which is well
//! above the precision of a modern transformation.
//!
//! Source: the P03 precession theory of Capitaine, Wallace and Chapront,
//! adopted by the International Astronomical Union in 2006, as implemented
//! in the IAU SOFA library and its ERFA derivative, functions `eraPfw06`,
//! `eraPmat06`, `eraFw2m`, `eraBi00` and `eraObl06`.
use apsis_math::Mat3;
use core::f64::consts::PI;
/// Converts arcseconds to radians.
const ARCSECOND_TO_RADIAN: f64 = PI / (180.0 * 3_600.0);
/// The frame bias in right ascension, in radians.
///
/// The offset between the origin of the International Celestial Reference
/// System and the dynamical equinox of J2000.0. The value is
/// minus 0.0146 arcseconds.
pub const FRAME_BIAS_RIGHT_ASCENSION: f64 = -0.014_6 * ARCSECOND_TO_RADIAN;
/// The frame bias in declination, in radians.
///
/// The offset between the pole of the International Celestial Reference
/// System and the mean pole of J2000.0. The value is minus 0.0068192
/// arcseconds.
pub const FRAME_BIAS_DECLINATION: f64 = -0.006_819_2 * ARCSECOND_TO_RADIAN;
/// The frame bias in longitude, in radians.
///
/// The offset between the two origins measured along the equator. The value
/// is minus 0.041775 arcseconds.
pub const FRAME_BIAS_LONGITUDE: f64 = -0.041_775 * ARCSECOND_TO_RADIAN;
/// The four angles that describe the precession from J2000.0 to a date.
///
/// They are applied in the order the reference uses: a rotation by
/// `gamma_bar` about the third axis, then by `phi_bar` about the first, then
/// by minus `psi_bar` about the third, then by minus `epsilon_a` about the
/// first.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PrecessionAngles {
/// The precession of the equator along the ecliptic, in radians.
pub gamma_bar: f64,
/// The inclination of the moving equator to the ecliptic, in radians.
pub phi_bar: f64,
/// The accumulated precession in longitude, in radians.
pub psi_bar: f64,
/// The mean obliquity at the date, in radians.
pub epsilon: f64,
}
/// Returns the mean obliquity of the ecliptic, in radians.
///
/// The polynomial is the IAU 2006 expression, evaluated for a time in Julian
/// centuries from J2000.0 on the Terrestrial Time scale. The constant term
/// is 84381.406 arcseconds, which is the obliquity at J2000.0.
///
/// Source: the IAU 2006 precession, as implemented by `eraObl06`.
pub fn mean_obliquity(t: f64) -> f64 {
let arcseconds = 84_381.406
+ t * (-46.836_769
+ t * (-0.000_183_1
+ t * (0.002_003_40 + t * (-0.000_000_576 + t * (-0.000_000_043_4)))));
arcseconds * ARCSECOND_TO_RADIAN
}
/// Returns the four precession angles for a time.
///
/// The time is in Julian centuries from J2000.0 on the Terrestrial Time
/// scale.
///
/// Source: the P03 precession, as implemented by `eraPfw06`.
pub fn precession_angles(t: f64) -> PrecessionAngles {
let gamma_bar = -0.052_928
+ t * (10.556_378
+ t * (0.493_204_4
+ t * (-0.000_312_38 + t * (-0.000_002_788 + t * (0.000_000_026_0)))));
let phi_bar = 84_381.412_819
+ t * (-46.811_016
+ t * (0.051_126_8
+ t * (0.000_532_89 + t * (-0.000_000_440 + t * (-0.000_000_017_6)))));
let psi_bar = -0.041_775
+ t * (5_038.481_484
+ t * (1.558_417_5
+ t * (-0.000_185_22 + t * (-0.000_026_452 + t * (-0.000_000_014_8)))));
PrecessionAngles {
gamma_bar: gamma_bar * ARCSECOND_TO_RADIAN,
phi_bar: phi_bar * ARCSECOND_TO_RADIAN,
psi_bar: psi_bar * ARCSECOND_TO_RADIAN,
epsilon: mean_obliquity(t),
}
}
/// Returns the frame bias rotation matrix.
///
/// The bias is a constant rotation from the International Celestial
/// Reference System to the dynamical mean equator and equinox of J2000.0.
/// Its three angles are far below a tenth of an arcsecond, but omitting it
/// would leave a systematic error that is larger than the precision of the
/// nutation series.
///
/// Source: as implemented by `eraBi00`, which gives the bias as three small
/// rotations rather than as a matrix.
pub fn frame_bias_matrix() -> Mat3 {
// The reference builds the bias from three elementary rotations, and
// the result is reproduced here in the same order.
rotation_z(FRAME_BIAS_LONGITUDE)
* rotation_x(FRAME_BIAS_DECLINATION)
* rotation_y(FRAME_BIAS_RIGHT_ASCENSION)
}
/// Returns the rotation from the celestial reference frame to the mean
/// equator and equinox of date.
///
/// This is the bias-precession matrix. Nutation is applied on top of it to
/// reach the true equator and equinox of date.
///
/// The time is in Julian centuries from J2000.0 on the Terrestrial Time
/// scale.
///
/// Source: as implemented by `eraPmat06`.
pub fn bias_precession_matrix(t: f64) -> Mat3 {
let angles = precession_angles(t);
frame_to_mean_of_date(&angles)
}
/// Returns the rotation from the celestial reference frame to the mean
/// equator and equinox of date, given precomputed precession angles.
///
/// Source: as implemented by `eraFw2m`, which applies the four elementary
/// rotations in sequence. The reference helper composes each new rotation
/// on the left of the accumulator, so the product here is written with the
/// last rotation applied first, in the same order.
pub fn frame_to_mean_of_date(angles: &PrecessionAngles) -> Mat3 {
rotation_x(-angles.epsilon)
* rotation_z(-angles.psi_bar)
* rotation_x(angles.phi_bar)
* rotation_z(angles.gamma_bar)
}
/// Returns a rotation of `angle` about the first axis, in radians.
///
/// The convention matches the reference helpers `eraRx` and `eraRz`: the
/// rotation is applied to the frame rather than to a vector, so the matrix
/// is the transpose of the one that would rotate a vector by the same angle.
/// Composing `a * b` therefore applies `b` first, and the products in this
/// module read in the order the rotations are physically applied.
pub(crate) fn rotation_x(angle: f64) -> Mat3 {
let (sine, cosine) = angle.sin_cos();
Mat3::from_rows([1.0, 0.0, 0.0], [0.0, cosine, sine], [0.0, -sine, cosine])
}
/// Returns a rotation of `angle` about the second axis, in radians.
///
/// The convention is the one of [`rotation_x`].
fn rotation_y(angle: f64) -> Mat3 {
let (sine, cosine) = angle.sin_cos();
Mat3::from_rows([cosine, 0.0, -sine], [0.0, 1.0, 0.0], [sine, 0.0, cosine])
}
/// Returns a rotation of `angle` about the third axis, in radians.
///
/// The convention is the one of [`rotation_x`].
pub fn rotation_z(angle: f64) -> Mat3 {
let (sine, cosine) = angle.sin_cos();
Mat3::from_rows([cosine, sine, 0.0], [-sine, cosine, 0.0], [0.0, 0.0, 1.0])
}

276
crates/physics/apsis-frames/src/iau2000/terrestrial.rs Normal file
View File

@ -0,0 +1,276 @@
//! The celestial-to-terrestrial transformation.
//!
//! This module composes the frame bias, precession, nutation, Earth rotation
//! angle, and polar motion into the rotation that carries a vector from the
//! inertial frame to the Earth-fixed frame.
//!
//! Two intermediate frames appear in the composition. The **celestial
//! intermediate frame** is the true equator and equinox of date; the
//! rotation from the inertial frame to it is the bias, precession, and
//! nutation. The **terrestrial intermediate frame** is that frame carried to
//! the Earth's prime meridian by the Earth rotation angle. The polar motion
//! is applied last, because it is the motion of the pole within the Earth
//! rather than of the Earth within space.
//!
//! The elementary rotations follow the convention of the reference: `rz`
//! rotates the frame about its third axis, and composing `a * b` applies `b`
//! first. That convention is what makes the products below read in the order
//! the rotations are physically applied.
//!
//! Source: the composition follows the IAU SOFA library and its ERFA
//! derivative. The bias and precession are `eraPmat06`, the nutation matrix
//! is `eraNumat`, the rotation angle is `eraEra00`, the sidereal time for
//! the equinox-based composition is `eraGmst06`, and the polar motion is
//! `eraPom00`, composed by `eraC2t00b`.
use super::nutation::nutation_2000a_with_2006_adjustment;
use super::precession::{frame_to_mean_of_date, mean_obliquity, rotation_z};
use apsis_math::Mat3;
use core::f64::consts::PI;
/// Converts arcseconds to radians.
const ARCSECOND_TO_RADIAN: f64 = PI / (180.0 * 3_600.0);
/// Returns the bias-precession-nutation matrix for a time.
///
/// The matrix carries a vector from the celestial reference frame to the
/// true equator and equinox of date. It is built from the four precession
/// angles with the nutation added to the last two, which is the form the
/// reference uses.
///
/// Source: as implemented by `eraPnm06a` and `eraFw2m`.
pub fn bias_precession_nutation_matrix(t: f64) -> Mat3 {
let angles = super::precession::precession_angles(t);
let nutation = nutation_2000a_with_2006_adjustment(t);
frame_to_mean_of_date(&super::precession::PrecessionAngles {
gamma_bar: angles.gamma_bar,
phi_bar: angles.phi_bar,
psi_bar: angles.psi_bar + nutation.longitude,
epsilon: angles.epsilon + nutation.obliquity,
})
}
/// Returns the complementary terms of the equation of the equinoxes.
///
/// The equation of the equinoxes has, besides the projection of the nutation
/// in longitude onto the true equator, a set of small periodic terms from
/// the same theory. They reach about twice the constant term of the
/// osculating nutation, which is a millisecond of time, so they matter for a
/// sidereal time that is quoted to that precision.
///
/// Source: the MHB2000 complementary terms, as implemented by `eraEect00`.
pub fn equation_of_the_equinoxes_complementary_terms(t: f64) -> f64 {
use super::eect_table::{CONSTANT_TERMS, LINEAR_TERMS};
// The complementary terms use the IERS 2003 forms of two arguments,
// which differ from the ones the nutation series was fitted with.
let arguments = [
super::fundamental::mean_anomaly_moon(t),
super::fundamental::mean_anomaly_sun_iers2003(t),
super::fundamental::mean_latitude_moon(t),
super::fundamental::mean_elongation_moon_sun_iers2003(t),
super::fundamental::mean_longitude_node(t),
super::fundamental::mean_longitude_venus(t),
super::fundamental::mean_longitude_earth(t),
super::fundamental::general_precession(t),
];
let mut constant = 0.0;
for term in CONSTANT_TERMS.iter() {
let argument: f64 = term
.argument
.iter()
.zip(arguments.iter())
.map(|(multiplier, value)| f64::from(*multiplier) * value)
.sum();
let (sine, cosine) = argument.sin_cos();
constant += term.coefficients[0] * sine + term.coefficients[1] * cosine;
}
let mut linear = 0.0;
for term in LINEAR_TERMS.iter() {
let argument: f64 = term
.argument
.iter()
.zip(arguments.iter())
.map(|(multiplier, value)| f64::from(*multiplier) * value)
.sum();
let (sine, cosine) = argument.sin_cos();
linear += term.coefficients[0] * sine + term.coefficients[1] * cosine;
}
// The tabulated coefficients are in arcseconds.
(constant + linear * t) * ARCSECOND_TO_RADIAN
}
/// Returns the equation of the equinoxes for a time.
///
/// The equation of the equinoxes is the difference between the true and mean
/// equinoxes measured along the true equator. It is the quantity a sidereal
/// time must be corrected by to become the hour angle of the true equinox.
///
/// The value is the projection of the nutation in longitude onto the true
/// equator, `dpsi cos(eps)`, plus the complementary terms of the theory.
pub fn equation_of_the_equinoxes(t: f64) -> f64 {
let nutation = nutation_2000a_with_2006_adjustment(t);
let obliquity = mean_obliquity(t);
nutation.longitude * obliquity.cos() + equation_of_the_equinoxes_complementary_terms(t)
}
/// The polar motion of the Earth, as two angles and the rotation.
///
/// The polar motion is the motion of the rotation pole within the Earth. It
/// is measured daily and published by the International Earth Rotation and
/// Reference Systems Service, so it cannot be computed from a model and is
/// supplied by the caller.
#[derive(Debug, Clone, Copy, PartialEq, Default)]
pub struct PolarMotion {
/// The displacement of the pole along the Greenwich meridian, in
/// radians.
pub x: f64,
/// The displacement of the pole along the meridian ninety degrees west,
/// in radians.
pub y: f64,
/// The rotation of the terrestrial frame about the pole, in radians.
///
/// This is the terrestrial intermediate origin locator, which is small
/// and is often taken as zero.
pub s: f64,
}
impl PolarMotion {
/// The zero polar motion, which is the approximation used when the
/// published values are not available.
///
/// The error from omitting the polar motion is under a tenth of an
/// arcsecond, well below the nutation but above the precision of a
/// modern transformation.
pub const NONE: Self = Self {
x: 0.0,
y: 0.0,
s: 0.0,
};
/// Creates a polar motion from its three angles.
pub const fn new(x: f64, y: f64, s: f64) -> Self {
Self { x, y, s }
}
/// Returns `true` when the angles are usable.
pub fn is_valid(self) -> bool {
self.x.is_finite() && self.y.is_finite() && self.s.is_finite()
}
/// Returns the rotation matrix of the polar motion.
///
/// The rotation is a turn about the third axis by the locator, then
/// about the second by minus the `x` displacement, then about the first
/// by minus the `y` displacement. The signs are the ones the reference
/// uses, which carry the terrestrial frame into alignment with the
/// celestial pole.
///
/// Source: as implemented by `eraPom00`, which applies the rotations in
/// sequence by composing each on the left of the accumulator.
pub fn matrix(self) -> Mat3 {
rotation_x(-self.y) * rotation_y(-self.x) * rotation_z(self.s)
}
}
/// Returns the intermediate-to-terrestrial matrix for a time.
///
/// The intermediate frame is the true equator and equinox of date, and the
/// terrestrial frame is the Earth-fixed frame at the same date. This uses
/// the equinox-based approach, in which the rotation from one to the other
/// is Greenwich apparent sidereal time.
///
/// The time arguments are the UT1 and TT Julian dates.
///
/// Source: the composition is the one `eraC2t06a` produces for the equinox
/// based chain, with the polar motion applied separately by the caller of
/// the full transformation.
pub fn intermediate_to_terrestrial_matrix(ut1_julian_date: f64, tt_julian_date: f64) -> Mat3 {
let apparent_sidereal = greenwich_apparent_sidereal_time(ut1_julian_date, tt_julian_date);
rotation_z(apparent_sidereal)
}
/// Returns Greenwich apparent sidereal time for a UT1 and a TT Julian date.
///
/// This is the mean sidereal time plus the equation of the equinoxes, so it
/// is the hour angle of the true equinox.
///
/// Source: as implemented by `eraGst06a`.
pub fn greenwich_apparent_sidereal_time(ut1_julian_date: f64, tt_julian_date: f64) -> f64 {
let t = (tt_julian_date - 2_451_545.0) / 36_525.0;
let mean = greenwich_mean_sidereal_time(ut1_julian_date, tt_julian_date);
(mean + equation_of_the_equinoxes(t)).rem_euclid(2.0 * PI)
}
/// Returns Greenwich mean sidereal time for a UT1 and a TT Julian date.
///
/// The relation is the Earth rotation angle plus the accumulated precession
/// in right ascension, which is the quantity that relates the equinox to the
/// non-rotating origin.
///
/// Source: as implemented by `eraGmst06`.
pub fn greenwich_mean_sidereal_time(ut1_julian_date: f64, tt_julian_date: f64) -> f64 {
let t = (tt_julian_date - 2_451_545.0) / 36_525.0;
// The accumulated precession in right ascension, in arcseconds.
let arcseconds = 0.014_506
+ t * (4_612.156_534
+ t * (1.391_581_7
+ t * (-0.000_000_44 + t * (-0.000_029_956 + t * (-0.000_000_036_8)))));
(earth_rotation_angle(ut1_julian_date) + arcseconds * ARCSECOND_TO_RADIAN).rem_euclid(2.0 * PI)
}
/// Returns the Earth rotation angle for a UT1 Julian date.
///
/// This is the linear expression of the International Astronomical Union,
/// which is exact to the precision the angle is known.
///
/// Source: IERS Conventions (2010), equation 8.32.
pub fn earth_rotation_angle(ut1_julian_date: f64) -> f64 {
let days = ut1_julian_date - 2_451_545.0;
let turns = 0.779_057_273_264_0 + 1.002_737_811_911_354_6 * days;
(turns.rem_euclid(1.0)) * 2.0 * PI
}
/// Returns the celestial-to-terrestrial matrix for a time.
///
/// This composes the whole chain: the frame bias, the precession, the
/// nutation, the Earth rotation, and the polar motion.
///
/// The time arguments are the UT1 and TT Julian dates, and the polar motion
/// is supplied by the caller because it is measured rather than modelled.
///
/// Source: the composition is the one `eraC2t06a` produces.
pub fn celestial_to_terrestrial_matrix(
ut1_julian_date: f64,
tt_julian_date: f64,
polar_motion: PolarMotion,
) -> Mat3 {
let t = (tt_julian_date - 2_451_545.0) / 36_525.0;
let bias_precession_nutation = bias_precession_nutation_matrix(t);
let sidereal = greenwich_apparent_sidereal_time(ut1_julian_date, tt_julian_date);
// The composition is the one `eraC2tcio` produces: the bias, precession,
// and nutation first, then the Earth rotation, then the polar motion,
// with each new rotation composed on the left of the accumulator.
polar_motion.matrix() * rotation_z(sidereal) * bias_precession_nutation
}
/// Returns a rotation of `angle` about the first axis, in radians.
///
/// The convention matches the reference helper `eraRx`: the rotation is
/// applied to the frame, so the matrix is the transpose of the one that
/// would rotate a vector by the same angle.
fn rotation_x(angle: f64) -> Mat3 {
let (sine, cosine) = angle.sin_cos();
Mat3::from_rows([1.0, 0.0, 0.0], [0.0, cosine, sine], [0.0, -sine, cosine])
}
/// Returns a rotation of `angle` about the second axis, in radians.
///
/// The convention is the one of [`rotation_x`].
fn rotation_y(angle: f64) -> Mat3 {
let (sine, cosine) = angle.sin_cos();
Mat3::from_rows([cosine, 0.0, -sine], [0.0, 1.0, 0.0], [sine, 0.0, cosine])
}

1
crates/physics/apsis-frames/src/lib.rs
View File

@ -46,6 +46,7 @@
#![deny(missing_docs)] #![deny(missing_docs)]
mod earth_rotation; mod earth_rotation;
pub mod iau2000;
pub use earth_rotation::{ pub use earth_rotation::{
ERA_J2000_TURNS, ERA_RATE_TURNS_PER_DAY, EarthRotationAngle, rotation_ecef_to_eci, ERA_J2000_TURNS, ERA_RATE_TURNS_PER_DAY, EarthRotationAngle, rotation_ecef_to_eci,

398
crates/physics/apsis-frames/tests/iau2000.rs Normal file
View File

@ -0,0 +1,398 @@
//! Tests for the IAU 2000/2006 Earth orientation models.
//!
//! The expected values are those of the IAU SOFA library and its ERFA
//! derivative for the same instants. They are recorded here so that a change
//! to a coefficient or a rotation convention shows up as a failing
//! assertion.
use apsis_frames::iau2000::{
PolarMotion, bias_precession_nutation_matrix, celestial_to_terrestrial_matrix,
equation_of_the_equinoxes, equation_of_the_equinoxes_complementary_terms, general_precession,
greenwich_apparent_sidereal_time, mean_obliquity, nutation_2000a,
nutation_2000a_with_2006_adjustment, precession_angles,
};
use core::f64::consts::PI;
/// Converts arcseconds to radians.
const ARCSECOND: f64 = PI / (180.0 * 3_600.0);
/// Returns the Terrestrial Time Julian date for a Modified Julian Date on
/// the UTC scale, using the offsets in force from 2017.
fn tt_julian_date(utc_mjd: f64) -> f64 {
utc_mjd + 2_400_000.5 + 69.184 / 86_400.0
}
#[test]
fn the_mean_obliquity_at_j2000_matches_the_reference() {
// The IAU 2006 obliquity at J2000.0 is 84381.406 arcseconds.
let obliquity = mean_obliquity(0.0);
assert!(
(obliquity / ARCSECOND - 84_381.406).abs() < 1e-9,
"{} arcseconds",
obliquity / ARCSECOND
);
}
#[test]
fn the_nutation_at_j2000_matches_the_reference() {
// The IAU 2000A nutation at J2000.0, from ERFA.
let t = 0.0;
let nutation = nutation_2000a(t);
assert!(
(nutation.longitude - -6.754422426417298e-05).abs() < 1e-16,
"longitude {}",
nutation.longitude
);
assert!(
(nutation.obliquity - -2.7970831192374137e-05).abs() < 1e-16,
"obliquity {}",
nutation.obliquity
);
}
#[test]
fn the_nutation_2006_adjustment_is_a_scaling() {
// The IAU 2006 adjustment scales the IAU 2000A nutation rather than
// adding an independent term, so the two must stay in proportion.
for t in [0.0, 0.5, 1.0, 2.0] {
let base = nutation_2000a(t);
let adjusted = nutation_2000a_with_2006_adjustment(t);
let expected_longitude = base.longitude * (1.0 + 0.469_7e-6 - 2.777_4e-6 * t);
let expected_obliquity = base.obliquity * (1.0 - 2.777_4e-6 * t);
assert!(
(adjusted.longitude - expected_longitude).abs() < 1e-20,
"t {t}: {} against {expected_longitude}",
adjusted.longitude
);
assert!(
(adjusted.obliquity - expected_obliquity).abs() < 1e-20,
"t {t}: {} against {expected_obliquity}",
adjusted.obliquity
);
}
}
#[test]
fn the_nutation_matches_the_reference_over_a_century() {
// Recorded from ERFA at each date. The values are in radians.
let cases = [
(40_587.0, 2.1896610891495316e-05, 4.070477585626008e-05),
(51_544.5, -6.754419567422628e-05, -2.797091091919891e-05),
(60_000.0, -4.496357815269203e-05, 3.753534415176348e-05),
];
for (mjd, longitude, obliquity) in cases {
let tt = tt_julian_date(mjd);
let t = (tt - 2_451_545.0) / 36_525.0;
let nutation = nutation_2000a(t);
assert!(
(nutation.longitude - longitude).abs() < 1e-14,
"MJD {mjd}: longitude {} against {longitude}",
nutation.longitude
);
assert!(
(nutation.obliquity - obliquity).abs() < 1e-14,
"MJD {mjd}: obliquity {} against {obliquity}",
nutation.obliquity
);
}
}
#[test]
fn the_precession_angles_at_j2000_match_the_reference() {
// The P03 precession angles at J2000.0.
let angles = precession_angles(0.0);
assert!(
(angles.gamma_bar - -0.052_928 * ARCSECOND).abs() < 1e-18,
"gamma bar {}",
angles.gamma_bar
);
assert!(
(angles.phi_bar - 84_381.412_819 * ARCSECOND).abs() < 1e-15,
"phi bar {}",
angles.phi_bar
);
assert!(
(angles.psi_bar - -0.041_775 * ARCSECOND).abs() < 1e-18,
"psi bar {}",
angles.psi_bar
);
assert!(
(angles.epsilon - 84_381.406 * ARCSECOND).abs() < 1e-15,
"epsilon {}",
angles.epsilon
);
}
#[test]
fn the_bias_precession_nutation_matrix_matches_the_reference() {
// The GCRS-to-true-of-date matrix at J2000.0 on the Terrestrial Time
// scale, from ERFA's `pnm06a`. The values are the first row.
let tt = 51_544.5 + 2_400_000.5 + 69.184 / 86_400.0;
let t = (tt - 2_451_545.0) / 36_525.0;
let matrix = bias_precession_nutation_matrix(t);
let expected = [
[
0.9999999977211409,
6.189934767298694e-05,
2.69478892058386e-05,
],
[
-6.1900102296513e-05,
0.999999997692101,
2.8003132867160385e-05,
],
[
-2.6946155767988307e-05,
-2.8004800880443437e-05,
0.9999999992448179,
],
];
for (row, expected_row) in expected.iter().enumerate() {
for (column, expected_value) in expected_row.iter().enumerate() {
let value = matrix.get(row, column).expect("the indices are in range");
assert!(
(value - expected_value).abs() < 1e-14,
"at ({row}, {column}): {value} against {expected_value}"
);
}
}
}
#[test]
fn the_bias_precession_nutation_matrix_is_a_rotation() {
for mjd in [40_587.0, 51_544.5, 60_000.0] {
let tt = tt_julian_date(mjd);
let t = (tt - 2_451_545.0) / 36_525.0;
let matrix = bias_precession_nutation_matrix(t);
assert!(matrix.is_rotation(1e-12), "MJD {mjd}");
assert!(
(matrix.determinant() - 1.0).abs() < 1e-12,
"MJD {mjd}: determinant {}",
matrix.determinant()
);
}
}
#[test]
fn the_celestial_to_terrestrial_matrix_is_a_rotation() {
for mjd in [40_587.0, 51_544.5, 60_000.0] {
let tt = tt_julian_date(mjd);
let matrix = celestial_to_terrestrial_matrix(tt, tt, PolarMotion::NONE);
assert!(matrix.is_rotation(1e-9), "MJD {mjd}");
assert!(
(matrix.determinant() - 1.0).abs() < 1e-9,
"MJD {mjd}: determinant {}",
matrix.determinant()
);
}
}
#[test]
fn the_celestial_to_terrestrial_matrix_matches_the_reference() {
// The matrix at J2000.0 on the Terrestrial Time scale, from ERFA with
// zero polar motion. Each entry is checked to a nanoradian.
let tt = 51_544.5 + 2_400_000.5 + 69.184 / 86_400.0;
let matrix = celestial_to_terrestrial_matrix(tt, tt, PolarMotion::NONE);
// The first row, from ERFA's `c2t06a` composed with the equinox based
// sidereal time.
let expected_first_row = [
0.18651846184811946,
-0.982451455739228,
-2.248740188309595e-05,
];
for (column, expected_value) in expected_first_row.iter().enumerate() {
let value = matrix.get(0, column).expect("the index is in range");
assert!(
(value - expected_value).abs() < 1e-9,
"at (0, {column}): {value} against {expected_value}"
);
}
}
#[test]
fn the_equation_of_the_equinoxes_includes_the_complementary_terms() {
// The complementary terms are a few thousandths of an arcsecond and must
// not be dropped: their omission would show as a millisecond of
// sidereal time.
let t = 0.0;
let complementary = equation_of_the_equinoxes_complementary_terms(t);
assert!(
complementary.abs() > 1e-9,
"the complementary terms are {complementary}"
);
// At J2000.0 the leading term is 2640.96 microarcseconds times the sine
// of the Moon's ascending node.
let expected = 1.021330096302465e-08;
assert!(
(complementary - expected).abs() < 1e-15,
"{complementary} against {expected}"
);
}
#[test]
fn the_general_precession_is_quadratic_in_time() {
// The argument is (0.024381750 + 0.00000538691 t) t, so at small times
// it grows linearly with the leading coefficient.
for t in [0.1, 0.5, 1.0] {
let value = general_precession(t);
let expected = (0.024_381_750 + 0.000_005_386_91 * t) * t;
assert!((value - expected).abs() < 1e-18, "t {t}");
}
}
#[test]
fn sidereal_time_advances_by_about_a_turn_per_day() {
// Greenwich apparent sidereal time must increase by a full turn plus a
// small amount each day. The small amount is about 3.94 minutes of
// time, which is 0.9856 degrees.
let first = greenwich_apparent_sidereal_time(60_000.0, 60_000.0);
let second = greenwich_apparent_sidereal_time(60_001.0, 60_001.0);
let advance = (second - first).rem_euclid(2.0 * PI);
let degrees = advance.to_degrees();
assert!(
(degrees - 0.985_6).abs() < 0.01,
"advance {degrees} degrees"
);
}
#[test]
fn zero_polar_motion_is_the_identity_rotation() {
let matrix = PolarMotion::NONE.matrix();
for row in 0..3 {
for column in 0..3 {
let value = matrix.get(row, column).expect("the indices are in range");
let expected = if row == column { 1.0 } else { 0.0 };
assert!((value - expected).abs() < 1e-15, "at ({row}, {column})");
}
}
assert!(PolarMotion::NONE.is_valid());
}
#[test]
fn a_polar_motion_offsets_the_rotation() {
// A polar motion of one arcsecond must tilt the matrix by about that
// much, which is 4.85 microradians.
let motion = PolarMotion::new(ARCSECOND, ARCSECOND, 0.0);
let matrix = motion.matrix();
assert!(matrix.is_rotation(1e-12));
// The off-diagonal term is the sine of the displacement.
let value = matrix.get(0, 2).expect("the index is in range");
assert!(
(value - ARCSECOND.sin()).abs() < 1e-15,
"at (0, 2): {value}"
);
assert!(motion.is_valid());
assert!(!PolarMotion::new(f64::NAN, 0.0, 0.0).is_valid());
}
#[test]
fn the_polar_motion_shifts_the_earth_fixed_frame() {
// Applying a polar motion must change the terrestrial position of a
// fixed inertial direction.
let tt = 60_000.0;
let without = celestial_to_terrestrial_matrix(tt, tt, PolarMotion::NONE);
let with = celestial_to_terrestrial_matrix(
tt,
tt,
PolarMotion::new(0.1 * ARCSECOND, -0.2 * ARCSECOND, 0.0),
);
let difference = (without.transpose() * with).transpose() * with;
// The two must differ, and the difference must be of the order of the
// applied angles.
let mut largest = 0.0_f64;
for row in 0..3 {
for column in 0..3 {
let a = without.get(row, column).expect("in range");
let b = difference.get(row, column).expect("in range");
largest = largest.max((a - b).abs());
}
}
assert!(largest > 1e-8, "the largest difference is {largest}");
assert!(largest < 1e-5, "the largest difference is {largest}");
}
#[test]
fn the_nutation_reaches_the_expected_magnitude() {
// The nutation in longitude reaches about seventeen arcseconds and the
// nutation in obliquity about nine, so over a range of dates the
// largest value found must be in that range.
let mut largest_longitude = 0.0_f64;
let mut largest_obliquity = 0.0_f64;
for step in 0..400 {
// A span of about fifty years in steps of a fortnight.
let t = f64::from(step) * 0.125;
let nutation = nutation_2000a(t);
largest_longitude = largest_longitude.max(nutation.longitude.abs());
largest_obliquity = largest_obliquity.max(nutation.obliquity.abs());
}
let longitude_arcseconds = largest_longitude / ARCSECOND;
let obliquity_arcseconds = largest_obliquity / ARCSECOND;
assert!(
(19.0..20.0).contains(&longitude_arcseconds),
"longitude {longitude_arcseconds} arcseconds"
);
assert!(
(8.0..10.0).contains(&obliquity_arcseconds),
"obliquity {obliquity_arcseconds} arcseconds"
);
}
#[test]
fn the_nutation_is_periodic() {
// The dominant terms have periods of about eighteen years and half a
// year, so the value must return close to its initial value after
// eighteen years and must oscillate in between.
let initial = nutation_2000a(0.0).longitude;
let after_eighteen_years = nutation_2000a(0.18).longitude;
assert!(
(after_eighteen_years - initial).abs() < 1.5e-05,
"{} against {initial}",
after_eighteen_years
);
// The sign must change at some point within that span.
let mut seen_positive = false;
let mut seen_negative = false;
for step in 0..200 {
let t = f64::from(step) * 0.001;
let value = nutation_2000a(t).longitude;
seen_positive |= value > 0.0;
seen_negative |= value < 0.0;
}
assert!(seen_positive && seen_negative);
}
#[test]
fn the_equation_of_the_equinoxes_follows_the_nutation() {
// The equation of the equinoxes is dominated by the nutation in
// longitude projected onto the true equator, so it shares its sign.
for t in [0.02, 0.05, 0.1, 0.2] {
let nutation = nutation_2000a_with_2006_adjustment(t);
let equation = equation_of_the_equinoxes(t);
assert!(
nutation.longitude.signum() == equation.signum(),
"t {t}: {} against {}",
nutation.longitude,
equation
);
}
}
#[test]
fn a_matrix_at_two_epochs_separated_by_a_century_differs_by_the_precession() {
// The precession over a century is about 1.4 degrees, so the two
// matrices must differ by an angle of that order.
let first = bias_precession_nutation_matrix(0.0);
let second = bias_precession_nutation_matrix(1.0);
let difference = first.transpose() * second;
// The trace of a rotation matrix is 1 + 2 cos(angle), so the angle
// follows from it.
let cosine = ((difference.trace() - 1.0) / 2.0).clamp(-1.0, 1.0);
let degrees = cosine.acos().to_degrees();
assert!(
(1.0..2.0).contains(&degrees),
"the rotation angle is {degrees} degrees"
);
}

8
validation/README.md
View File

@ -82,6 +82,10 @@ cargo run -p apsis-rs --features rf --example dump_p839 \
# ITU-R P.676-12 精确斜路径（分层积分 + 折射率射线追踪） # ITU-R P.676-12 精确斜路径（分层积分 + 折射率射线追踪）
cargo run -p apsis-rs --features rf --example dump_gas_exact \ cargo run -p apsis-rs --features rf --example dump_gas_exact \
> validation/data/expected/rust_gas_exact.csv > validation/data/expected/rust_gas_exact.csv
# IAU 2000/2006 地球定向（岁差、章动、框架偏移、GCRS→ITRF）
cargo run -p apsis-rs --features physics --example dump_iau2000 \
> validation/data/expected/rust_iau2000.csv
``` ```
Python 侧读取同一批输入，用参考实现重算并输出误差表： Python 侧读取同一批输入，用参考实现重算并输出误差表：
@ -97,6 +101,7 @@ pixi run python validation/python/atmosphere/check_rain_path.py
pixi run python validation/python/link/check_s465.py 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_p839.py
pixi run python validation/python/atmosphere/check_gas_exact.py pixi run python validation/python/atmosphere/check_gas_exact.py
pixi run python validation/python/time/check_iau2000.py
``` ```
两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。 两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。
@ -128,12 +133,13 @@ Earth orientation 数据自动下载；本 crate 涉及的换算不需要该数
| `link/check_s465.py` | 参考天线方向图，3 种口径 × 8 个离轴角 | ITU-R S.465 方程（Python 转写） | 4.0e-10 dB | | `link/check_s465.py` | 参考天线方向图，3 种口径 × 8 个离轴角 | ITU-R S.465 方程（Python 转写） | 4.0e-10 dB |
| `atmosphere/check_p839.py` | 零度等温线高度，15 个全球点 | `itur` 的 ITU-R P.839-4 实现 | 3.3e-13 km | | `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（相对） | | `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 |
### 关于射频层对拍的说明 ### 关于射频层对拍的说明
射频层分两类：有权威参考实现的模型与纯闭式算式。 射频层分两类：有权威参考实现的模型与纯闭式算式。
- 有参考实现的，直接对拍：大气与雨衰用 `itur`，它是 ITU-R 研究组维护的建议书实现； - 有参考实现的，直接对拍：地球定向用 `erfa`（IAU SOFA 的派生实现），大气与雨衰用 `itur`；
链路预算的每一步都是初等运算，`check_link.py` 用 Python 独立重写同一方程。 链路预算的每一步都是初等运算，`check_link.py` 用 Python 独立重写同一方程。
- 其余以解析性质断言为主：各 crate 的单元测试验证模型必须满足的性质，例如 3 dB 点、 - 其余以解析性质断言为主：各 crate 的单元测试验证模型必须满足的性质，例如 3 dB 点、
平方律下降、Friis 级联行为、P.838-3 系数在表列频率处的取值。 平方律下降、Friis 级联行为、P.838-3 系数在表列频率处的取值。

156
validation/python/time/check_iau2000.py Normal file
View File

@ -0,0 +1,156 @@
"""Cross-check the ApsisRS IAU 2000/2006 Earth orientation against ERFA.
Reads the quantities written by ``cargo run -p apsis-rs --features physics
--example dump_iau2000`` and recomputes each one with ERFA, which is the
reference implementation distributed by the International Astronomical Union
through the SOFA library and its ERFA derivative.
The script needs no network access.
Usage:
pixi run python validation/python/time/check_iau2000.py [path]
Exit status is zero when every quantity agrees within tolerance and one
otherwise.
"""
from __future__ import annotations
import csv
import sys
from pathlib import Path
import erfa
# Tolerances.
#
# The nutation, obliquity, and precession are reproduced at the level of the
# last bits of a double, so they are compared at 1e-14.
#
# The sidereal time and the matrices accumulate the double precision error of
# several rotations, so they are allowed more.
#
# Greenwich apparent sidereal time gets its own tolerance because ERFA's
# `gst06a` is built on the non-rotating origin while the Rust value is built
# on the equinox. The two formulations agree to about two microarcseconds,
# which is 1e-11 radians, and that is what is allowed here.
SCALAR_TOLERANCE = 1e-14
SIDEREAL_TOLERANCE = 1e-10
MATRIX_TOLERANCE = 1e-9
def tolerance_for(quantity: str) -> float:
"""Returns the tolerance that applies to a quantity."""
if quantity == "gast":
return SIDEREAL_TOLERANCE
if quantity in ("bpn", "c2t"):
return MATRIX_TOLERANCE
return SCALAR_TOLERANCE
DEFAULT_INPUT = Path("validation/data/expected/rust_iau2000.csv")
# The Terrestrial Time offset from UTC, in seconds, matching the Rust side.
TT_MINUS_UTC_SECONDS = 32.184 + 37.0
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 tt_julian_date(utc_mjd: float) -> float:
"""Returns the Terrestrial Time Julian date for a UTC Modified Julian Date."""
return utc_mjd + 2400000.5 + TT_MINUS_UTC_SECONDS / 86400.0
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_iau2000 > " + 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
worst: dict[str, float] = {}
worst_at: dict[str, object] = {}
compared = 0
for row in rows:
mjd = float(row["mjd"])
quantity = row["quantity"]
index = int(row["index"])
actual = float(row["value"])
tt = tt_julian_date(mjd)
d1 = 2400000.5
d2 = tt - d1
if quantity == "obliquity":
reference = erfa.obl06(d1, d2)
elif quantity == "nutation_longitude_a":
reference = erfa.nut00a(d1, d2)[0]
elif quantity == "nutation_obliquity_a":
reference = erfa.nut00a(d1, d2)[1]
elif quantity == "nutation_longitude_06a":
reference = erfa.nut06a(d1, d2)[0]
elif quantity == "nutation_obliquity_06a":
reference = erfa.nut06a(d1, d2)[1]
elif quantity == "precession_gamma_bar":
reference = erfa.pfw06(d1, d2)[0]
elif quantity == "precession_phi_bar":
reference = erfa.pfw06(d1, d2)[1]
elif quantity == "precession_psi_bar":
reference = erfa.pfw06(d1, d2)[2]
elif quantity == "precession_epsilon":
reference = erfa.pfw06(d1, d2)[3]
elif quantity == "bpn":
matrix = erfa.pnm06a(d1, d2)
reference = matrix[index // 3][index % 3]
elif quantity == "gast":
reference = erfa.gst06a(d1, d2, d1, d2)
elif quantity == "c2t":
# The equinox based chain: the bias, precession and nutation,
# then Greenwich apparent sidereal time, then the polar motion.
matrix = erfa.c2t06a(d1, d2, d1, d2, 0.0, 0.0)
reference = matrix[index // 3][index % 3]
else:
continue
error = abs(actual - reference)
if error > worst.get(quantity, 0.0):
worst[quantity] = error
worst_at[quantity] = (mjd, index)
compared += 1
print(f"values compared: {compared}")
for quantity in sorted(worst):
tolerance = tolerance_for(quantity)
status = "ok" if worst[quantity] <= tolerance else "FAIL"
print(
f" {quantity:26} worst {worst[quantity]:.3e} "
f"(tolerance {tolerance:.0e}) at {worst_at[quantity]} {status}"
)
passed = all(
value <= tolerance_for(quantity) for quantity, value in worst.items()
)
if passed and compared > 0:
print("result: pass")
return 0
print("result: FAIL")
return 1
if __name__ == "__main__":
raise SystemExit(main(sys.argv))

98
validation/tools/generate_eect_table.py Normal file
View File

@ -0,0 +1,98 @@
#!/usr/bin/env python3
"""Generate the Rust source for the IAU 2000 complementary terms.
The complementary terms are the small periodic additions to the equation of
the equinoxes that the nutation in longitude does not contain. They come
from the same theory as the nutation and are needed for the equation of the
equinoxes to be consistent with the nutation series.
The reference is the ERFA library, a derivative of the IAU SOFA library.
Fetch the reference with:
curl -O https://raw.githubusercontent.com/liberfa/erfa/master/src/eect00.c
Usage:
python validation/tools/generate_eect_table.py path/to/eect00.c
"""
from __future__ import annotations
import re
import sys
from pathlib import Path
def extract_array(source: str, name: str) -> str:
"""Returns the body of a C array initializer."""
start = source.index(name + "[] = {")
brace = source.index("{", start)
depth = 0
index = brace
while index < len(source):
if source[index] == "{":
depth += 1
elif source[index] == "}":
depth -= 1
if depth == 0:
break
index += 1
return source[brace : index + 1]
def parse_terms(body: str) -> list[tuple[str, str]]:
"""Returns the multiplier and coefficient groups of each term."""
body = re.sub(r"/\*.*?\*/", "", body, flags=re.S)
return re.findall(r"\{\{([^{}]*)\},\s*([^}]*)\}", body)
def emit(name: str, doc: str, terms: list[tuple[str, str]]) -> None:
"""Prints one Rust constant."""
print(doc)
print(f"pub static {name}: [ComplementaryTerm; {len(terms)}] = [")
for multipliers, coefficients in terms:
arguments = ", ".join(value.strip() for value in multipliers.split(","))
values = ", ".join(value.strip() for value in coefficients.split(","))
print(
f" ComplementaryTerm {{ argument: [{arguments}], "
f"coefficients: [{values}] }},"
)
print("];")
print()
def main(argv: list[str]) -> int:
if len(argv) < 2:
print(__doc__, file=sys.stderr)
return 2
source = Path(argv[1]).read_text()
constant = parse_terms(extract_array(source, "e0"))
linear = parse_terms(extract_array(source, "e1"))
print("//! The IAU 2000 complementary terms of the equation of the equinoxes.")
print("//!")
print("//! The table is generated from the reference implementation by")
print("//! `validation/tools/generate_eect_table.py` and is reproduced without")
print("//! change.")
print("//!")
print("//! Source: the MHB2000 model, as distributed with the IAU SOFA library")
print("//! and implemented by ERFA.")
print()
print("use super::ComplementaryTerm;")
print()
emit(
"CONSTANT_TERMS",
"/// The terms that do not depend on time.\n///\n/// The eight multipliers apply to the arguments `l`, `l'`, `F`, `D`, `Om`,\n/// the mean longitude of Venus, the mean longitude of the Earth, and the\n/// general accumulated precession, in that order.",
constant,
)
emit(
"LINEAR_TERMS",
"/// The terms whose amplitude grows linearly with time.\n///\n/// The arguments are the same as for the constant terms.",
linear,
)
return 0
if __name__ == "__main__":
raise SystemExit(main(sys.argv))

143
validation/tools/generate_nutation_table.py Normal file
View File

@ -0,0 +1,143 @@
#!/usr/bin/env python3
"""Generate the Rust source for the IAU 2000A nutation series.
The series is the MHB2000 model. It has a luni-solar part and a planetary
part, both long tables of integer argument multipliers and real
coefficients, so they are transcribed mechanically from the reference
rather than by hand.
The reference is the ERFA library, a derivative of the IAU SOFA library
maintained by the International Astronomical Union. Its source is the
authoritative statement of the series.
Field order is the one the reference uses, which is not the order the
coefficients are listed in the paper:
luni-solar: multipliers l, l', F, D, Om; then sp, spt, cp, ce, cet, se
planetary: multipliers l, F, D, Om, Me, Ve, E, Ma, Ju, Sa, Ur, Ne, pa;
then sp, cp, se, ce
The two series use different argument sets and different coefficient
groupings for the obliquity, so the generated Rust types keep them separate.
Fetch the reference with:
curl -O https://raw.githubusercontent.com/liberfa/erfa/master/src/nut00a.c
Usage:
python validation/tools/generate_nutation_table.py path/to/nut00a.c
"""
from __future__ import annotations
import re
import sys
from pathlib import Path
def extract_array(source: str, name: str) -> str:
"""Returns the body of a C array initializer."""
start = source.index(name + "[] = {")
brace = source.index("{", start)
depth = 0
index = brace
while index < len(source):
if source[index] == "{":
depth += 1
elif source[index] == "}":
depth -= 1
if depth == 0:
break
index += 1
return source[brace : index + 1]
def parse_terms(body: str) -> list[list[str]]:
"""Returns the numeric terms of an array of brace groups."""
body = re.sub(r"/\*.*?\*/", "", body, flags=re.S)
terms = []
for group in re.findall(r"\{([^{}]*)\}", body):
values = group.replace("\n", " ").strip()
if values:
terms.append([value.strip() for value in values.split(",")])
return terms
def main(argv: list[str]) -> int:
if len(argv) < 2:
print(__doc__, file=sys.stderr)
return 2
source = Path(argv[1]).read_text()
lunisolar = parse_terms(extract_array(source, "xls"))
planetary = parse_terms(extract_array(source, "xpl"))
if any(len(term) != 11 for term in lunisolar):
print("unexpected luni-solar row width", file=sys.stderr)
return 1
if any(len(term) != 17 for term in planetary):
print("unexpected planetary row width", file=sys.stderr)
return 1
print("//! The IAU 2000A nutation series.")
print("//!")
print("//! The table is generated from the reference implementation of the")
print("//! model by `validation/tools/generate_nutation_table.py` and is")
print("//! reproduced without change.")
print("//!")
print("//! The field order follows the reference, which is not the order the")
print("//! coefficients are listed in the paper.")
print("//!")
print("//! Source: the MHB2000 model of Mathews, Herring and Buffett, adopted by")
print("//! the International Astronomical Union in 2000 and distributed as part")
print("//! of the IAU SOFA library, as implemented by ERFA.")
print()
print("use super::{PlanetaryTerm, SeriesTerm};")
print()
print("/// The luni-solar terms of the series.")
print("///")
print("/// The multipliers apply to the arguments `l`, `l'`, `F`, `D`, and `Om`")
print("/// in that order. The coefficients are the longitude sine, its time")
print("/// derivative, the longitude cosine, the obliquity cosine, its time")
print("/// derivative, and the obliquity sine, in units of 0.1 microarcsecond.")
print(f"pub static LUNI_SOLAR: [SeriesTerm; {len(lunisolar)}] = [")
for term in lunisolar:
multipliers = ", ".join(term[:5])
coefficients = ", ".join(
value if "." in value or "e" in value.lower() else value + ".0"
for value in term[5:]
)
print(
f" SeriesTerm {{ argument: [{multipliers}], "
f"coefficients: [{coefficients}] }},"
)
print("];")
print()
print("/// The planetary terms of the series.")
print("///")
print("/// The multipliers apply to the arguments `l`, `F`, `D`, `Om`, the mean")
print("/// longitudes of Mercury through Uranus, the mean longitude of Neptune,")
print("/// and the general accumulated precession, in that order. The")
print("/// coefficients are the longitude sine, the longitude cosine, the")
print("/// obliquity sine, and the obliquity cosine, in units of 0.1")
print("/// microarcsecond.")
print(f"pub static PLANETARY: [PlanetaryTerm; {len(planetary)}] = [")
for term in planetary:
multipliers = ", ".join(term[:13])
coefficients = ", ".join(
value if "." in value or "e" in value.lower() else value + ".0"
for value in term[13:]
)
print(
f" PlanetaryTerm {{ argument: [{multipliers}], "
f"coefficients: [{coefficients}] }},"
)
print("];")
return 0
if __name__ == "__main__":
raise SystemExit(main(sys.argv))