proj_core/
ellipsoid.rs

1//! Reference ellipsoid definitions.
2
3use crate::error::{Error, Result};
4use std::f64::consts::PI;
5
6/// A reference ellipsoid.
7#[derive(Debug, Clone, Copy)]
8pub struct Ellipsoid {
9    /// Semi-major axis (equatorial radius) in meters.
10    a: f64,
11    /// Flattening (f = (a - b) / a).
12    f: f64,
13}
14
15impl Ellipsoid {
16    /// Create an ellipsoid from semi-major axis and inverse flattening.
17    pub fn from_a_rf(a: f64, rf: f64) -> Result<Self> {
18        if !rf.is_finite() || rf <= 1.0 {
19            return Err(Error::InvalidDefinition(
20                "inverse flattening must be a finite number greater than 1".into(),
21            ));
22        }
23        Self::from_a_f(a, 1.0 / rf)
24    }
25
26    /// Create an ellipsoid from semi-major axis and flattening.
27    pub fn from_a_f(a: f64, f: f64) -> Result<Self> {
28        if !a.is_finite() || a <= 0.0 {
29            return Err(Error::InvalidDefinition(
30                "semi-major axis must be a finite positive number".into(),
31            ));
32        }
33        if !f.is_finite() || !(0.0..1.0).contains(&f) {
34            return Err(Error::InvalidDefinition(
35                "flattening must be finite and in the range [0, 1)".into(),
36            ));
37        }
38        Ok(Self { a, f })
39    }
40
41    /// Create a sphere (flattening = 0).
42    pub fn sphere(radius: f64) -> Result<Self> {
43        Self::from_a_f(radius, 0.0)
44    }
45
46    const fn from_a_rf_unchecked(a: f64, rf: f64) -> Self {
47        Self { a, f: 1.0 / rf }
48    }
49
50    const fn sphere_unchecked(radius: f64) -> Self {
51        Self { a: radius, f: 0.0 }
52    }
53
54    /// Semi-major axis (equatorial radius) in meters.
55    pub const fn semi_major_axis(&self) -> f64 {
56        self.a
57    }
58
59    /// Flattening (f = (a - b) / a).
60    pub const fn flattening(&self) -> f64 {
61        self.f
62    }
63
64    /// Inverse flattening, or 0 for a sphere.
65    pub fn inverse_flattening(&self) -> f64 {
66        if self.f == 0.0 {
67            0.0
68        } else {
69            1.0 / self.f
70        }
71    }
72
73    /// Semi-minor axis (polar radius) in meters.
74    pub fn b(&self) -> f64 {
75        self.a * (1.0 - self.f)
76    }
77
78    /// First eccentricity squared: e² = 2f - f².
79    pub fn e2(&self) -> f64 {
80        2.0 * self.f - self.f * self.f
81    }
82
83    /// First eccentricity.
84    pub fn e(&self) -> f64 {
85        self.e2().sqrt()
86    }
87
88    /// Third flattening: n = f / (2 - f) = (a - b) / (a + b).
89    pub fn n(&self) -> f64 {
90        self.f / (2.0 - self.f)
91    }
92
93    /// Second eccentricity squared: e'² = e² / (1 - e²).
94    pub fn ep2(&self) -> f64 {
95        let e2 = self.e2();
96        e2 / (1.0 - e2)
97    }
98
99    /// Radius of curvature in the prime vertical at latitude phi (radians).
100    pub fn n_radius(&self, phi: f64) -> f64 {
101        let sin_phi = phi.sin();
102        self.a / (1.0 - self.e2() * sin_phi * sin_phi).sqrt()
103    }
104
105    /// Radius of curvature in the meridian at latitude phi (radians).
106    pub fn m_radius(&self, phi: f64) -> f64 {
107        let e2 = self.e2();
108        let sin_phi = phi.sin();
109        let denom = 1.0 - e2 * sin_phi * sin_phi;
110        self.a * (1.0 - e2) / denom.powf(1.5)
111    }
112}
113
114// Well-known ellipsoids.
115
116/// WGS 84 ellipsoid.
117pub const WGS84: Ellipsoid = Ellipsoid::from_a_rf_unchecked(6378137.0, 298.257223563);
118
119/// GRS 1980 ellipsoid (used by NAD83).
120pub const GRS80: Ellipsoid = Ellipsoid::from_a_rf_unchecked(6378137.0, 298.257222101);
121
122/// Clarke 1866 ellipsoid (used by NAD27).
123pub const CLARKE1866: Ellipsoid = Ellipsoid::from_a_rf_unchecked(6378206.4, 294.978698214);
124
125/// International 1924 (Hayford) ellipsoid.
126pub const INTL1924: Ellipsoid = Ellipsoid::from_a_rf_unchecked(6378388.0, 297.0);
127
128/// Bessel 1841 ellipsoid.
129pub const BESSEL1841: Ellipsoid = Ellipsoid::from_a_rf_unchecked(6377397.155, 299.1528128);
130
131/// Krassowsky 1940 ellipsoid.
132pub const KRASSOWSKY: Ellipsoid = Ellipsoid::from_a_rf_unchecked(6378245.0, 298.3);
133
134/// Airy 1830 ellipsoid (used by OSGB36).
135pub const AIRY1830: Ellipsoid = Ellipsoid::from_a_rf_unchecked(6377563.396, 299.3249646);
136
137/// Unit sphere (radius = 1).
138pub const UNIT_SPHERE: Ellipsoid = Ellipsoid::sphere_unchecked(1.0);
139
140/// Authalic sphere matching WGS84 surface area.
141pub const WGS84_SPHERE: Ellipsoid = Ellipsoid::sphere_unchecked(6371007.181);
142
143/// Convert degrees to radians.
144pub fn deg_to_rad(deg: f64) -> f64 {
145    deg * PI / 180.0
146}
147
148/// Convert radians to degrees.
149pub fn rad_to_deg(rad: f64) -> f64 {
150    rad * 180.0 / PI
151}
152
153#[cfg(test)]
154mod tests {
155    use super::*;
156
157    #[test]
158    fn wgs84_basics() {
159        let e = WGS84;
160        assert!((e.semi_major_axis() - 6378137.0).abs() < 1e-6);
161        assert!((e.b() - 6356752.314245).abs() < 1.0);
162        assert!((e.e2() - 0.00669437999014).abs() < 1e-11);
163    }
164
165    #[test]
166    fn sphere_has_zero_eccentricity() {
167        let s = UNIT_SPHERE;
168        assert_eq!(s.e2(), 0.0);
169        assert_eq!(s.b(), 1.0);
170    }
171
172    #[test]
173    fn deg_rad_roundtrip() {
174        let deg = 45.0;
175        assert!((rad_to_deg(deg_to_rad(deg)) - deg).abs() < 1e-12);
176    }
177
178    #[test]
179    fn reject_invalid_ellipsoid_dimensions() {
180        assert!(Ellipsoid::from_a_rf(f64::NAN, 298.257223563).is_err());
181        assert!(Ellipsoid::from_a_rf(6378137.0, 0.0).is_err());
182        assert!(Ellipsoid::from_a_f(6378137.0, 1.0).is_err());
183        assert!(Ellipsoid::sphere(-1.0).is_err());
184    }
185}
proj_core/ellipsoid.rs

proj_core/
ellipsoid.rs
