From 74df22d3daa5a7141244b31ae25851d108dc2612 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 13:27:05 +0800 Subject: [PATCH] physics: add frame carrying position types --- crates/physics/apsis-coord/Cargo.toml | 4 + crates/physics/apsis-coord/src/frames.rs | 121 +++++++++++++++++++++ crates/physics/apsis-coord/src/lib.rs | 52 ++++++++- crates/physics/apsis-coord/tests/frames.rs | 93 ++++++++++++++++ 4 files changed, 269 insertions(+), 1 deletion(-) create mode 100644 crates/physics/apsis-coord/src/frames.rs create mode 100644 crates/physics/apsis-coord/tests/frames.rs diff --git a/crates/physics/apsis-coord/Cargo.toml b/crates/physics/apsis-coord/Cargo.toml index 04efaf0..32c0476 100644 --- a/crates/physics/apsis-coord/Cargo.toml +++ b/crates/physics/apsis-coord/Cargo.toml @@ -11,3 +11,7 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-math.workspace = true +apsis-geodesy.workspace = true +apsis-frames.workspace = true diff --git a/crates/physics/apsis-coord/src/frames.rs b/crates/physics/apsis-coord/src/frames.rs new file mode 100644 index 0000000..5794756 --- /dev/null +++ b/crates/physics/apsis-coord/src/frames.rs @@ -0,0 +1,121 @@ +//! Positions that carry their reference frame in the type. + +use apsis_frames::{EarthRotationAngle, rotation_ecef_to_eci, rotation_eci_to_ecef}; +use apsis_geodesy::{ + Ellipsoid, Geodetic as EllipsoidalPosition, WGS84, cartesian_from_geodetic, + geodetic_from_cartesian, +}; +use apsis_math::Vec3; + +/// A Cartesian position in the Earth-centered inertial frame. +/// +/// The frame is aligned with the mean equator and equinox of J2000.0 and is +/// treated as inertial, so it does not rotate with the Earth. Coordinates +/// are in meters. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Eci(Vec3); + +/// A Cartesian position 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. Coordinates are in meters. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Ecef(Vec3); + +/// A position in the Earth-fixed frame given as latitude, longitude, and +/// height above the WGS 84 ellipsoid. +/// +/// Latitude and longitude are in radians and height is in meters. Latitude +/// is geodetic, measured between the ellipsoid normal and the equatorial +/// plane. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Geodetic(EllipsoidalPosition); + +impl Eci { + /// Creates an inertial position from a vector in meters. + pub const fn new(vector: Vec3) -> Self { + Self(vector) + } + + /// Returns the position in meters. + pub const fn vector(self) -> Vec3 { + self.0 + } + + /// Converts this inertial position to the Earth-fixed frame at the given + /// Earth Rotation Angle. + pub fn to_ecef(self, angle: EarthRotationAngle) -> Ecef { + Ecef(rotation_eci_to_ecef(angle) * self.0) + } +} + +impl Ecef { + /// Creates an Earth-fixed position from a vector in meters. + pub const fn new(vector: Vec3) -> Self { + Self(vector) + } + + /// Returns the position in meters. + pub const fn vector(self) -> Vec3 { + self.0 + } + + /// Converts this Earth-fixed position to the inertial frame at the given + /// Earth Rotation Angle. + pub fn to_eci(self, angle: EarthRotationAngle) -> Eci { + Eci(rotation_ecef_to_eci(angle) * self.0) + } + + /// Converts this Earth-fixed position to latitude, longitude, and height + /// above the given ellipsoid. + /// + /// Returns `None` when the position lies on the polar axis, where no + /// longitude is defined. + pub fn to_geodetic(self, ellipsoid: Ellipsoid) -> Option { + geodetic_from_cartesian(self.0, ellipsoid).map(Geodetic) + } +} + +impl Geodetic { + /// Creates a geodetic position from latitude, longitude, and height. + /// + /// Latitude and longitude are in radians and height is in meters. + pub const fn new(latitude: f64, longitude: f64, height: f64) -> Self { + Self(EllipsoidalPosition::new(latitude, longitude, height)) + } + + /// Returns the geodetic latitude in radians. + pub const fn latitude(self) -> f64 { + self.0.latitude + } + + /// Returns the longitude in radians, positive east. + pub const fn longitude(self) -> f64 { + self.0.longitude + } + + /// Returns the height above the ellipsoid in meters. + pub const fn height(self) -> f64 { + self.0.height + } + + /// Returns the underlying ellipsoidal position. + pub const fn ellipsoidal(self) -> EllipsoidalPosition { + self.0 + } + + /// Converts this position to the Earth-fixed frame on the WGS 84 + /// ellipsoid. + /// + /// This is the common case. Use [`Geodetic::to_ecef_on`] for a different + /// ellipsoid. + pub fn to_ecef(self) -> Ecef { + self.to_ecef_on(WGS84) + } + + /// Converts this position to the Earth-fixed frame on the given + /// ellipsoid. + pub fn to_ecef_on(self, ellipsoid: Ellipsoid) -> Ecef { + Ecef(cartesian_from_geodetic(self.0, ellipsoid)) + } +} diff --git a/crates/physics/apsis-coord/src/lib.rs b/crates/physics/apsis-coord/src/lib.rs index a562cea..c235775 100644 --- a/crates/physics/apsis-coord/src/lib.rs +++ b/crates/physics/apsis-coord/src/lib.rs @@ -1 +1,51 @@ -//! Coordinate frame markers and Earth-centered inertial, Earth-centered Earth-fixed, and geodetic transformations. +//! Coordinate frames and the transformations between them. +//! +//! This crate defines positions that carry their reference frame in the +//! type, so a position expressed in one frame cannot be passed where another +//! is expected. The two frames are the Earth-centered inertial (ECI) frame +//! and the Earth-centered Earth-fixed (ECEF) frame. +//! +//! Three position representations are offered: +//! +//! - [`Eci`], an inertial Cartesian position in meters +//! - [`Ecef`], an Earth-fixed Cartesian position in meters +//! - [`Geodetic`], a latitude, longitude, and height above the WGS 84 +//! ellipsoid +//! +//! The transformation between [`Eci`] and [`Ecef`] is the rotation defined +//! by the Earth Rotation Angle, which `apsis-frames` computes. The +//! transformation between [`Ecef`] and [`Geodetic`] uses the WGS 84 +//! ellipsoid from `apsis-geodesy`. +//! +//! ``` +//! use apsis_coord::{Ecef, Eci, Geodetic, WGS84}; +//! use apsis_frames::EarthRotationAngle; +//! use apsis_math::Vec3; +//! +//! // A station on the equator at zero longitude. +//! let station = Geodetic::new(0.0, 0.0, 0.0); +//! let fixed = station.to_ecef(); +//! assert!((fixed.vector().x - WGS84.semi_major_axis()).abs() < 1e-6); +//! +//! // Rotating it into the inertial frame preserves its length. +//! let inertial = fixed.to_eci(EarthRotationAngle::from_ut1_days(0.0)); +//! assert!((inertial.vector().norm() - fixed.vector().norm()).abs() < 1e-6); +//! let _ = Eci::new(Vec3::X); +//! ``` +//! +//! Sources: +//! +//! - Earth Rotation Angle: *IERS Conventions (2010)*, equation 8.32, as +//! implemented by `apsis-frames`. +//! - Geodetic conversions and the WGS 84 ellipsoid: NIMA TR8350.2, third +//! edition, 2000, as implemented by `apsis-geodesy`. + +#![forbid(unsafe_code)] +#![deny(missing_docs)] + +mod frames; + +pub use frames::{Ecef, Eci, Geodetic}; + +pub use apsis_frames::EarthRotationAngle; +pub use apsis_geodesy::{Ellipsoid, WGS84}; diff --git a/crates/physics/apsis-coord/tests/frames.rs b/crates/physics/apsis-coord/tests/frames.rs new file mode 100644 index 0000000..a0537be --- /dev/null +++ b/crates/physics/apsis-coord/tests/frames.rs @@ -0,0 +1,93 @@ +//! Tests for the frame carrying position types. + +use apsis_coord::{EarthRotationAngle, Ecef, Eci, Geodetic, WGS84}; +use apsis_math::Vec3; + +const EPS: f64 = 1e-9; + +#[test] +fn constructors_round_trip_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 geodetic_fields_round_trip() { + let position = Geodetic::new(0.5, -1.0, 250.0); + assert_eq!(position.latitude(), 0.5); + assert_eq!(position.longitude(), -1.0); + assert_eq!(position.height(), 250.0); + assert_eq!(position.ellipsoidal().latitude, 0.5); +} + +#[test] +fn geodetic_to_ecef_uses_wgs84_by_default() { + let position = Geodetic::new(0.35, -1.2, 123.0); + assert_eq!(position.to_ecef(), position.to_ecef_on(WGS84)); +} + +#[test] +fn geodetic_to_ecef_and_back_round_trips() { + let original = Geodetic::new(0.921, -1.234, 12_345.0); + let recovered = original + .to_ecef() + .to_geodetic(WGS84) + .expect("off-axis points have a longitude"); + assert!((recovered.latitude() - original.latitude()).abs() < 1e-11); + assert!((recovered.longitude() - original.longitude()).abs() < 1e-11); + assert!((recovered.height() - original.height()).abs() < 1e-4); +} + +#[test] +fn polar_axis_has_no_defined_longitude() { + let on_axis = Ecef::new(Vec3::new(0.0, 0.0, WGS84.semi_minor_axis())); + assert_eq!(on_axis.to_geodetic(WGS84), None); + assert_eq!(Ecef::new(Vec3::ZERO).to_geodetic(WGS84), None); +} + +#[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() < EPS); +} + +#[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 a_ground_point_rotates_with_the_earth() { + // A point on the equator at zero longitude in the Earth-fixed frame is + // rotated by 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 radians = angle.radians_positive(); + let expected = Vec3::new( + 6_378_137.0 * radians.cos(), + 6_378_137.0 * radians.sin(), + 0.0, + ); + assert!((inertial.vector() - expected).norm() < 1e-6); +} + +#[test] +fn a_geodetic_station_round_trips_through_both_frames() { + // A station is fixed to the ground, so converting it to the inertial + // frame and back at the same angle returns the original position. + let station = Geodetic::new(0.662, 0.175, 80.0); + let angle = EarthRotationAngle::from_ut1_julian_date(2_460_000.5); + let inertial = station.to_ecef().to_eci(angle); + let back = inertial.to_ecef(angle).to_geodetic(WGS84).expect("station"); + assert!((back.latitude() - station.latitude()).abs() < 1e-11); + assert!((back.longitude() - station.longitude()).abs() < 1e-11); + assert!((back.height() - station.height()).abs() < 1e-4); +}