Added time-keeping and kernel downloading
This commit is contained in:
1 parent
c5e75e3b0a
commit
0bb6d48c68
7 files changed
+634
-11
No files matched your search
@@ -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.
|
||||
|
Can't render this file because it contains an unexpected character in line 9 and column 63.
|
@@ -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:
|
||||
|
||||
Executable
+54
@@ -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
|
||||
+51
-10
@@ -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<T: Scalar> Add for Contribution<T> {
|
||||
/// 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<T: Scalar>(&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.
|
||||
|
||||
+245
@@ -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<Kernel> {
|
||||
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<PathBuf> {
|
||||
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<Almanac, KernelError> {
|
||||
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<PathBuf> {
|
||||
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}"
|
||||
);
|
||||
}
|
||||
}
|
||||
+4
-1
@@ -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;
|
||||
}
|
||||
+247
@@ -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<Epoch> 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<f64> for UniformEpoch {
|
||||
type Output = Self;
|
||||
|
||||
fn add(self, seconds: f64) -> Self {
|
||||
Self(self.0 + Duration::from_seconds(seconds))
|
||||
}
|
||||
}
|
||||
|
||||
impl AddAssign<f64> for UniformEpoch {
|
||||
fn add_assign(&mut self, seconds: f64) {
|
||||
*self = *self + seconds;
|
||||
}
|
||||
}
|
||||
|
||||
/// Adds a duration, which carries no scale of its own.
|
||||
impl Add<Duration> 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
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user