From efb69c77bf5bb8f0d8fddf480a31e9897dc1156e Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 13:27:05 +0800 Subject: [PATCH] physics: return plain vectors from geodetic conversions --- crates/physics/apsis-geodesy/src/geodetic.rs | 74 ++++++---------- crates/physics/apsis-geodesy/src/lib.rs | 6 +- .../physics/apsis-geodesy/tests/geodetic.rs | 84 ++++++++----------- 3 files changed, 65 insertions(+), 99 deletions(-) diff --git a/crates/physics/apsis-geodesy/src/geodetic.rs b/crates/physics/apsis-geodesy/src/geodetic.rs index 93dcb08..b0592da 100644 --- a/crates/physics/apsis-geodesy/src/geodetic.rs +++ b/crates/physics/apsis-geodesy/src/geodetic.rs @@ -1,5 +1,11 @@ -//! Geodetic positions and their conversion to Earth-centered, -//! Earth-fixed Cartesian coordinates. +//! Geodetic positions and their conversion to Earth-fixed Cartesian +//! coordinates. +//! +//! A Cartesian position here is a plain [`Vec3`] in meters. This crate +//! deliberately does not name the frame: the same conversion applies to any +//! frame whose third axis is the ellipsoid's minor axis, and the choice of +//! frame is the responsibility of the caller. `apsis-coord` supplies +//! position types that carry the frame. use crate::ellipsoid::Ellipsoid; use apsis_math::Vec3; @@ -31,13 +37,13 @@ impl Geodetic { } } - /// Converts this geodetic position to Earth-centered, Earth-fixed - /// Cartesian coordinates on the given ellipsoid. + /// Converts this geodetic position to Earth-fixed Cartesian coordinates + /// on the given ellipsoid. /// - /// The frame has its third axis along the ellipsoid's minor axis - /// through the north pole, its first axis through the intersection of - /// the equator and the prime meridian, and its second axis completing - /// a right-handed set. The result is in meters. + /// The frame has its third axis along the ellipsoid's minor axis through + /// the north pole, its first axis through the intersection of the + /// equator and the prime meridian, and its second axis completing a + /// right-handed set. The result is in meters. /// /// The closed form is /// @@ -51,8 +57,8 @@ impl Geodetic { /// where `N` is the prime vertical radius of curvature. This form is /// exact everywhere, including at the poles, where the cosine factors /// vanish. - pub fn to_ecef(self, ellipsoid: Ellipsoid) -> Ecef { - ecef_from_geodetic(self, ellipsoid) + pub fn to_cartesian(self, ellipsoid: Ellipsoid) -> Vec3 { + cartesian_from_geodetic(self, ellipsoid) } /// Returns `true` when every field is finite and the latitude and @@ -69,40 +75,11 @@ impl Geodetic { } } -/// An Earth-centered, Earth-fixed Cartesian position. +/// Converts a geodetic position to Earth-fixed Cartesian coordinates. /// -/// The frame is the one described by [`Geodetic::to_ecef`]. The wrapped -/// [`Vec3`] holds the position in meters. -#[derive(Debug, Clone, Copy, PartialEq)] -pub struct Ecef(pub Vec3); - -impl Ecef { - /// Creates an Earth-centered, Earth-fixed position from a vector in - /// meters. - pub const fn new(position: Vec3) -> Self { - Self(position) - } - - /// Returns the position in meters. - pub const fn position(self) -> Vec3 { - self.0 - } - - /// Converts this position to a geodetic position on the given - /// ellipsoid. - /// - /// Returns `None` when the position is not finite, or when it lies on - /// the polar axis so that no longitude is defined. - pub fn to_geodetic(self, ellipsoid: Ellipsoid) -> Option { - geodetic_from_ecef(self, ellipsoid) - } -} - -/// Converts a geodetic position to Earth-centered, Earth-fixed Cartesian -/// coordinates. -/// -/// See [`Geodetic::to_ecef`] for the frame definition and the closed form. -pub fn ecef_from_geodetic(geodetic: Geodetic, ellipsoid: Ellipsoid) -> Ecef { +/// See [`Geodetic::to_cartesian`] for the frame definition and the closed +/// form. +pub fn cartesian_from_geodetic(geodetic: Geodetic, ellipsoid: Ellipsoid) -> Vec3 { let Geodetic { latitude, longitude, @@ -117,11 +94,10 @@ pub fn ecef_from_geodetic(geodetic: Geodetic, ellipsoid: Ellipsoid) -> Ecef { let x = horizontal * cos_latitude * cos_longitude; let y = horizontal * cos_latitude * sin_longitude; let z = (prime_vertical * (1.0 - eccentricity_squared) + height) * sin_latitude; - Ecef(Vec3::new(x, y, z)) + Vec3::new(x, y, z) } -/// Converts Earth-centered, Earth-fixed Cartesian coordinates to a geodetic -/// position. +/// Converts Earth-fixed Cartesian coordinates to a geodetic position. /// /// The latitude and height are recovered by Bowring's method, which solves /// the meridian equation in closed form and then refines the result with a @@ -131,9 +107,9 @@ pub fn ecef_from_geodetic(geodetic: Geodetic, ellipsoid: Ellipsoid) -> Ecef { /// /// Returns `None` when the position is not finite, or when it lies on the /// polar axis so that no longitude is defined. -pub fn geodetic_from_ecef(ecef: Ecef, ellipsoid: Ellipsoid) -> Option { - let Vec3 { x, y, z } = ecef.0; - if !ecef.0.is_finite() { +pub fn geodetic_from_cartesian(cartesian: Vec3, ellipsoid: Ellipsoid) -> Option { + let Vec3 { x, y, z } = cartesian; + if !cartesian.is_finite() { return None; } diff --git a/crates/physics/apsis-geodesy/src/lib.rs b/crates/physics/apsis-geodesy/src/lib.rs index aa802e9..7534c10 100644 --- a/crates/physics/apsis-geodesy/src/lib.rs +++ b/crates/physics/apsis-geodesy/src/lib.rs @@ -20,8 +20,8 @@ //! //! // The equator has zero height and a radius equal to the semi-major axis. //! let point = Geodetic::new(0.0, 0.0, 0.0); -//! let cartesian = point.to_ecef(WGS84); -//! assert!((cartesian.position().norm() - WGS84.semi_major_axis()).abs() < 1e-9); +//! let cartesian = point.to_cartesian(WGS84); +//! assert!((cartesian.norm() - WGS84.semi_major_axis()).abs() < 1e-9); //! ``` #![forbid(unsafe_code)] @@ -31,4 +31,4 @@ mod ellipsoid; mod geodetic; pub use ellipsoid::{Ellipsoid, WGS84, WGS84_GRAVITATIONAL_PARAMETER}; -pub use geodetic::{Ecef, Geodetic, ecef_from_geodetic, geodetic_from_ecef}; +pub use geodetic::{Geodetic, cartesian_from_geodetic, geodetic_from_cartesian}; diff --git a/crates/physics/apsis-geodesy/tests/geodetic.rs b/crates/physics/apsis-geodesy/tests/geodetic.rs index 21a9eda..b1dc986 100644 --- a/crates/physics/apsis-geodesy/tests/geodetic.rs +++ b/crates/physics/apsis-geodesy/tests/geodetic.rs @@ -1,6 +1,6 @@ -//! Tests for geodetic and Earth-centered, Earth-fixed conversions. +//! Tests for geodetic and Cartesian conversions. -use apsis_geodesy::{Ecef, Geodetic, WGS84}; +use apsis_geodesy::{Geodetic, WGS84, geodetic_from_cartesian}; use apsis_math::Vec3; /// Tolerance for position comparisons, in meters. @@ -14,35 +14,35 @@ const ANGLE_TOL: f64 = 1e-11; #[test] fn equatorial_point_lies_on_the_semi_major_axis() { - let point = Geodetic::new(0.0, 0.0, 0.0).to_ecef(WGS84); - assert!((point.position().x - WGS84.semi_major_axis()).abs() < POSITION_TOL); - assert!(point.position().y.abs() < POSITION_TOL); - assert!(point.position().z.abs() < POSITION_TOL); + let point = Geodetic::new(0.0, 0.0, 0.0).to_cartesian(WGS84); + assert!((point.x - WGS84.semi_major_axis()).abs() < POSITION_TOL); + assert!(point.y.abs() < POSITION_TOL); + assert!(point.z.abs() < POSITION_TOL); } #[test] fn pole_lies_on_the_semi_minor_axis() { - let point = Geodetic::new(core::f64::consts::FRAC_PI_2, 0.0, 0.0).to_ecef(WGS84); - assert!(point.position().x.abs() < POSITION_TOL); - assert!(point.position().y.abs() < POSITION_TOL); - assert!((point.position().z - WGS84.semi_minor_axis()).abs() < POSITION_TOL); + let point = Geodetic::new(core::f64::consts::FRAC_PI_2, 0.0, 0.0).to_cartesian(WGS84); + assert!(point.x.abs() < POSITION_TOL); + assert!(point.y.abs() < POSITION_TOL); + assert!((point.z - WGS84.semi_minor_axis()).abs() < POSITION_TOL); } #[test] fn longitude_quarter_turn_places_the_point_on_the_second_axis() { - let point = Geodetic::new(0.0, core::f64::consts::FRAC_PI_2, 0.0).to_ecef(WGS84); - assert!(point.position().x.abs() < POSITION_TOL); - assert!((point.position().y - WGS84.semi_major_axis()).abs() < POSITION_TOL); + let point = Geodetic::new(0.0, core::f64::consts::FRAC_PI_2, 0.0).to_cartesian(WGS84); + assert!(point.x.abs() < POSITION_TOL); + assert!((point.y - WGS84.semi_major_axis()).abs() < POSITION_TOL); } #[test] fn height_moves_the_point_along_the_ellipsoid_normal() { - let surface = Geodetic::new(0.6, 0.3, 0.0).to_ecef(WGS84); - let elevated = Geodetic::new(0.6, 0.3, 500.0).to_ecef(WGS84); - assert!((surface.position().distance(elevated.position()) - 500.0).abs() < POSITION_TOL); + let surface = Geodetic::new(0.6, 0.3, 0.0).to_cartesian(WGS84); + let elevated = Geodetic::new(0.6, 0.3, 500.0).to_cartesian(WGS84); + assert!((surface.distance(elevated) - 500.0).abs() < POSITION_TOL); // The normal at a non-polar latitude is not vertical, so the height // changes the horizontal components as well as the third one. - assert!(elevated.position().x > surface.position().x); + assert!(elevated.x > surface.x); } #[test] @@ -55,9 +55,8 @@ fn conversion_round_trips_over_a_global_grid() { longitude_degrees.to_radians(), height, ); - let ecef = original.to_ecef(WGS84); - let recovered = ecef - .to_geodetic(WGS84) + let cartesian = original.to_cartesian(WGS84); + let recovered = geodetic_from_cartesian(cartesian, WGS84) .expect("off-axis points have a longitude"); assert!( (recovered.latitude - original.latitude).abs() < ANGLE_TOL, @@ -80,28 +79,27 @@ fn conversion_round_trips_over_a_global_grid() { #[test] fn round_trip_preserves_the_cartesian_position() { - let original = Geodetic::new(0.921, -1.234, 12_345.0).to_ecef(WGS84); - let recovered = original - .to_geodetic(WGS84) + let original = Geodetic::new(0.921, -1.234, 12_345.0).to_cartesian(WGS84); + let recovered = geodetic_from_cartesian(original, WGS84) .expect("off-axis points have a longitude") - .to_ecef(WGS84); - assert!(original.position().distance(recovered.position()) < POSITION_TOL); + .to_cartesian(WGS84); + assert!(original.distance(recovered) < POSITION_TOL); } #[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); + let on_axis = Vec3::new(0.0, 0.0, WGS84.semi_minor_axis()); + assert_eq!(geodetic_from_cartesian(on_axis, WGS84), None); // The origin is also on the axis. - assert_eq!(Ecef::new(Vec3::ZERO).to_geodetic(WGS84), None); + assert_eq!(geodetic_from_cartesian(Vec3::ZERO, WGS84), None); } #[test] fn non_finite_positions_are_rejected() { - let invalid = Ecef::new(Vec3::new(f64::NAN, 0.0, 6_000_000.0)); - assert_eq!(invalid.to_geodetic(WGS84), None); - let infinite = Ecef::new(Vec3::new(f64::INFINITY, 0.0, 0.0)); - assert_eq!(infinite.to_geodetic(WGS84), None); + let invalid = Vec3::new(f64::NAN, 0.0, 6_000_000.0); + assert_eq!(geodetic_from_cartesian(invalid, WGS84), None); + let infinite = Vec3::new(f64::INFINITY, 0.0, 0.0); + assert_eq!(geodetic_from_cartesian(infinite, WGS84), None); } #[test] @@ -114,21 +112,21 @@ fn known_point_matches_an_independent_computation() { -121.333_333_333_333_33_f64.to_radians(), 0.0, ) - .to_ecef(WGS84); + .to_cartesian(WGS84); assert!( - (point.position().x + 2_634_612.181_141_605_6).abs() < 1e-6, + (point.x + 2_634_612.181_141_605_6).abs() < 1e-6, "x = {}", - point.position().x + point.x ); assert!( - (point.position().y + 4_327_502.253_728_065).abs() < 1e-6, + (point.y + 4_327_502.253_728_065).abs() < 1e-6, "y = {}", - point.position().y + point.y ); assert!( - (point.position().z - 3_861_564.098_809_606_4).abs() < 1e-6, + (point.z - 3_861_564.098_809_606_4).abs() < 1e-6, "z = {}", - point.position().z + point.z ); } @@ -140,11 +138,3 @@ fn validity_is_reported() { assert!(!Geodetic::new(0.0, 4.0, 0.0).is_valid()); assert!(!Geodetic::new(0.0, 0.0, f64::NAN).is_valid()); } - -#[test] -fn ecef_wrapper_exposes_its_vector() { - let vector = Vec3::new(1.0, 2.0, 3.0); - let ecef = Ecef::new(vector); - assert_eq!(ecef.position(), vector); - assert_eq!(ecef, Ecef(vector)); -}