Skip to main content

geometry_cs/
spheroid.rs

1//! Reference ellipsoid for geographic coordinate systems.
2//!
3//! Mirrors `boost::geometry::srs::spheroid<RadiusType>` from
4//! `boost/geometry/srs/spheroid.hpp`. Boost stores the spheroid as a
5//! pair of radii (`m_a`, `m_b`) and exposes them via a
6//! `get_radius<I>()` template; we store the *flattening* alongside the
7//! equatorial radius instead, because every geodesic formula
8//! (Andoyer, Vincenty, Thomas) needs `f` directly — keeping it as the
9//! primary stored field avoids the round-trip through `(a − b) / a`
10//! that Boost's strategies do at every call site.
11//!
12//! `WGS84` is the only reference we ship in this task; other ellipsoids
13//! (GRS80, IUGG67, …) land alongside the geographic strategies in T42.
14
15/// Reference ellipsoid for geographic coordinate systems.
16///
17/// Field semantics match Boost's `srs::spheroid<T>`
18/// (`boost/geometry/srs/spheroid.hpp:49-112`), except the second
19/// stored value is the *flattening* instead of the polar radius — the
20/// two forms are interconvertible via [`polar_radius`] /
21/// [`flattening`] but `f` is what the geodesic formulas actually take
22/// as input.
23///
24/// [`polar_radius`]: Spheroid::polar_radius
25/// [`flattening`]: Spheroid::flattening
26///
27/// # Examples
28///
29/// ```
30/// use geometry_cs::Spheroid;
31/// let s = Spheroid::WGS84;
32/// assert_eq!(s.equatorial_radius, 6_378_137.0);
33/// // The WGS84 polar radius is ~6 356 752.3 m.
34/// assert!((s.polar_radius() - 6_356_752.314_245).abs() < 0.01);
35/// ```
36#[derive(Debug, Clone, Copy, PartialEq)]
37pub struct Spheroid {
38    /// Equatorial (semi-major) radius `a`, in metres.
39    pub equatorial_radius: f64,
40    /// Flattening `f = (a − b) / a`, dimensionless.
41    pub flattening: f64,
42}
43
44impl Spheroid {
45    /// The WGS84 reference ellipsoid.
46    ///
47    /// Equatorial radius and flattening per the WGS84 defining
48    /// parameters; matches the default-constructed
49    /// `srs::spheroid<RadiusType>` in
50    /// `boost/geometry/srs/spheroid.hpp:62-69`, which seeds
51    /// `m_a = 6_378_137.0` and `m_b = 6_356_752.314_245_179_3`.
52    pub const WGS84: Self = Self {
53        equatorial_radius: 6_378_137.0,
54        flattening: 1.0 / 298.257_223_563,
55    };
56
57    /// Semi-minor axis `b = a · (1 − f)`, in metres.
58    ///
59    /// Counterpart to `srs::spheroid<T>::get_radius<2>()` in
60    /// `boost/geometry/srs/spheroid.hpp:78-92`.
61    #[inline]
62    #[must_use]
63    pub fn polar_radius(&self) -> f64 {
64        self.equatorial_radius * (1.0 - self.flattening)
65    }
66
67    /// First eccentricity squared `e² = 2f − f²`.
68    ///
69    /// Used by the geodesic strategies (Andoyer / Vincenty / Thomas)
70    /// in `boost/geometry/strategies/geographic/*.hpp`. The identity
71    /// `e² = 2f − f²` follows from `b = a(1 − f)` and
72    /// `e² = (a² − b²) / a²`.
73    #[inline]
74    #[must_use]
75    pub fn eccentricity_squared(&self) -> f64 {
76        2.0 * self.flattening - self.flattening * self.flattening
77    }
78}