From 0bb6d48c6849527104ffba8a5ae6ef508766799e Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Wed, 30 Sep 2026 17:13:31 -0400 Subject: [PATCH] Added time-keeping and kernel downloading --- data/kernels.tsv | 15 +++ readme.md | 18 +++ scripts/fetch-kernels.sh | 54 +++++++++ src/dynamics/mod.rs | 61 ++++++++-- src/kernels.rs | 245 ++++++++++++++++++++++++++++++++++++++ src/lib.rs | 5 +- src/time.rs | 247 +++++++++++++++++++++++++++++++++++++++ 7 files changed, 634 insertions(+), 11 deletions(-) create mode 100644 data/kernels.tsv create mode 100755 scripts/fetch-kernels.sh create mode 100644 src/kernels.rs create mode 100644 src/time.rs diff --git a/data/kernels.tsv b/data/kernels.tsv new file mode 100644 index 0000000..c8b4257 --- /dev/null +++ b/data/kernels.tsv @@ -0,0 +1,15 @@ +# ANISE kernels, pinned by URL and SHA-256. Fetch with scripts/fetch-kernels.sh. +# +# Not committed: together they are about 38 MB, and a checksum pins a file more +# precisely than a copy of it does. The date is when the checksum was taken, so a +# later mismatch can be told apart from a transfer error. Tab-separated; read by +# both the fetch script and src/kernels.rs, so there is one list rather than two +# that can disagree. +# +# NAIF's earth_latest_high_prec.bpc is deliberately NOT used. "Latest" is not a +# thing a reproducible result can cite, so the dated snapshot covering the same +# span is pinned instead. +# +# name sha256 downloaded url purpose +de440s.bsp c1c7feeab882263fc493a9d5a5b2ddd71b54826cdf65d8d17a76126b260a49f2 2026-09-30 https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440s.bsp Planetary and lunar ephemerides, 1849 to 2150. Sun and Moon positions for third-body and solar radiation pressure. +earth_000101_261226_260929.bpc df5510d1995ec56af62c2682affc38781e3aa4fda37722d2fa0a10e531b90b1d 2026-09-30 https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_000101_261226_260929.bpc Earth orientation, 2000-01-01 to 2026-12-26, reconstructed through 2026-09-29. Needed for inertial to Earth-fixed transformations. diff --git a/readme.md b/readme.md index da79a6f..bd66750 100644 --- a/readme.md +++ b/readme.md @@ -76,6 +76,24 @@ are measured here rather than cited. | Full-force LEO, 1 day | not yet measured | — | — | | Ensemble of N states, scaling vs. cores | not yet measured | — | — | +## Kernels + +Ephemerides and Earth orientation come from NAIF kernels, which are too large to +commit. `data/kernels.tsv` pins each one by URL and SHA-256: + +```sh +scripts/fetch-kernels.sh # into data/kernels +scripts/fetch-kernels.sh /somewhere # or wherever you keep them +``` + +The script verifies checksums and leaves an already-correct file alone. Point +`ASTROKIT_KERNELS` at the directory if it is not `data/kernels`. Tests that need +kernels say so and skip when they are absent. + +`earth_latest_high_prec.bpc` is deliberately not used: a result that cites "the +latest Earth orientation kernel" cannot be reproduced, so the dated snapshot +covering the same span is pinned instead. + ## Development Enable the pre-commit hook once per clone: diff --git a/scripts/fetch-kernels.sh b/scripts/fetch-kernels.sh new file mode 100755 index 0000000..f924cf5 --- /dev/null +++ b/scripts/fetch-kernels.sh @@ -0,0 +1,54 @@ +#!/bin/sh +# Fetches the ANISE kernels listed in data/kernels.tsv and verifies their +# checksums. Idempotent: a file that is already present and correct is left alone. +# +# Usage: scripts/fetch-kernels.sh [directory] +# ASTROKIT_KERNELS=/somewhere scripts/fetch-kernels.sh + +set -e + +root=$(CDPATH= cd -- "$(dirname -- "$0")/.." && pwd) +manifest="$root/data/kernels.tsv" +target=${1:-${ASTROKIT_KERNELS:-$root/data/kernels}} + +if [ ! -f "$manifest" ]; then + echo "no manifest at $manifest" >&2 + exit 1 +fi + +mkdir -p "$target" +echo "kernels -> $target" + +status=0 +while IFS="$(printf '\t')" read -r name sha downloaded url purpose; do + case "$name" in '' | \#*) continue ;; esac + + path="$target/$name" + if [ -f "$path" ] && [ "$(sha256sum "$path" | cut -d' ' -f1)" = "$sha" ]; then + echo " ok $name" + continue + fi + + echo " fetching $name (pinned $downloaded)" + if ! curl -sfL --retry 3 -o "$path.partial" "$url"; then + echo " FAILED $name: could not download $url" >&2 + rm -f "$path.partial" + status=1 + continue + fi + + got=$(sha256sum "$path.partial" | cut -d' ' -f1) + if [ "$got" != "$sha" ]; then + echo " FAILED $name: checksum mismatch" >&2 + echo " expected $sha" >&2 + echo " got $got" >&2 + rm -f "$path.partial" + status=1 + continue + fi + + mv "$path.partial" "$path" + echo " verified $name" +done < "$manifest" + +exit $status diff --git a/src/dynamics/mod.rs b/src/dynamics/mod.rs index fd68eaa..7583b91 100644 --- a/src/dynamics/mod.rs +++ b/src/dynamics/mod.rs @@ -42,7 +42,7 @@ //! Assembly sums them as `acceleration + force / mass`, which mirrors what the //! rotational side will do with torque and the inertia tensor. -use hifitime::{Duration, Epoch}; +use hifitime::Epoch; use nalgebra::Vector3; use ordinary_diffeq::scalar::Scalar; use serde::{Deserialize, Serialize}; @@ -51,6 +51,7 @@ use std::ops::{Add, AddAssign}; use crate::body::Body; use crate::frames::Frame; use crate::state::{HasPosition, HasVelocity, Layout, State}; +use crate::time::UniformEpoch; pub mod models; @@ -129,9 +130,9 @@ impl Add for Contribution { /// What a [`DynamicsModel`] needs that is not in the state vector. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Context { - /// The epoch the propagation started from. The `t` a model receives is - /// seconds after this. - pub epoch: Epoch, + /// The epoch the propagation started from, in uniform time. The `t` a model + /// receives is seconds after this. + pub epoch: UniformEpoch, /// The frame the state is expressed in, and so the central body a model can /// assume without being told. pub frame: Frame, @@ -139,19 +140,36 @@ pub struct Context { impl Context { /// A context anchored at `epoch`, in `frame`. + /// + /// The epoch is normalized to the uniform scale a propagation runs in, so it + /// may be given in any scale. pub fn new(epoch: Epoch, frame: Frame) -> Self { - Self { epoch, frame } + Self { + epoch: UniformEpoch::new(epoch), + frame, + } } /// The epoch `t` seconds after the start of the propagation. /// - /// `t` is projected to its real part: an epoch is a real quantity, and there - /// is no such thing as a dual date. A propagation differentiated with respect - /// to its initial state is unaffected, since the time grid does not depend on - /// the state. It does mean sensitivities with respect to epoch itself are not + /// `t` is uniform seconds, which is what the integrator advances, so this goes + /// through [`crate::time::advance`] rather than adding a `Duration` directly + /// + /// `t` is projected to its real part: an epoch is a real quantity, and there is + /// no such thing as a dual date. A propagation differentiated with respect to + /// its initial state is unaffected, since the time grid does not depend on the + /// state. It does mean sensitivities with respect to epoch itself are not /// available this way. pub fn epoch_at(&self, t: T) -> Epoch { - self.epoch + Duration::from_seconds(t.re()) + (self.epoch + t.re()).epoch() + } + + /// The propagation time, in uniform seconds, at which `epoch` occurs. + /// + /// The inverse of [`Context::epoch_at`], for turning a measurement schedule + /// into integrator stop times. + pub fn seconds_at(&self, epoch: Epoch) -> f64 { + (UniformEpoch::new(epoch) - self.epoch).to_seconds() } } @@ -287,6 +305,29 @@ mod tests { context.epoch_at(3600.0), Epoch::from_gregorian_utc_hms(2026, 1, 1, 1, 0, 0) ); + assert_eq!(context.seconds_at(context.epoch_at(3600.0)), 3600.0); + } + + /// The integrator advances uniform seconds, so an arc spanning a leap second + /// has to land where the integrator thinks it did, not where a wall clock says. + #[test] + fn test_context_time_is_uniform_across_a_leap_second() { + let context = Context::new( + Epoch::from_gregorian_utc_hms(2016, 12, 31, 12, 0, 0), + Frame::EARTH_J2000, + ); + + let after_a_day = context.epoch_at(86400.0); + assert_eq!( + after_a_day.to_tai_seconds() - context.epoch.epoch().to_tai_seconds(), + 86400.0 + ); + // Which is one second before the next day's noon UTC, the leap second + // having been inserted in between. + assert_eq!( + context.seconds_at(Epoch::from_gregorian_utc_hms(2017, 1, 1, 12, 0, 0)), + 86401.0 + ); } /// A model bounded on a capability works for any layout that has it. diff --git a/src/kernels.rs b/src/kernels.rs new file mode 100644 index 0000000..a1123c2 --- /dev/null +++ b/src/kernels.rs @@ -0,0 +1,245 @@ +//! The ANISE kernels this crate reads, and loading them. +//! +//! Kernels are large binary files that are not committed. `data/kernels.tsv` pins +//! each one by URL and SHA-256, `scripts/fetch-kernels.sh` fetches and verifies +//! them, and this module reads the same manifest so there is one list rather than +//! two that can drift apart. +//! +//! NAIF's `earth_latest_high_prec.bpc` is deliberately not used. A result that +//! cites "the latest Earth orientation kernel" cannot be reproduced, so the dated +//! snapshot covering the same span is pinned instead. +//! +//! A kernel's location is the caller's to decide. [`directory`] reads +//! `ASTROKIT_KERNELS` for convenience, but nothing here guesses. + +use std::path::{Path, PathBuf}; + +use anise::almanac::Almanac; +use thiserror::Error; + +/// The environment variable [`directory`] reads. +pub const DIRECTORY_VARIABLE: &str = "ASTROKIT_KERNELS"; + +/// What went wrong loading kernels. +#[derive(Debug, Error)] +pub enum KernelError { + /// A manifest kernel is not in the directory. + #[error( + "kernel {name} is not in {}: run scripts/fetch-kernels.sh", + directory.display() + )] + Missing { + /// The kernel's file name. + name: String, + /// Where it was looked for. + directory: PathBuf, + }, + + /// ANISE would not read a file that is present. + #[error("ANISE could not load {}: {message}", path.display())] + Unreadable { + /// The file that failed. + path: PathBuf, + /// What ANISE said. + message: String, + }, +} + +/// One kernel: what it is, where it came from, and what it hashes to. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct Kernel { + /// File name, which is also the name in the manifest. + pub name: &'static str, + /// Expected SHA-256, lowercase hex. + pub sha256: &'static str, + /// The date the checksum was taken, `YYYY-MM-DD`. + /// + /// Lets a later mismatch be told apart from a transfer error: if the file at + /// the URL has been revised since this date, that is a different problem from a + /// corrupted download. + pub downloaded: &'static str, + /// Where it was fetched from. + pub url: &'static str, + /// What it is for. + pub purpose: &'static str, +} + +/// The manifest, embedded at compile time so it cannot go missing at runtime. +const MANIFEST: &str = include_str!("../data/kernels.tsv"); + +/// The kernels this crate expects. +/// +/// Suitable for a provenance record: every entry names a file and the hash that +/// identifies its exact contents. +pub fn manifest() -> Vec { + MANIFEST + .lines() + .map(str::trim_end) + .filter(|line| !line.is_empty() && !line.starts_with('#')) + .filter_map(|line| { + let mut fields = line.split('\t'); + Some(Kernel { + name: fields.next()?, + sha256: fields.next()?, + downloaded: fields.next()?, + url: fields.next()?, + purpose: fields.next().unwrap_or(""), + }) + }) + .collect() +} + +/// The kernel directory named by `ASTROKIT_KERNELS`, if it is set. +pub fn directory() -> Option { + std::env::var_os(DIRECTORY_VARIABLE).map(PathBuf::from) +} + +/// Loads every kernel in the manifest from `directory`. +/// +/// All or nothing: a missing kernel is an error rather than a quietly smaller +/// almanac, because a propagation that silently lost its ephemerides would look +/// like a physics problem rather than a setup one. +pub fn load(directory: &Path) -> Result { + let mut almanac = Almanac::default(); + + for kernel in manifest() { + let path = directory.join(kernel.name); + if !path.is_file() { + return Err(KernelError::Missing { + name: kernel.name.to_string(), + directory: directory.to_path_buf(), + }); + } + + let as_str = path.to_string_lossy().into_owned(); + almanac = almanac + .load(&as_str) + .map_err(|error| KernelError::Unreadable { + path: path.clone(), + message: error.to_string(), + })?; + } + + Ok(almanac) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Where the tests look for kernels: the variable if set, else the repository's + /// own directory. Returns `None` when they have not been fetched. + fn available() -> Option { + let candidate = directory() + .unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels")); + manifest() + .iter() + .all(|kernel| candidate.join(kernel.name).is_file()) + .then_some(candidate) + } + + #[test] + fn test_manifest_parses() { + let kernels = manifest(); + assert_eq!(kernels.len(), 2, "{kernels:#?}"); + + for kernel in &kernels { + assert!(!kernel.name.is_empty()); + assert_eq!(kernel.sha256.len(), 64, "{} sha256", kernel.name); + assert!(kernel.sha256.chars().all(|c| c.is_ascii_hexdigit())); + assert!(kernel.url.starts_with("https://"), "{}", kernel.url); + assert_eq!(kernel.downloaded.len(), 10, "{} date", kernel.name); + assert!( + kernel.downloaded.split('-').count() == 3, + "{} date is not YYYY-MM-DD: {}", + kernel.name, + kernel.downloaded + ); + assert!(!kernel.purpose.is_empty(), "{} purpose", kernel.name); + } + } + + /// "Latest" kernels are not reproducible, so the manifest must not name one. + #[test] + fn test_manifest_pins_no_moving_target() { + for kernel in manifest() { + assert!( + !kernel.name.contains("latest") && !kernel.url.contains("latest"), + "{} is a moving target", + kernel.name + ); + } + } + + #[test] + fn test_missing_kernel_says_what_to_do() { + // `Almanac` is not `Debug`, so match rather than unwrap_err. + let message = match load(Path::new("/nonexistent")) { + Err(error) => error.to_string(), + Ok(_) => panic!("loaded kernels from a directory that does not exist"), + }; + assert!(message.contains("fetch-kernels.sh"), "{message}"); + } + + #[test] + fn test_loads_the_manifest() { + let Some(directory) = available() else { + eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh"); + return; + }; + + let almanac = load(&directory).unwrap(); + assert_eq!(almanac.num_loaded_spk(), 1); + assert_eq!(almanac.num_loaded_bpc(), 1); + } + + /// The ephemerides have to actually answer a query, not merely parse. + #[test] + fn test_ephemerides_place_the_sun_about_one_astronomical_unit_away() { + let Some(directory) = available() else { + eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh"); + return; + }; + + use anise::constants::frames::{EARTH_J2000, SUN_J2000}; + use hifitime::Epoch; + + let almanac = load(&directory).unwrap(); + let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0); + let state = almanac + .transform(SUN_J2000, EARTH_J2000, epoch, None) + .unwrap(); + + // ANISE works in kilometers; one astronomical unit to within Earth's + // orbital eccentricity, which peaks near perihelion in early January. + let distance_m = state.radius_km.norm() * 1e3; + assert!( + (1.47e11..1.53e11).contains(&distance_m), + "Sun at {distance_m:e} m, which is not about an astronomical unit" + ); + } + + /// And the Earth orientation kernel has to answer a rotation query, which is + /// what the frame transformations will need. + #[test] + fn test_earth_orientation_rotates_within_its_covered_span() { + let Some(directory) = available() else { + eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh"); + return; + }; + + use anise::constants::frames::{EARTH_ITRF93, EARTH_J2000}; + use hifitime::Epoch; + + let almanac = load(&directory).unwrap(); + let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0); + let rotation = almanac.rotate(EARTH_J2000, EARTH_ITRF93, epoch).unwrap(); + + // A rotation, so orthonormal: the determinant of the 3x3 part is one. + let determinant = rotation.rot_mat.determinant(); + assert!( + (determinant - 1.0).abs() < 1e-9, + "rotation determinant {determinant}" + ); + } +} diff --git a/src/lib.rs b/src/lib.rs index 9eaf64e..90f6d51 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -15,8 +15,10 @@ pub mod dynamics; pub mod equinoctial; pub mod frames; pub mod kepler; +pub mod kernels; pub mod propagation; pub mod state; +pub mod time; /// The types you want in scope to set up and run a propagation. pub mod prelude { @@ -34,7 +36,8 @@ pub mod prelude { Full, HasAngularRate, HasAttitude, HasEnergy, HasMass, HasPosition, HasVelocity, Layout, Orbit, OrbitMass, State, }; - pub use hifitime::{Duration, Epoch, Unit}; + pub use crate::time; + pub use hifitime::{Duration, Epoch, TimeScale, Unit}; pub use ordinary_diffeq::prelude::{Method, MethodName, Tolerances}; pub use ordinary_diffeq::scalar::Scalar; } diff --git a/src/time.rs b/src/time.rs new file mode 100644 index 0000000..a20b780 --- /dev/null +++ b/src/time.rs @@ -0,0 +1,247 @@ +//! Epochs, and the uniform time a propagation runs in. +//! +//! # Propagation time is uniform, and that is not automatic +//! +//! An integrator advances a state by seconds of *uniform* time. Calendar time is +//! not uniform: UTC inserts a leap second occasionally, so two UTC timestamps a +//! nominal day apart can be 86401 seconds of real elapsed time apart. +//! +//! `hifitime` does exactly what it is asked here, and what it is asked depends on +//! the epoch's own time scale: +//! +//! ``` +//! use hifitime::{Duration, Epoch, TimeScale}; +//! +//! // The 2016 leap second falls inside this day. +//! let utc = Epoch::from_gregorian_utc_hms(2016, 12, 31, 12, 0, 0); +//! let tai = utc.to_time_scale(TimeScale::TAI); +//! +//! let advanced_in_utc = utc + Duration::from_seconds(86400.0); +//! let advanced_in_tai = tai + Duration::from_seconds(86400.0); +//! +//! // Adding to a UTC epoch adds UTC seconds, which is 86401 seconds of real time. +//! assert_eq!(advanced_in_utc.to_tai_seconds() - utc.to_tai_seconds(), 86401.0); +//! // Adding to a TAI epoch adds uniform seconds, which is what an integrator did. +//! assert_eq!(advanced_in_tai.to_tai_seconds() - tai.to_tai_seconds(), 86400.0); +//! ``` +//! +//! A one-second error puts a low Earth orbit about 7.7 km out of place, so this is +//! not a rounding concern. Everything in this crate that converts between a +//! propagation's `t` and an epoch goes through [`advance`] and +//! [`elapsed_seconds`], which convert to TAI first. Pass epochs in whatever scale +//! is convenient; they are normalized on the way in. + +use core::fmt; +use core::ops::{Add, AddAssign, Sub}; + +use hifitime::{Duration, Epoch, TimeScale}; + +/// The time scale a propagation's `t` is measured in. +/// +/// TAI: uniform, no leap seconds, no relativistic terms. An integrator's step size +/// means the same thing at every epoch. +pub const PROPAGATION_SCALE: TimeScale = TimeScale::TAI; + +/// An epoch in the uniform scale a propagation runs in. +/// +/// Adding seconds to one of these adds seconds. Adding seconds to a bare [`Epoch`] +/// adds seconds of whatever scale that epoch carries, which across a leap second is +/// a different thing. Wrapping the epoch makes the right arithmetic the only +/// arithmetic available, rather than something to remember. +#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub struct UniformEpoch(Epoch); + +impl UniformEpoch { + /// `epoch` converted to the uniform scale. Any scale in, uniform out. + pub fn new(epoch: Epoch) -> Self { + Self(uniform(epoch)) + } + + /// The underlying epoch, for display or for handing to ANISE. + pub fn epoch(self) -> Epoch { + self.0 + } +} + +impl From for UniformEpoch { + fn from(epoch: Epoch) -> Self { + Self::new(epoch) + } +} + +impl fmt::Display for UniformEpoch { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "{}", self.0) + } +} + +/// Adds uniform seconds. +impl Add for UniformEpoch { + type Output = Self; + + fn add(self, seconds: f64) -> Self { + Self(self.0 + Duration::from_seconds(seconds)) + } +} + +impl AddAssign for UniformEpoch { + fn add_assign(&mut self, seconds: f64) { + *self = *self + seconds; + } +} + +/// Adds a duration, which carries no scale of its own. +impl Add for UniformEpoch { + type Output = Self; + + fn add(self, duration: Duration) -> Self { + Self(self.0 + duration) + } +} + +/// The interval between two uniform epochs, later minus earlier, as `hifitime` +/// does it for [`Epoch`]. +/// +/// `.to_seconds()` on the result is the interval an integrator has to cover. +impl Sub for UniformEpoch { + type Output = Duration; + + fn sub(self, earlier: Self) -> Duration { + self.0 - earlier.0 + } +} + +/// `epoch` re-expressed in the uniform scale a propagation runs in. +/// +/// Cheap, and idempotent on an epoch that is already uniform. +pub fn uniform(epoch: Epoch) -> Epoch { + epoch.to_time_scale(PROPAGATION_SCALE) +} + +/// The epoch `seconds` of uniform time after `epoch`. +/// +/// Not `epoch + Duration::from_seconds(seconds)`: that adds seconds of the epoch's +/// own scale, which is wrong for a UTC epoch spanning a leap second. +pub fn advance(epoch: Epoch, seconds: f64) -> Epoch { + uniform(epoch) + Duration::from_seconds(seconds) +} + +/// The uniform seconds elapsed from `from` to `to`, negative if `to` is earlier. +/// +/// Counts leap seconds as the seconds they are, so this is the interval an +/// integrator would have to cover, not the difference of two wall clocks. +pub fn elapsed_seconds(from: Epoch, to: Epoch) -> f64 { + (uniform(to) - uniform(from)).to_seconds() +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + /// Noon to noon across the 2016 leap second: 86401 uniform seconds, not 86400. + fn before_leap_second() -> Epoch { + Epoch::from_gregorian_utc_hms(2016, 12, 31, 12, 0, 0) + } + + fn after_leap_second() -> Epoch { + Epoch::from_gregorian_utc_hms(2017, 1, 1, 12, 0, 0) + } + + /// The newtype exists so that plain addition is the correct addition. + #[test] + fn test_uniform_epoch_arithmetic_is_uniform() { + let start = UniformEpoch::new(before_leap_second()); + let advanced = start + 86400.0; + + assert_relative_eq!( + advanced.epoch().to_tai_seconds() - start.epoch().to_tai_seconds(), + 86400.0, + epsilon = 1e-9 + ); + assert_relative_eq!((advanced - start).to_seconds(), 86400.0, epsilon = 1e-9); + // Later minus earlier is positive, matching hifitime's own convention. + assert_relative_eq!((start - advanced).to_seconds(), -86400.0, epsilon = 1e-9); + + // Which agrees with the free function, and is idempotent under conversion. + assert_eq!(advanced.epoch(), advance(before_leap_second(), 86400.0)); + assert_eq!(UniformEpoch::new(advanced.epoch()), advanced); + } + + #[test] + fn test_uniform_epoch_accepts_any_scale() { + let utc = before_leap_second(); + for epoch in [utc, uniform(utc), utc.to_time_scale(TimeScale::TDB)] { + assert_eq!(UniformEpoch::new(epoch), UniformEpoch::from(utc)); + } + } + + #[test] + fn test_advance_adds_uniform_seconds_across_a_leap_second() { + let start = before_leap_second(); + let advanced = advance(start, 86400.0); + + // Exactly the seconds asked for, in real elapsed time. + assert_relative_eq!( + advanced.to_tai_seconds() - start.to_tai_seconds(), + 86400.0, + epsilon = 1e-9 + ); + + // Which is one second short of the next day's noon, because a leap second + // was inserted in between. + assert_relative_eq!( + elapsed_seconds(advanced, after_leap_second()), + 1.0, + epsilon = 1e-9 + ); + } + + /// The naive form, kept as a test so the difference is recorded rather than + /// remembered. + #[test] + fn test_naive_addition_is_wrong_across_a_leap_second() { + let start = before_leap_second(); + let naive = start + Duration::from_seconds(86400.0); + + assert_relative_eq!( + naive.to_tai_seconds() - start.to_tai_seconds(), + 86401.0, + epsilon = 1e-9 + ); + assert_relative_eq!( + advance(start, 86400.0).to_tai_seconds(), + naive.to_tai_seconds() - 1.0, + epsilon = 1e-9 + ); + } + + #[test] + fn test_elapsed_seconds_counts_the_leap_second() { + assert_relative_eq!( + elapsed_seconds(before_leap_second(), after_leap_second()), + 86401.0, + epsilon = 1e-9 + ); + // And is antisymmetric. + assert_relative_eq!( + elapsed_seconds(after_leap_second(), before_leap_second()), + -86401.0, + epsilon = 1e-9 + ); + } + + #[test] + fn test_elapsed_seconds_is_scale_agnostic() { + let utc = before_leap_second(); + let tai = uniform(utc); + let tdb = utc.to_time_scale(TimeScale::TDB); + for epoch in [utc, tai, tdb] { + assert_relative_eq!( + elapsed_seconds(epoch, after_leap_second()), + 86401.0, + epsilon = 1e-9 + ); + } + } +}