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}