rustial_engine/math/

ellipsoid.rs

//! Reference ellipsoid definitions and derived geometric constants.
//!
//! # What is a reference ellipsoid?
//!
//! A reference ellipsoid is a mathematically defined surface that
//! approximates the shape of the Earth (or another body).  It is
//! defined by two independent parameters:
//!
//! - **Semi-major axis `a`** (equatorial radius), in meters.
//! - **Inverse flattening `1/f`** (the standard geodetic convention).
//!
//! All other geometric constants (`b`, `e2`, `e'2`) are derived from
//! these two.  See [`Ellipsoid::from_a_inv_f`] for the formulae.
//!
//! # Pre-built constants
//!
//! | Constant | Standard | Usage |
//! |----------|----------|-------|
//! | [`Ellipsoid::WGS84`] | NIMA TR8350.2 / EPSG:4326 | GPS, global mapping |
//! | [`Ellipsoid::GRS80`] | IUGG 1980 / EPSG:7019 | NAD83, ETRS89 |
//!
//! # Consumers
//!
//! This type is used by [`geodesic_distance_on`](crate::geodesic_distance_on),
//! [`geodesic_destination_on`](crate::geodesic_destination_on),
//! [`WebMercator`](crate::WebMercator), and
//! [`Equirectangular`](crate::Equirectangular).

/// A reference ellipsoid defined by its semi-major axis and flattening.
///
/// All four fields are stored for fast access, but only two are
/// independent (`a` and `f`).  The derived fields `b` and `e2` are
/// computed at construction time (in `const` context) to avoid
/// repeated work in hot paths.
///
/// # Field summary
///
/// | Field | Formula | Description |
/// |-------|---------|-------------|
/// | `a` | (input) | Semi-major axis (equatorial radius), meters |
/// | `f` | `1 / inv_f` | Flattening |
/// | `b` | `a * (1 - f)` | Semi-minor axis (polar radius), meters |
/// | `e2` | `2f - f^2` | First eccentricity squared |
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Ellipsoid {
    /// Semi-major axis (equatorial radius) in meters.
    pub a: f64,
    /// Flattening: `f = (a - b) / a`.
    pub f: f64,
    /// Semi-minor axis (polar radius) in meters: `b = a * (1 - f)`.
    pub b: f64,
    /// First eccentricity squared: `e2 = 2f - f^2`.
    pub e2: f64,
}

impl Ellipsoid {
    /// Construct an ellipsoid from semi-major axis `a` (meters) and
    /// inverse flattening `inv_f` (dimensionless).
    ///
    /// Inverse flattening is the standard geodetic convention because
    /// `f` itself is a very small number (~0.003 for Earth) and
    /// `1/f` (~298.257) is easier to read and compare across standards.
    ///
    /// # Derived quantities
    ///
    /// ```text
    /// f  = 1 / inv_f
    /// b  = a * (1 - f)
    /// e2 = 2f - f^2
    /// ```
    #[inline]
    pub const fn from_a_inv_f(a: f64, inv_f: f64) -> Self {
        let f = 1.0 / inv_f;
        let b = a * (1.0 - f);
        let e2 = 2.0 * f - f * f;
        Self { a, f, b, e2 }
    }

    /// WGS-84 ellipsoid (EPSG:4326, NIMA TR8350.2).
    ///
    /// The World Geodetic System 1984 is the reference frame used by
    /// GPS and the default for all rustial coordinate operations.
    ///
    /// - `a  = 6 378 137.0 m`
    /// - `1/f = 298.257 223 563`
    pub const WGS84: Self = Self::from_a_inv_f(6_378_137.0, 298.257_223_563);

    /// GRS-80 ellipsoid (EPSG:7019, IUGG 1980).
    ///
    /// The Geodetic Reference System 1980 is used by NAD83 (North
    /// America) and ETRS89 (Europe).  It is nearly identical to WGS-84
    /// (difference in `1/f` is ~0.0001) but is a distinct definition.
    ///
    /// - `a  = 6 378 137.0 m`
    /// - `1/f = 298.257 222 101`
    pub const GRS80: Self = Self::from_a_inv_f(6_378_137.0, 298.257_222_101);

    /// Second eccentricity squared: `e'2 = (a^2 - b^2) / b^2`.
    ///
    /// Used in Vincenty geodesic formulae and some map projections.
    #[inline]
    pub fn second_eccentricity_squared(&self) -> f64 {
        (self.a * self.a - self.b * self.b) / (self.b * self.b)
    }

    /// Radius of curvature in the prime vertical (N) at geodetic
    /// latitude `lat_rad` (radians).
    ///
    /// `N = a / sqrt(1 - e2 * sin^2(lat))`
    ///
    /// This is the distance from the surface to the polar axis along
    /// the ellipsoid normal.  Used in geographic-to-ECEF conversion
    /// and UTM projection.
    #[inline]
    pub fn radius_of_curvature_n(&self, lat_rad: f64) -> f64 {
        let sin_lat = lat_rad.sin();
        self.a / (1.0 - self.e2 * sin_lat * sin_lat).sqrt()
    }

    /// Whether this ellipsoid is effectively a sphere (`f == 0`).
    ///
    /// Can be used to select simplified (spherical) code paths where
    /// the ellipsoidal correction is negligible.
    #[inline]
    pub fn is_sphere(&self) -> bool {
        self.f == 0.0
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    // -- WGS-84 reference values ------------------------------------------

    #[test]
    fn wgs84_constants() {
        let e = Ellipsoid::WGS84;
        assert!((e.a - 6_378_137.0).abs() < 1e-6);
        assert!((e.b - 6_356_752.314_245_179).abs() < 1e-3);
        assert!((e.f - 1.0 / 298.257_223_563).abs() < 1e-15);
        assert!((e.e2 - 0.006_694_379_990_141_316).abs() < 1e-12);
    }

    #[test]
    fn second_eccentricity() {
        let e = Ellipsoid::WGS84;
        let ep2 = e.second_eccentricity_squared();
        assert!((ep2 - 0.006_739_496_742_276_434).abs() < 1e-12);
    }

    // -- GRS-80 reference values ------------------------------------------

    #[test]
    fn grs80_constants() {
        let e = Ellipsoid::GRS80;
        assert!((e.a - 6_378_137.0).abs() < 1e-6);
        assert!((e.f - 1.0 / 298.257_222_101).abs() < 1e-15);
        // GRS-80 and WGS-84 share the same `a` but differ in `f`.
        assert!(e.f != Ellipsoid::WGS84.f);
        // Difference is tiny (~1e-10).
        assert!((e.f - Ellipsoid::WGS84.f).abs() < 1e-8);
    }

    // -- Custom ellipsoid -------------------------------------------------

    #[test]
    fn custom_ellipsoid_roundtrip() {
        let e = Ellipsoid::from_a_inv_f(6_400_000.0, 300.0);
        assert!((e.a - 6_400_000.0).abs() < 1e-6);
        assert!((e.f - 1.0 / 300.0).abs() < 1e-15);
        assert!((e.b - e.a * (1.0 - e.f)).abs() < 1e-6);
        assert!((e.e2 - (2.0 * e.f - e.f * e.f)).abs() < 1e-15);
    }

    // -- Radius of curvature N --------------------------------------------

    #[test]
    fn radius_n_at_equator() {
        let e = Ellipsoid::WGS84;
        // At the equator (lat=0), N = a.
        let n = e.radius_of_curvature_n(0.0);
        assert!((n - e.a).abs() < 1e-6);
    }

    #[test]
    fn radius_n_at_pole() {
        let e = Ellipsoid::WGS84;
        // At the pole (lat=90 deg), N = a / sqrt(1 - e2) = a^2 / b.
        let n = e.radius_of_curvature_n(std::f64::consts::FRAC_PI_2);
        let expected = e.a * e.a / e.b;
        assert!((n - expected).abs() < 1e-3);
    }

    // -- Sphere detection -------------------------------------------------

    #[test]
    fn wgs84_is_not_sphere() {
        assert!(!Ellipsoid::WGS84.is_sphere());
    }

    #[test]
    fn perfect_sphere() {
        // inv_f = infinity => f = 0 => sphere.
        let s = Ellipsoid::from_a_inv_f(6_371_000.0, f64::INFINITY);
        assert!(s.is_sphere());
        assert!((s.a - s.b).abs() < 1e-6);
        assert!(s.e2.abs() < 1e-15);
    }
}