physics: add frame markers and Earth rotation angle

crates/physics/apsis-frames/Cargo.toml

@@ -11,3 +11,5 @@ publish = false
 workspace = true
 
 [dependencies]
+apsis-core.workspace = true
+apsis-math.workspace = true

crates/physics/apsis-frames/src/earth_rotation.rs

@@ -0,0 +1,143 @@
+//! The Earth Rotation Angle and the rotation it defines between the
+//! inertial and Earth-fixed frames.
+
+use crate::frame::{Ecef, Eci};
+use apsis_math::Mat3;
+use core::f64::consts::PI;
+
+/// The Earth Rotation Angle at the J2000.0 epoch, in turns.
+///
+/// The epoch value is zero by the definition of the angle: the International
+/// Astronomical Union defines the Earth Rotation Angle to be zero at
+/// 2000-01-01T12:00:00 UT1.
+pub const ERA_J2000_TURNS: f64 = 0.0;
+
+/// The rate of change of the Earth Rotation Angle, in turns per day of UT1.
+///
+/// This is `1.002 737 811 911 354 6`, which is `360.985 612 288 088`
+/// degrees per day. It exceeds one turn per day because the Earth turns once
+/// in a sidereal day, which is shorter than a solar day.
+///
+/// Source: IERS Conventions (2010), equation 8.32.
+pub const ERA_RATE_TURNS_PER_DAY: f64 = 1.002_737_811_911_354_6;
+
+/// The number of days from the J2000.0 epoch to the given Julian Date on the
+/// UT1 scale.
+///
+/// The J2000.0 epoch is Julian Date 2451545.0, which is noon on
+/// 2000-01-01.
+pub fn ut1_days_from_j2000(ut1_julian_date: f64) -> f64 {
+    ut1_julian_date - 2_451_545.0
+}
+
+/// The Earth Rotation Angle at a given instant.
+///
+/// The angle is the rotation from the celestial intermediate origin to the
+/// terrestrial intermediate origin, measured about the celestial
+/// intermediate pole. It replaces Greenwich Mean Sidereal Time in modern
+/// practice and, unlike sidereal time, is a linear function of UT1.
+///
+/// The value is held in turns. A turn is one full revolution, so the angle
+/// is periodic with a period of one. The value returned is not reduced to
+/// any interval, so differences between two angles are meaningful even when
+/// the interval is longer than a day.
+#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
+pub struct EarthRotationAngle {
+    turns: f64,
+}
+
+impl EarthRotationAngle {
+    /// Creates an angle from a value in turns.
+    pub const fn from_turns(turns: f64) -> Self {
+        Self { turns }
+    }
+
+    /// Creates an angle from a value in radians.
+    pub fn from_radians(radians: f64) -> Self {
+        Self {
+            turns: radians / (2.0 * PI),
+        }
+    }
+
+    /// Returns the angle for a Julian Date on the UT1 scale.
+    pub fn from_ut1_julian_date(ut1_julian_date: f64) -> Self {
+        Self::from_ut1_days(ut1_days_from_j2000(ut1_julian_date))
+    }
+
+    /// Returns the angle for a number of days of UT1 since the J2000.0
+    /// epoch.
+    ///
+    /// The expression is linear: `theta = 0.779 057 273 2640 + 1.002 737
+    /// 909 350 795 * d`, where `d` is the interval in days. The constant
+    /// term of `0.779 057` turns is the angle at the epoch itself, which is
+    /// about 280.46 degrees.
+    ///
+    /// Source: IERS Conventions (2010), equation 8.32.
+    pub fn from_ut1_days(days_since_j2000: f64) -> Self {
+        const EPOCH_TURNS: f64 = 0.779_057_273_264_0;
+        Self {
+            turns: EPOCH_TURNS + ERA_RATE_TURNS_PER_DAY * days_since_j2000,
+        }
+    }
+
+    /// Returns the angle in turns.
+    ///
+    /// The value is not reduced to `[0, 1)`.
+    pub const fn turns(self) -> f64 {
+        self.turns
+    }
+
+    /// Returns the angle in turns, reduced to `[0, 1)`.
+    pub fn turns_positive(self) -> f64 {
+        self.turns.rem_euclid(1.0)
+    }
+
+    /// Returns the angle in radians.
+    ///
+    /// The value is not reduced to `[0, 2 pi)`.
+    pub fn radians(self) -> f64 {
+        self.turns * 2.0 * PI
+    }
+
+    /// Returns the angle in radians, reduced to `[0, 2 pi)`.
+    pub fn radians_positive(self) -> f64 {
+        self.turns_positive() * 2.0 * PI
+    }
+}
+
+/// Returns the rotation that maps a vector in the inertial frame to the
+/// Earth-fixed frame.
+///
+/// The rotation is a clockwise turn about the third axis by the Earth
+/// Rotation Angle. Applying it to an inertial vector yields the Earth-fixed
+/// vector, so a vector that is fixed in inertial space acquires a negative
+/// rotation as the Earth turns beneath it.
+pub fn rotation_eci_to_ecef(angle: EarthRotationAngle) -> Mat3 {
+    let (sin, cos) = angle.radians_positive().sin_cos();
+    Mat3::from_rows([cos, sin, 0.0], [-sin, cos, 0.0], [0.0, 0.0, 1.0])
+}
+
+/// Returns the rotation that maps a vector in the Earth-fixed frame to the
+/// inertial frame.
+///
+/// This is the transpose of [`rotation_eci_to_ecef`], which for a rotation
+/// matrix is also its inverse.
+pub fn rotation_ecef_to_eci(angle: EarthRotationAngle) -> Mat3 {
+    rotation_eci_to_ecef(angle).transpose()
+}
+
+impl Eci {
+    /// Converts this inertial vector to the Earth-fixed frame at the given
+    /// Earth Rotation Angle.
+    pub fn to_ecef(self, angle: EarthRotationAngle) -> Ecef {
+        Ecef::new(rotation_eci_to_ecef(angle) * self.0)
+    }
+}
+
+impl Ecef {
+    /// Converts this Earth-fixed vector to the inertial frame at the given
+    /// Earth Rotation Angle.
+    pub fn to_eci(self, angle: EarthRotationAngle) -> Eci {
+        Eci::new(rotation_ecef_to_eci(angle) * self.0)
+    }
+}

crates/physics/apsis-frames/src/frame.rs

@@ -0,0 +1,41 @@
+//! Frame marker types.
+
+use apsis_math::Vec3;
+
+/// A position or vector expressed in the Earth-centered inertial frame.
+///
+/// The frame is aligned with the mean equator and equinox of J2000.0. It is
+/// treated as inertial here, so it does not rotate with the Earth.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Eci(pub Vec3);
+
+/// A position or vector expressed in the Earth-centered Earth-fixed frame.
+///
+/// The frame rotates with the Earth, so a point fixed to the ground has
+/// constant coordinates in it.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Ecef(pub Vec3);
+
+impl Eci {
+    /// Creates an inertial frame vector from its components.
+    pub const fn new(vector: Vec3) -> Self {
+        Self(vector)
+    }
+
+    /// Returns the underlying vector.
+    pub const fn vector(self) -> Vec3 {
+        self.0
+    }
+}
+
+impl Ecef {
+    /// Creates an Earth-fixed frame vector from its components.
+    pub const fn new(vector: Vec3) -> Self {
+        Self(vector)
+    }
+
+    /// Returns the underlying vector.
+    pub const fn vector(self) -> Vec3 {
+        self.0
+    }
+}

crates/physics/apsis-frames/src/lib.rs

@@ -1 +1,52 @@
-//! Reference frame definitions and transformation primitives.
+//! Reference frame definitions and the transformations between them.
+//!
+//! A reference frame fixes the orientation and origin that a vector is
+//! expressed relative to. This crate defines the two frames used for
+//! near-Earth work and the rotation between them.
+//!
+//! The inertial frame is the Earth-centered inertial (ECI) frame aligned
+//! with the mean equator and equinox of J2000.0. Its axes are fixed in
+//! inertial space, so an orbit propagated in it is not affected by the
+//! rotation of the Earth.
+//!
+//! The rotating frame is the Earth-centered Earth-fixed (ECEF) frame. Its
+//! third axis is the Earth's rotation axis and its first axis follows the
+//! prime meridian, so a point fixed to the ground has constant
+//! coordinates in it.
+//!
+//! The two are related by a rotation of [`EarthRotationAngle`] about the
+//! common third axis. This crate models the Earth's rotation with the
+//! Earth Rotation Angle, which is a function of UT1 only. It does not
+//! apply precession or nutation, because an angle accurate to a few
+//! hundredths of a degree is sufficient for pointing and coverage work and
+//! the full model needs the orientation of the celestial pole, which is a
+//! separate computation.
+//!
+//! Source: the Earth Rotation Angle expression is from the International
+//! Astronomical Union, as given in the *IERS Conventions (2010)*, equation
+//! 8.32.
+//!
+//! ```
+//! use apsis_frames::{EarthRotationAngle, rotation_eci_to_ecef};
+//!
+//! // The angle advances by roughly 360.9856 degrees per day.
+//! let first = EarthRotationAngle::from_ut1_days(0.0);
+//! let second = EarthRotationAngle::from_ut1_days(1.0);
+//! let advance = (second.radians() - first.radians()).to_degrees();
+//! assert!((advance - 360.985_612_288_088).abs() < 1e-6);
+//!
+//! // The rotation is ready to apply to a vector.
+//! let _ = rotation_eci_to_ecef(first);
+//! ```
+
+#![forbid(unsafe_code)]
+#![deny(missing_docs)]
+
+mod earth_rotation;
+mod frame;
+
+pub use earth_rotation::{
+    ERA_J2000_TURNS, ERA_RATE_TURNS_PER_DAY, EarthRotationAngle, rotation_ecef_to_eci,
+    rotation_eci_to_ecef, ut1_days_from_j2000,
+};
+pub use frame::{Ecef, Eci};

crates/physics/apsis-frames/tests/earth_rotation.rs

@@ -0,0 +1,181 @@
+//! Tests for the Earth Rotation Angle and the frame rotations.
+
+use apsis_frames::{
+    ERA_RATE_TURNS_PER_DAY, EarthRotationAngle, Ecef, Eci, rotation_ecef_to_eci,
+    rotation_eci_to_ecef, ut1_days_from_j2000,
+};
+use apsis_math::Vec3;
+use core::f64::consts::PI;
+
+const EPS: f64 = 1e-12;
+
+#[test]
+fn rate_matches_the_published_value() {
+    // 1.002 737 811 911 354 6 turns per day is 360.985 612 288 088 degrees
+    // per day.
+    assert!((ERA_RATE_TURNS_PER_DAY - 1.002_737_811_911_354_6).abs() < 1e-15);
+    let degrees = ERA_RATE_TURNS_PER_DAY * 360.0;
+    assert!(
+        (degrees - 360.985_612_288_088).abs() < 1e-9,
+        "got {degrees}"
+    );
+}
+
+#[test]
+fn angle_advances_by_the_rate_each_day() {
+    let first = EarthRotationAngle::from_ut1_days(0.0);
+    let second = EarthRotationAngle::from_ut1_days(1.0);
+    assert!((second.turns() - first.turns() - ERA_RATE_TURNS_PER_DAY).abs() < 1e-15);
+    // The angle advances a little under one degree per day more than a
+    // full turn.
+    let advance_degrees = (second.turns() - first.turns()) * 360.0 - 360.0;
+    assert!((advance_degrees - 0.985_612_288_088).abs() < 1e-9);
+}
+
+#[test]
+fn epoch_angle_matches_the_published_constant() {
+    // At J2000.0 the angle is 0.779 057 273 2640 turns, about 280.46
+    // degrees.
+    let epoch = EarthRotationAngle::from_ut1_days(0.0);
+    assert!((epoch.turns() - 0.779_057_273_264).abs() < 1e-15);
+    let degrees = epoch.turns_positive() * 360.0;
+    assert!((degrees - 280.460_618_375).abs() < 1e-9, "got {degrees}");
+}
+
+#[test]
+fn julian_date_and_day_count_agree() {
+    // J2000.0 is Julian Date 2451545.0.
+    assert_eq!(ut1_days_from_j2000(2_451_545.0), 0.0);
+    assert_eq!(ut1_days_from_j2000(2_451_546.0), 1.0);
+    let via_date = EarthRotationAngle::from_ut1_julian_date(2_451_546.0);
+    let via_days = EarthRotationAngle::from_ut1_days(1.0);
+    assert!((via_date.turns() - via_days.turns()).abs() < 1e-15);
+}
+
+#[test]
+fn radians_and_turns_agree() {
+    let angle = EarthRotationAngle::from_ut1_days(3.25);
+    assert!((angle.radians() - angle.turns() * 2.0 * PI).abs() < 1e-12);
+    let from_radians = EarthRotationAngle::from_radians(1.25);
+    assert!((from_radians.turns() - 1.25 / (2.0 * PI)).abs() < 1e-15);
+}
+
+#[test]
+fn positive_variants_reduce_to_the_principal_interval() {
+    for days in [-5.0, -1.0, 0.0, 1.0, 7.5, 100.0] {
+        let angle = EarthRotationAngle::from_ut1_days(days);
+        let turns = angle.turns_positive();
+        assert!((0.0..1.0).contains(&turns), "turns out of range: {turns}");
+        // Reduction removes whole turns, so the angle is unchanged modulo
+        // one turn.
+        let difference = angle.turns() - turns;
+        assert!((difference - difference.round()).abs() < 1e-9);
+
+        let radians = angle.radians_positive();
+        assert!((0.0..2.0 * PI).contains(&radians));
+    }
+}
+
+#[test]
+fn rotation_is_orthonormal_with_unit_determinant() {
+    let rotation = rotation_eci_to_ecef(EarthRotationAngle::from_ut1_days(12.5));
+    assert!(rotation.is_rotation(1e-12));
+    assert!((rotation.determinant() - 1.0).abs() < EPS);
+}
+
+#[test]
+fn a_positive_angle_drifts_inertial_directions_westward() {
+    // The Earth turns eastward, so a direction fixed in inertial space
+    // drifts westward when read in the Earth-fixed frame. In the
+    // Earth-fixed frame that is toward negative second axis.
+    let angle = EarthRotationAngle::from_radians(0.25);
+    let rotation = rotation_eci_to_ecef(angle);
+    let rotated = rotation * Vec3::X;
+    assert!(rotated.y < 0.0, "expected negative y, got {}", rotated.y);
+    assert!((rotated.norm() - 1.0).abs() < EPS);
+}
+
+#[test]
+fn opposite_rotations_cancel() {
+    let angle = EarthRotationAngle::from_ut1_days(4.75);
+    let forward = rotation_eci_to_ecef(angle);
+    let backward = rotation_ecef_to_eci(angle);
+    let identity = forward * backward;
+    assert!(identity.is_rotation(1e-12));
+    // Applying one after the other returns the original vector.
+    let vector = Vec3::new(1.5, -2.5, 3.5);
+    assert!((backward * (forward * vector) - vector).norm() < EPS);
+}
+
+#[test]
+fn earth_rotation_turns_the_third_axis_into_itself() {
+    let angle = EarthRotationAngle::from_ut1_days(2.0);
+    let rotation = rotation_eci_to_ecef(angle);
+    assert!((rotation * Vec3::Z - Vec3::Z).norm() < EPS);
+}
+
+#[test]
+fn rotation_carries_the_second_axis_toward_the_first() {
+    // The companion to the westward drift test: the second axis of the
+    // inertial frame maps toward the first axis of the Earth-fixed frame.
+    let angle = EarthRotationAngle::from_radians(0.25);
+    let rotation = rotation_eci_to_ecef(angle);
+    let rotated = rotation * Vec3::Y;
+    assert!(rotated.x > 0.0, "expected positive x, got {}", rotated.x);
+    assert!((rotated.norm() - 1.0).abs() < EPS);
+}
+
+#[test]
+fn frame_conversions_round_trip() {
+    let angle = EarthRotationAngle::from_ut1_days(9.25);
+    let original = Eci::new(Vec3::new(7_000.0, -1_500.0, 3_200.0));
+    let ecef = original.to_ecef(angle);
+    let back = ecef.to_eci(angle);
+    assert!((back.vector() - original.vector()).norm() < 1e-9);
+}
+
+#[test]
+fn frame_conversion_preserves_length() {
+    let angle = EarthRotationAngle::from_ut1_days(-3.5);
+    let vector = Vec3::new(6_878_137.0, 0.0, 0.0);
+    let converted = Eci::new(vector).to_ecef(angle);
+    assert!((converted.vector().norm() - vector.norm()).abs() < 1e-6);
+}
+
+#[test]
+fn rotating_background_returns_a_ground_point_to_the_same_longitude() {
+    // A point on the equator at zero longitude in the Earth-fixed frame is
+    // rotated by minus the Earth Rotation Angle in the inertial frame.
+    let angle = EarthRotationAngle::from_ut1_days(1.5);
+    let ground = Ecef::new(Vec3::new(6_378_137.0, 0.0, 0.0));
+    let inertial = ground.to_eci(angle);
+    let expected = Vec3::new(
+        6_378_137.0 * angle.radians_positive().cos(),
+        6_378_137.0 * angle.radians_positive().sin(),
+        0.0,
+    );
+    assert!((inertial.vector() - expected).norm() < 1e-6);
+}
+
+#[test]
+fn a_full_turn_returns_to_the_same_rotation() {
+    // The rotation depends on the angle modulo one turn, so angles that
+    // differ by whole turns give the same matrix.
+    let first = rotation_eci_to_ecef(EarthRotationAngle::from_turns(0.25));
+    let second = rotation_eci_to_ecef(EarthRotationAngle::from_turns(2.25));
+    let difference = first * second.transpose();
+    assert!(difference.is_rotation(1e-12));
+    let vector = Vec3::new(1.0, 2.0, 3.0);
+    assert!((first * vector - second * vector).norm() < EPS);
+}
+
+#[test]
+fn identity_angle_leaves_inertial_vectors_unchanged() {
+    // An angle of zero, constructed directly rather than from the epoch
+    // constant, is the identity rotation.
+    let rotation = rotation_eci_to_ecef(EarthRotationAngle::from_turns(0.0));
+    let vector = Vec3::new(1.0, 2.0, 3.0);
+    assert!((rotation * vector - vector).norm() < EPS);
+    assert!(rotation.is_rotation(1e-12));
+    assert!((rotation.determinant() - 1.0).abs() < EPS);
+}

crates/physics/apsis-frames/tests/frame.rs

@@ -0,0 +1,19 @@
+//! Tests for the frame marker types.
+
+use apsis_frames::{Ecef, Eci};
+use apsis_math::Vec3;
+
+#[test]
+fn wrappers_preserve_the_vector() {
+    let vector = Vec3::new(1.0, 2.0, 3.0);
+    assert_eq!(Eci::new(vector).vector(), vector);
+    assert_eq!(Ecef::new(vector).vector(), vector);
+}
+
+#[test]
+fn wrappers_compare_by_value() {
+    let vector = Vec3::new(-4.0, 5.5, 0.0);
+    assert_eq!(Eci::new(vector), Eci(vector));
+    assert_eq!(Ecef::new(vector), Ecef(vector));
+    assert_ne!(Eci::new(vector), Eci::new(Vec3::ZERO));
+}