From 33204914560af433d00272e4b083371c57f0ef99 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 13:25:12 +0800 Subject: [PATCH] physics: add frame markers and Earth rotation angle --- crates/physics/apsis-frames/Cargo.toml | 2 + .../apsis-frames/src/earth_rotation.rs | 143 ++++++++++++++ crates/physics/apsis-frames/src/frame.rs | 41 ++++ crates/physics/apsis-frames/src/lib.rs | 53 ++++- .../apsis-frames/tests/earth_rotation.rs | 181 ++++++++++++++++++ crates/physics/apsis-frames/tests/frame.rs | 19 ++ 6 files changed, 438 insertions(+), 1 deletion(-) create mode 100644 crates/physics/apsis-frames/src/earth_rotation.rs create mode 100644 crates/physics/apsis-frames/src/frame.rs create mode 100644 crates/physics/apsis-frames/tests/earth_rotation.rs create mode 100644 crates/physics/apsis-frames/tests/frame.rs diff --git a/crates/physics/apsis-frames/Cargo.toml b/crates/physics/apsis-frames/Cargo.toml index 2f21ce1..6d9c852 100644 --- a/crates/physics/apsis-frames/Cargo.toml +++ b/crates/physics/apsis-frames/Cargo.toml @@ -11,3 +11,5 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-math.workspace = true diff --git a/crates/physics/apsis-frames/src/earth_rotation.rs b/crates/physics/apsis-frames/src/earth_rotation.rs new file mode 100644 index 0000000..262c11e --- /dev/null +++ b/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) + } +} diff --git a/crates/physics/apsis-frames/src/frame.rs b/crates/physics/apsis-frames/src/frame.rs new file mode 100644 index 0000000..cccda08 --- /dev/null +++ b/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 + } +} diff --git a/crates/physics/apsis-frames/src/lib.rs b/crates/physics/apsis-frames/src/lib.rs index 9172e2a..8f92d15 100644 --- a/crates/physics/apsis-frames/src/lib.rs +++ b/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}; diff --git a/crates/physics/apsis-frames/tests/earth_rotation.rs b/crates/physics/apsis-frames/tests/earth_rotation.rs new file mode 100644 index 0000000..492f774 --- /dev/null +++ b/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); +} diff --git a/crates/physics/apsis-frames/tests/frame.rs b/crates/physics/apsis-frames/tests/frame.rs new file mode 100644 index 0000000..7e61136 --- /dev/null +++ b/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)); +}