apsis-rs: add validation examples for geodesy and time

Cargo.lock

@@ -25,6 +25,12 @@ version = "0.1.0"
 [[package]]
 name = "apsis-coord"
 version = "0.1.0"
+dependencies = [
+ "apsis-core",
+ "apsis-frames",
+ "apsis-geodesy",
+ "apsis-math",
+]
 
 [[package]]
 name = "apsis-core"
@@ -49,10 +55,18 @@ version = "0.1.0"
 [[package]]
 name = "apsis-frames"
 version = "0.1.0"
+dependencies = [
+ "apsis-core",
+ "apsis-math",
+]
 
 [[package]]
 name = "apsis-geodesy"
 version = "0.1.0"
+dependencies = [
+ "apsis-core",
+ "apsis-math",
+]
 
 [[package]]
 name = "apsis-handover"
@@ -195,6 +209,10 @@ version = "0.1.0"
 [[package]]
 name = "apsis-time"
 version = "0.1.0"
+dependencies = [
+ "apsis-core",
+ "apsis-leap",
+]
 
 [[package]]
 name = "apsis-tle"

crates/apsis-rs/Cargo.toml

@@ -11,6 +11,17 @@ publish = false
 name = "apsis"
 path = "src/lib.rs"
 
+# Examples that produce validation input. Each needs the layers it reads
+# from, so it is built only when the matching feature is enabled.
+
+[[example]]
+name = "dump_geodesy"
+required-features = ["physics"]
+
+[[example]]
+name = "dump_time"
+required-features = ["physics"]
+
 [lints]
 workspace = true
 

crates/apsis-rs/examples/dump_geodesy.rs

@@ -0,0 +1,36 @@
+//! Writes a grid of geodetic positions and the Earth-fixed Cartesian
+//! coordinates this workspace computes for them.
+//!
+//! The output is CSV on standard output, with a header row. It is consumed
+//! by `validation/python/coord/check_geodetic.py`, which recomputes the same
+//! positions with a reference implementation and compares the two.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features physics --example dump_geodesy \
+//!     > validation/data/expected/rust_geodetic.csv
+//! ```
+
+use apsis_coord::{Geodetic, WGS84};
+
+fn main() {
+    println!("latitude_degrees,longitude_degrees,height_m,x_m,y_m,z_m");
+    let latitudes = [-89.0_f64, -60.0, -30.0, 0.0, 30.0, 45.0, 60.0, 89.0];
+    let longitudes = [-179.0_f64, -90.0, -45.0, 0.0, 45.0, 90.0, 179.0];
+    let heights = [-1_000.0_f64, 0.0, 100.0, 3_000.0, 35_786_000.0];
+
+    for latitude in latitudes {
+        for longitude in longitudes {
+            for height in heights {
+                let point = Geodetic::new(latitude.to_radians(), longitude.to_radians(), height)
+                    .to_ecef_on(WGS84)
+                    .vector();
+                println!(
+                    "{latitude},{longitude},{height},{:.9},{:.9},{:.9}",
+                    point.x, point.y, point.z
+                );
+            }
+        }
+    }
+}

crates/apsis-rs/examples/dump_time.rs

@@ -0,0 +1,80 @@
+//! Writes the leap second offset and the time scale conversions this
+//! workspace computes over a range of instants.
+//!
+//! The output is CSV on standard output, with a header row. It is consumed
+//! by `validation/python/time/check_time.py`, which recomputes each instant
+//! with a reference implementation.
+//!
+//! The TAI instant is written in two parts, an integer Modified Julian Date
+//! and a fraction of a day, because a single `f64` cannot hold a Modified
+//! Julian Date near 60 000 to better than about 40 microseconds. The two
+//! parts let the checker compare instants to well below a microsecond.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features physics --example dump_time \
+//!     > validation/data/expected/rust_time.csv
+//! ```
+
+use apsis_leap::{LEAP_SECONDS, tai_minus_utc_at_mjd};
+use apsis_time::{Calendar, Epoch, TT_MINUS_TAI, TimeScale};
+
+/// Splits a Modified Julian Date into its integer day and its fraction.
+fn split(mjd: f64) -> (i64, f64) {
+    let day = mjd.floor();
+    (day as i64, mjd - day)
+}
+
+fn emit(calendar: Calendar, epoch: Epoch) {
+    let offset = tai_minus_utc_at_mjd(
+        epoch
+            .mjd(TimeScale::Utc)
+            .expect("every emitted instant is in UTC"),
+    )
+    .expect("every emitted instant is covered by the table");
+    let (tai_day, tai_frac) = split(epoch.tai_mjd());
+    println!(
+        "{:04},{:02},{:02},{:02},{:02},{:.3},{:.1},{tai_day},{tai_frac:.17}",
+        calendar.year,
+        calendar.month,
+        calendar.day,
+        calendar.hour,
+        calendar.minute,
+        calendar.second,
+        offset
+    );
+}
+
+fn main() {
+    println!("year,month,day,hour,minute,second,tai_minus_utc,tai_day,tai_frac");
+
+    // Sample every leap second boundary at several times of day, so that
+    // the step behaviour and the interior of each interval are covered.
+    for entry in LEAP_SECONDS {
+        for (hour, minute, second) in [(0, 0, 0.0), (6, 0, 0.0), (12, 0, 0.0), (18, 0, 0.0)] {
+            let calendar = Calendar::new(entry.year, entry.month, entry.day, hour, minute, second);
+            if let Some(epoch) = Epoch::from_calendar(calendar, TimeScale::Utc) {
+                emit(calendar, epoch);
+            }
+        }
+    }
+
+    // The last second before a leap second and the first second after it.
+    for (year, month, day, hour, minute, second) in [
+        (2016, 12, 31, 23, 59, 59.0),
+        (2017, 1, 1, 0, 0, 0.0),
+        (2016, 12, 31, 12, 0, 0.0),
+        (1998, 12, 31, 23, 59, 59.0),
+        (1999, 1, 1, 0, 0, 0.0),
+    ] {
+        let calendar = Calendar::new(year, month, day, hour, minute, second);
+        if let Some(epoch) = Epoch::from_calendar(calendar, TimeScale::Utc) {
+            emit(calendar, epoch);
+        }
+    }
+
+    // The Terrestrial Time offset is a defining constant, written here so
+    // that the checker can verify it without recomputing the calendar.
+    println!("# tt_minus_tai_seconds={TT_MINUS_TAI:.3}");
+}