physics: add leap second table and UTC TAI conversion

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,
+};

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
+}

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<f64> {
+    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<f64> {
+    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<f64> {
+    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)
+}

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);
+}

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);
+}