From 8b5213747aa11e9e29000158d927b2996af894e4 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 13:15:52 +0800 Subject: [PATCH] physics: add leap second table and UTC TAI conversion --- crates/physics/apsis-leap/src/lib.rs | 46 ++++- crates/physics/apsis-leap/src/table.rs | 241 +++++++++++++++++++++++ crates/physics/apsis-leap/src/time.rs | 108 ++++++++++ crates/physics/apsis-leap/tests/table.rs | 80 ++++++++ crates/physics/apsis-leap/tests/time.rs | 117 +++++++++++ 5 files changed, 591 insertions(+), 1 deletion(-) create mode 100644 crates/physics/apsis-leap/src/table.rs create mode 100644 crates/physics/apsis-leap/src/time.rs create mode 100644 crates/physics/apsis-leap/tests/table.rs create mode 100644 crates/physics/apsis-leap/tests/time.rs diff --git a/crates/physics/apsis-leap/src/lib.rs b/crates/physics/apsis-leap/src/lib.rs index 035efc7..7af03c1 100644 --- a/crates/physics/apsis-leap/src/lib.rs +++ b/crates/physics/apsis-leap/src/lib.rs @@ -1 +1,45 @@ -//! Embedded leap second table and UTC/TAI time scale conversion. +//! Embedded leap second table for the UTC time scale. +//! +//! Coordinated Universal Time (UTC) is kept within 0.9 s of UT1 by +//! inserting or removing leap seconds. Because UTC is not a uniform time +//! scale, converting between UTC and the uniform scales TAI, TT, and GPS +//! requires the offset `TAI - UTC` in force at the instant of interest. +//! +//! This crate embeds that offset as a compile time table and answers +//! queries against it. Dates are handled as Modified Julian Date (MJD), +//! which is the Julian Date minus 2 400 000.5. An MJD is a count of days +//! since 1858-11-17T00:00:00, so it carries no time scale of its own; the +//! caller states the scale by choosing which function to call. +//! +//! The table begins at 1972-01-01, when UTC adopted the current integer +//! second offset scheme. Offsets before that date were fractional and +//! followed a different rule, so queries earlier than the first entry +//! return `None` rather than a misleading value. +//! +//! Source: IERS Bulletin C, which announces leap seconds after +//! confirmation by the International Earth Rotation and Reference Systems +//! Service. The most recent entry in this table is the leap second of +//! 2017-01-01; no leap second has been announced since. +//! +//! ``` +//! use apsis_leap::{tai_minus_utc_at_mjd, LEAP_SECONDS}; +//! +//! // The last entry in the embedded table fixes the current offset. +//! let last = LEAP_SECONDS.last().expect("table is not empty"); +//! assert_eq!(last.tai_minus_utc, 37.0); +//! assert_eq!(tai_minus_utc_at_mjd(60_000.0), Some(37.0)); +//! ``` + +#![forbid(unsafe_code)] +#![deny(missing_docs)] + +mod table; +mod time; + +pub use table::{ + LEAP_SECONDS, LeapSecond, TABLE_COVERS_FROM_MJD, TAI_MINUS_UTC_AFTER_LAST_ENTRY, mjd_from_ymd, +}; +pub use time::{ + MJD_AT_UNIX_EPOCH, calendar_mjd, tai_minus_utc_at_mjd, tai_mjd_from_utc_mjd, + utc_mjd_from_tai_mjd, +}; diff --git a/crates/physics/apsis-leap/src/table.rs b/crates/physics/apsis-leap/src/table.rs new file mode 100644 index 0000000..8d385c4 --- /dev/null +++ b/crates/physics/apsis-leap/src/table.rs @@ -0,0 +1,241 @@ +//! The compile time leap second table and the calendar helpers that place +//! its entries on the Modified Julian Date scale. + +/// One leap second announcement. +/// +/// The entry records the date on which a new value of `TAI - UTC` took +/// effect. It does not record the size of the inserted second, because in +/// every announcement to date the inserted second was exactly one second +/// and the new offset is one greater than the previous one. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct LeapSecond { + /// Year of the effective date, in the proleptic Gregorian calendar. + pub year: i32, + /// Month of the effective date, in `1..=12`. + pub month: u32, + /// Day of the effective date, in `1..=31`. + pub day: u32, + /// The value of `TAI - UTC` in seconds from this date onward. + pub tai_minus_utc: f64, +} + +impl LeapSecond { + /// Returns the Modified Julian Date of the effective instant. + /// + /// Leap seconds take effect at 00:00:00 UTC, so the effective instant + /// is the start of the given calendar day. + pub const fn effective_mjd(self) -> f64 { + mjd_from_ymd(self.year, self.month, self.day) + } +} + +/// The leap second table, ordered by effective date. +/// +/// Source: IERS Bulletin C. The first entry is 1972-01-01, when the +/// current scheme began with an initial offset of 10 s inherited from +/// earlier adjustments. +pub const LEAP_SECONDS: &[LeapSecond] = &[ + LeapSecond { + year: 1972, + month: 1, + day: 1, + tai_minus_utc: 10.0, + }, + LeapSecond { + year: 1972, + month: 7, + day: 1, + tai_minus_utc: 11.0, + }, + LeapSecond { + year: 1973, + month: 1, + day: 1, + tai_minus_utc: 12.0, + }, + LeapSecond { + year: 1974, + month: 1, + day: 1, + tai_minus_utc: 13.0, + }, + LeapSecond { + year: 1975, + month: 1, + day: 1, + tai_minus_utc: 14.0, + }, + LeapSecond { + year: 1976, + month: 1, + day: 1, + tai_minus_utc: 15.0, + }, + LeapSecond { + year: 1977, + month: 1, + day: 1, + tai_minus_utc: 16.0, + }, + LeapSecond { + year: 1978, + month: 1, + day: 1, + tai_minus_utc: 17.0, + }, + LeapSecond { + year: 1979, + month: 1, + day: 1, + tai_minus_utc: 18.0, + }, + LeapSecond { + year: 1980, + month: 1, + day: 1, + tai_minus_utc: 19.0, + }, + LeapSecond { + year: 1981, + month: 7, + day: 1, + tai_minus_utc: 20.0, + }, + LeapSecond { + year: 1982, + month: 7, + day: 1, + tai_minus_utc: 21.0, + }, + LeapSecond { + year: 1983, + month: 7, + day: 1, + tai_minus_utc: 22.0, + }, + LeapSecond { + year: 1985, + month: 7, + day: 1, + tai_minus_utc: 23.0, + }, + LeapSecond { + year: 1988, + month: 1, + day: 1, + tai_minus_utc: 24.0, + }, + LeapSecond { + year: 1990, + month: 1, + day: 1, + tai_minus_utc: 25.0, + }, + LeapSecond { + year: 1991, + month: 1, + day: 1, + tai_minus_utc: 26.0, + }, + LeapSecond { + year: 1992, + month: 7, + day: 1, + tai_minus_utc: 27.0, + }, + LeapSecond { + year: 1993, + month: 7, + day: 1, + tai_minus_utc: 28.0, + }, + LeapSecond { + year: 1994, + month: 7, + day: 1, + tai_minus_utc: 29.0, + }, + LeapSecond { + year: 1996, + month: 1, + day: 1, + tai_minus_utc: 30.0, + }, + LeapSecond { + year: 1997, + month: 7, + day: 1, + tai_minus_utc: 31.0, + }, + LeapSecond { + year: 1999, + month: 1, + day: 1, + tai_minus_utc: 32.0, + }, + LeapSecond { + year: 2006, + month: 1, + day: 1, + tai_minus_utc: 33.0, + }, + LeapSecond { + year: 2009, + month: 1, + day: 1, + tai_minus_utc: 34.0, + }, + LeapSecond { + year: 2012, + month: 7, + day: 1, + tai_minus_utc: 35.0, + }, + LeapSecond { + year: 2015, + month: 7, + day: 1, + tai_minus_utc: 36.0, + }, + LeapSecond { + year: 2017, + month: 1, + day: 1, + tai_minus_utc: 37.0, + }, +]; + +/// The Modified Julian Date of the first table entry, 1972-01-01. +/// +/// Queries earlier than this instant are rejected because UTC did not yet +/// use integer second offsets. +pub const TABLE_COVERS_FROM_MJD: f64 = mjd_from_ymd( + LEAP_SECONDS[0].year, + LEAP_SECONDS[0].month, + LEAP_SECONDS[0].day, +); + +/// The value of `TAI - UTC` on and after the last table entry. +/// +/// This is the offset that remains in force until a new leap second is +/// announced. +pub const TAI_MINUS_UTC_AFTER_LAST_ENTRY: f64 = LEAP_SECONDS[LEAP_SECONDS.len() - 1].tai_minus_utc; + +/// Converts a proleptic Gregorian calendar date to a Modified Julian Date. +/// +/// The result is the MJD of 00:00:00 on the given day. The algorithm is the +/// Fliegel-Van Flandern integer conversion, which is exact for every date +/// the Gregorian calendar is defined on and does not need a lookup table. +/// +/// The computation is done in `i64` and converted once at the end, so the +/// result is exact for the range of dates in which a leap second can +/// plausibly occur. +pub const fn mjd_from_ymd(year: i32, month: u32, day: u32) -> f64 { + let a = (14 - month as i64) / 12; + let y = year as i64 + 4800 - a; + let m = month as i64 + 12 * a - 3; + let jdn = day as i64 + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045; + // The Julian Day Number counts days from noon, and MJD 0 is the + // midnight starting 1858-11-17, whose Julian Day Number is 2 400 001. + (jdn - 2_400_001) as f64 +} diff --git a/crates/physics/apsis-leap/src/time.rs b/crates/physics/apsis-leap/src/time.rs new file mode 100644 index 0000000..d5d3a58 --- /dev/null +++ b/crates/physics/apsis-leap/src/time.rs @@ -0,0 +1,108 @@ +//! Conversions between UTC and TAI expressed on the Modified Julian Date +//! scale. + +use crate::table::{ + LEAP_SECONDS, TABLE_COVERS_FROM_MJD, TAI_MINUS_UTC_AFTER_LAST_ENTRY, mjd_from_ymd, +}; + +/// The Modified Julian Date of the Unix epoch, 1970-01-01T00:00:00. +/// +/// MJD carries no time scale of its own, so this constant is the same +/// number whether the instant is read in UTC, TAI, or any other scale. +pub const MJD_AT_UNIX_EPOCH: f64 = 40_587.0; + +/// The number of seconds in a day. +const SECONDS_PER_DAY: f64 = 86_400.0; + +/// Returns the value of `TAI - UTC` in seconds at the given UTC Modified +/// Julian Date. +/// +/// The offset is a step function: it changes at 00:00:00 UTC on the +/// effective date of each leap second and is constant in between. A query +/// exactly on an effective date returns the new value. +/// +/// Returns `None` when the instant precedes 1972-01-01, before which UTC +/// used fractional offsets and this table does not apply. +pub fn tai_minus_utc_at_mjd(utc_mjd: f64) -> Option { + if !utc_mjd.is_finite() || utc_mjd < TABLE_COVERS_FROM_MJD { + return None; + } + let mut offset = TAI_MINUS_UTC_AFTER_LAST_ENTRY; + for entry in LEAP_SECONDS { + if utc_mjd >= entry.effective_mjd() { + offset = entry.tai_minus_utc; + } else { + break; + } + } + Some(offset) +} + +/// Converts a UTC Modified Julian Date to the corresponding TAI Modified +/// Julian Date. +/// +/// The conversion adds the offset in force at that instant. Because the +/// offset is a step function evaluated on the UTC side, the conversion is +/// defined for every UTC instant; the inverse conversion is not, because a +/// leap second makes some TAI instants correspond to no UTC instant. +/// +/// Returns `None` when the instant precedes 1972-01-01. +pub fn tai_mjd_from_utc_mjd(utc_mjd: f64) -> Option { + let offset = tai_minus_utc_at_mjd(utc_mjd)?; + Some(utc_mjd + offset / SECONDS_PER_DAY) +} + +/// Converts a TAI Modified Julian Date to the corresponding UTC Modified +/// Julian Date. +/// +/// The offset that applies depends on the UTC instant being sought, so the +/// value is recovered by a monotone iteration. The iteration starts from +/// the smallest possible answer, obtained with the largest offset in the +/// table, and increases until the offset stops changing. Because the offset +/// is a non-decreasing step function of UTC, the sequence is non-decreasing +/// and reaches the fixed point in as many steps as there are offset changes +/// between the seed and the answer, which is a small number. +/// +/// Returns `None` when the instant precedes the start of the table, or when +/// the result would fall before that boundary. +pub fn utc_mjd_from_tai_mjd(tai_mjd: f64) -> Option { + if !tai_mjd.is_finite() { + return None; + } + // Seed with the largest offset, which gives the smallest candidate UTC. + let mut utc_mjd = tai_mjd - TAI_MINUS_UTC_AFTER_LAST_ENTRY / SECONDS_PER_DAY; + for _ in 0..64 { + let offset = offset_for_iteration(utc_mjd); + let next = tai_mjd - offset / SECONDS_PER_DAY; + if next == utc_mjd { + break; + } + utc_mjd = next; + } + if utc_mjd < TABLE_COVERS_FROM_MJD { + return None; + } + Some(utc_mjd) +} + +/// Returns the offset in force near `utc_mjd`, using the first table entry +/// for instants below the start of the table. +/// +/// This is used only while iterating in [`utc_mjd_from_tai_mjd`], where a +/// candidate may transiently precede the boundary even though the answer +/// does not. Callers that need to distinguish an out of range instant from +/// a valid one must use [`tai_minus_utc_at_mjd`]. +fn offset_for_iteration(utc_mjd: f64) -> f64 { + tai_minus_utc_at_mjd(utc_mjd).unwrap_or(LEAP_SECONDS[0].tai_minus_utc) +} + +/// Converts a proleptic Gregorian calendar date at 00:00:00 to a Modified +/// Julian Date. +/// +/// This is a convenience re-export of the calendar conversion used to place +/// the leap second table entries. It is provided so that callers do not +/// have to reproduce the conversion, and because the calibration scripts in +/// `validation/` need the same function on both sides of a comparison. +pub const fn calendar_mjd(year: i32, month: u32, day: u32) -> f64 { + mjd_from_ymd(year, month, day) +} diff --git a/crates/physics/apsis-leap/tests/table.rs b/crates/physics/apsis-leap/tests/table.rs new file mode 100644 index 0000000..b21f2e7 --- /dev/null +++ b/crates/physics/apsis-leap/tests/table.rs @@ -0,0 +1,80 @@ +//! Tests for the leap second table. + +use apsis_leap::{ + LEAP_SECONDS, TABLE_COVERS_FROM_MJD, TAI_MINUS_UTC_AFTER_LAST_ENTRY, mjd_from_ymd, +}; + +#[test] +fn table_is_ordered_and_monotonic() { + let mut previous = f64::NEG_INFINITY; + for entry in LEAP_SECONDS { + assert!( + (1..=12).contains(&entry.month), + "month out of range: {entry:?}" + ); + assert!((1..=31).contains(&entry.day), "day out of range: {entry:?}"); + let mjd = entry.effective_mjd(); + assert!(mjd > previous, "table is not ordered at {entry:?}"); + previous = mjd; + } +} + +#[test] +fn every_entry_increases_the_offset_by_one_second() { + // The table starts at the inherited offset of ten seconds and every + // later announcement inserted exactly one second. + assert_eq!(LEAP_SECONDS[0].tai_minus_utc, 10.0); + for pair in LEAP_SECONDS.windows(2) { + let (earlier, later) = (pair[0], pair[1]); + assert!( + (later.tai_minus_utc - earlier.tai_minus_utc - 1.0).abs() < f64::EPSILON, + "offset step is not one second between {earlier:?} and {later:?}" + ); + } +} + +#[test] +fn table_extremes_match_the_announced_values() { + assert_eq!(LEAP_SECONDS[0].effective_mjd(), 41_317.0); + let last = LEAP_SECONDS[LEAP_SECONDS.len() - 1]; + assert_eq!(last.year, 2017); + assert_eq!(last.month, 1); + assert_eq!(last.day, 1); + assert_eq!(last.tai_minus_utc, 37.0); + assert_eq!(last.effective_mjd(), 57_754.0); +} + +#[test] +fn table_contains_every_announcement_since_1972() { + // Twenty seven leap seconds were introduced between 1972 and 2017 + // inclusive of the initial 1972-01-01 entry. + assert_eq!(LEAP_SECONDS.len(), 28); +} + +#[test] +fn coverage_constants_match_the_table_boundaries() { + assert_eq!(TABLE_COVERS_FROM_MJD, LEAP_SECONDS[0].effective_mjd()); + assert_eq!(TAI_MINUS_UTC_AFTER_LAST_ENTRY, 37.0); +} + +#[test] +fn known_calendar_dates_map_to_known_modified_julian_dates() { + // MJD 0 is the start of 1858-11-17. + assert_eq!(mjd_from_ymd(1858, 11, 17), 0.0); + // The calendar day of J2000.0 is MJD 51544. + assert_eq!(mjd_from_ymd(2000, 1, 1), 51_544.0); + // The Unix epoch is MJD 40587. + assert_eq!(mjd_from_ymd(1970, 1, 1), 40_587.0); + // The last leap second. + assert_eq!(mjd_from_ymd(2017, 1, 1), 57_754.0); + assert_eq!(mjd_from_ymd(2024, 1, 1), 60_310.0); +} + +#[test] +fn calendar_conversion_handles_leap_days_and_month_ends() { + // 2000 is a leap year, 1900 is not. + assert_eq!(mjd_from_ymd(2000, 3, 1) - mjd_from_ymd(2000, 2, 1), 29.0); + assert_eq!(mjd_from_ymd(1900, 3, 1) - mjd_from_ymd(1900, 2, 1), 28.0); + // Consecutive month ends are one day apart. + assert_eq!(mjd_from_ymd(2024, 3, 1) - mjd_from_ymd(2024, 2, 29), 1.0); +} diff --git a/crates/physics/apsis-leap/tests/time.rs b/crates/physics/apsis-leap/tests/time.rs new file mode 100644 index 0000000..570fc32 --- /dev/null +++ b/crates/physics/apsis-leap/tests/time.rs @@ -0,0 +1,117 @@ +//! Tests for UTC and TAI conversions on the Modified Julian Date scale. + +use apsis_leap::{ + MJD_AT_UNIX_EPOCH, calendar_mjd, tai_minus_utc_at_mjd, tai_mjd_from_utc_mjd, + utc_mjd_from_tai_mjd, +}; + +const SECONDS_PER_DAY: f64 = 86_400.0; + +#[test] +fn unix_epoch_constant_matches_the_calendar_conversion() { + assert_eq!(MJD_AT_UNIX_EPOCH, calendar_mjd(1970, 1, 1)); +} + +#[test] +fn offset_is_rejected_before_the_table_begins() { + let before = calendar_mjd(1971, 12, 31); + assert_eq!(tai_minus_utc_at_mjd(before), None); + assert_eq!(tai_mjd_from_utc_mjd(before), None); + assert_eq!(tai_minus_utc_at_mjd(f64::NAN), None); + assert_eq!(tai_minus_utc_at_mjd(f64::INFINITY), None); +} + +#[test] +fn first_entry_supplies_the_initial_offset() { + let start = calendar_mjd(1972, 1, 1); + assert_eq!(tai_minus_utc_at_mjd(start), Some(10.0)); + // Just before midnight one day earlier the table does not apply. + assert_eq!(tai_minus_utc_at_mjd(start - 1.0), None); +} + +#[test] +fn offset_changes_exactly_on_the_effective_date() { + // The 2017-01-01 leap second raised the offset from 36 to 37. + let effective = calendar_mjd(2017, 1, 1); + assert_eq!(tai_minus_utc_at_mjd(effective - 1.0), Some(36.0)); + assert_eq!(tai_minus_utc_at_mjd(effective), Some(37.0)); + assert_eq!(tai_minus_utc_at_mjd(effective + 0.5), Some(37.0)); +} + +#[test] +fn every_entry_is_reachable_at_its_effective_instant() { + for entry in apsis_leap::LEAP_SECONDS { + let at_entry = tai_minus_utc_at_mjd(entry.effective_mjd()); + assert_eq!(at_entry, Some(entry.tai_minus_utc), "entry {entry:?}"); + } +} + +#[test] +fn offset_is_constant_between_entries() { + // Mid 2016 sits between the 2015-07-01 and 2017-01-01 entries. + let instant = calendar_mjd(2016, 6, 1); + assert_eq!(tai_minus_utc_at_mjd(instant), Some(36.0)); + assert_eq!(tai_minus_utc_at_mjd(instant + 0.25), Some(36.0)); +} + +#[test] +fn future_instants_keep_the_last_offset() { + let far_future = calendar_mjd(2100, 1, 1); + assert_eq!(tai_minus_utc_at_mjd(far_future), Some(37.0)); +} + +#[test] +fn utc_to_tai_adds_the_offset() { + let utc = calendar_mjd(2017, 6, 1) + 0.5; + let tai = tai_mjd_from_utc_mjd(utc).expect("instant is inside the table"); + // An MJD near 58000 has a resolution of about 7e-12 days, so the sum + // can only be compared to well below a microsecond. + assert!((tai - utc - 37.0 / SECONDS_PER_DAY).abs() < 1e-11); +} + +#[test] +fn utc_to_tai_round_trips() { + for (year, month, day) in [ + (1972, 1, 1), + (1980, 1, 6), + (1999, 1, 1), + (2015, 7, 1), + (2017, 1, 1), + (2024, 6, 15), + ] { + let utc = calendar_mjd(year, month, day) + 0.375; + let tai = tai_mjd_from_utc_mjd(utc).expect("instant is inside the table"); + let back = utc_mjd_from_tai_mjd(tai).expect("instant is inside the table"); + assert!( + (back - utc).abs() < 1e-10, + "round trip failed at {year}-{month}-{day}: {utc} vs {back}" + ); + } +} + +#[test] +fn tai_to_utc_recovers_the_offset_on_each_side_of_a_step() { + let effective = calendar_mjd(2017, 1, 1); + // One hour before the step the offset is 36 s, one hour after it is 37 s. + let before_utc = effective - 1.0 / 24.0; + let after_utc = effective + 1.0 / 24.0; + let before_tai = tai_mjd_from_utc_mjd(before_utc).expect("inside the table"); + let after_tai = tai_mjd_from_utc_mjd(after_utc).expect("inside the table"); + assert!((utc_mjd_from_tai_mjd(before_tai).expect("inside") - before_utc).abs() < 1e-10); + assert!((utc_mjd_from_tai_mjd(after_tai).expect("inside") - after_utc).abs() < 1e-10); +} + +#[test] +fn tai_to_utc_rejects_non_finite_input() { + assert_eq!(utc_mjd_from_tai_mjd(f64::NAN), None); + assert_eq!(utc_mjd_from_tai_mjd(f64::NEG_INFINITY), None); +} + +#[test] +fn tai_to_utc_at_the_boundary_stays_inside_the_table() { + // The first TAI instant that the table covers is 1972-01-01 00:00:10 + // TAI, because the offset there is ten seconds. + let boundary = calendar_mjd(1972, 1, 1) + 10.0 / SECONDS_PER_DAY; + let utc = utc_mjd_from_tai_mjd(boundary).expect("boundary is covered"); + assert!((utc - calendar_mjd(1972, 1, 1)).abs() < 1e-10); +}