physics: add frame carrying position types

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

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> {
+        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))
+    }
+}

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};

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);
+}