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

physics: add leap second table and UTC TAI conversion

This commit is contained in:
ApsisRS contributors 2026-09-26 13:15:52 +08:00
parent 0b22758a8d
commit 8b5213747a
5 changed files with 591 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
crates/physics/apsis-leap/src/lib.rs
View File

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

241
crates/physics/apsis-leap/src/table.rs Normal file
View File

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

108
crates/physics/apsis-leap/src/time.rs Normal file
View File

@ -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)
}

80
crates/physics/apsis-leap/tests/table.rs Normal file
View File

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

117
crates/physics/apsis-leap/tests/time.rs Normal file
View File

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