jord/ellipsoidal/
ellipsoid.rs

1use crate::{
2    surface::Surface, Angle, Cartesian3DVector, GeocentricPosition, GeodeticPosition, Length,
3    NVector, Vec3,
4};
5
6/// An ellipsoid.
7#[derive(PartialEq, Clone, Copy, Debug, Default)]
8#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] // codecov:ignore:this
9pub struct Ellipsoid {
10    equatorial_radius: Length,
11    polar_radius: Length,
12    eccentricity: f64,
13    flattening: f64,
14}
15
16impl Ellipsoid {
17    /// [World Geodetic](https://en.wikipedia.org/wiki/World_Geodetic_System) 84 Ellipsoid.
18    pub const WGS84: Ellipsoid = Ellipsoid {
19        equatorial_radius: Length::from_metres(6_378_137.0f64),
20        polar_radius: Length::from_metres(6_356_752.314245179f64),
21        eccentricity: 0.08181919084262157f64,
22        flattening: 0.0033528106647474805f64,
23    };
24
25    /// Geodetic Reference System 1980 Ellipsoid.
26    pub const GRS80: Ellipsoid = Ellipsoid {
27        equatorial_radius: Length::from_metres(6_378_137.0f64),
28        polar_radius: Length::from_metres(6_356_752.314140356f64),
29        eccentricity: 0.08181919104281514f64,
30        flattening: 0.003352810681182319f64,
31    };
32
33    /// [World Geodetic](https://en.wikipedia.org/wiki/World_Geodetic_System) 72 Ellipsoid.
34    pub const WGS72: Ellipsoid = Ellipsoid {
35        equatorial_radius: Length::from_metres(6_378_135.0f64),
36        polar_radius: Length::from_metres(6_356_750.520016094f64),
37        eccentricity: 0.08181881066274845f64,
38        flattening: 0.003352779454167505,
39    };
40
41    /// [Mars Orbiter Laser Altimeter Ellipsoid](https://tharsis.gsfc.nasa.gov/geodesy.html).
42    pub const MOLA: Ellipsoid = Ellipsoid {
43        equatorial_radius: Length::from_metres(3_396_200f64),
44        polar_radius: Length::from_metres(3_376_198.822143698f64),
45        eccentricity: 0.10836918094475001f64,
46        flattening: 0.005889281507656065f64,
47    };
48
49    /// Creates a new ellipsoid from the given equatorial radius (semi-major axis A) and
50    /// inverse (or reciprocal) flattening.
51    pub fn new(equatorial_radius: Length, inverse_flattening: f64) -> Self {
52        let a = equatorial_radius.as_metres();
53        let f = 1.0 / inverse_flattening;
54        let b = a * (1.0 - f);
55        let e = (1.0 - (b * b) / (a * a)).sqrt();
56        Ellipsoid {
57            equatorial_radius,
58            polar_radius: Length::from_metres(b),
59            eccentricity: e,
60            flattening: f,
61        }
62    }
63
64    /// Returns the equatorial radius (or semi-major axis A) of this ellipsoid.
65    #[inline]
66    pub fn equatorial_radius(&self) -> Length {
67        self.equatorial_radius
68    }
69
70    /// Returns the polar radius (or semi-minor axis B) of this ellipsoid.
71    #[inline]
72    pub fn polar_radius(&self) -> Length {
73        self.polar_radius
74    }
75
76    /// Returns the eccentricity of this ellipsoid.
77    #[inline]
78    pub fn eccentricity(&self) -> f64 {
79        self.eccentricity
80    }
81
82    /// Returns the flattening of this ellipsoid.
83    #[inline]
84    pub fn flattening(&self) -> f64 {
85        self.flattening
86    }
87
88    /// Returns the geocentric radius at the given geodetic latitude: the distance from the Earth's center
89    /// to a position on the spheroid surface at geodetic latitude.
90    ///
91    /// See: [Location-dependent radii](https://en.wikipedia.org/wiki/Earth_radius#Location-dependent_radii)
92    ///
93    /// # Examples
94    ///
95    /// ```
96    /// use jord::{Angle, Length};
97    /// use jord::ellipsoidal::Ellipsoid;
98    ///
99    /// assert_eq!(Ellipsoid::WGS84.equatorial_radius(), Ellipsoid::WGS84.geocentric_radius(Angle::ZERO));
100    /// assert_eq!(Ellipsoid::WGS84.polar_radius(), Ellipsoid::WGS84.geocentric_radius(Angle::from_degrees(90.0)));
101    /// assert_eq!(
102    ///     Length::from_metres(6_367_490.0),
103    ///     Ellipsoid::WGS84.geocentric_radius(Angle::from_degrees(45.0)).round_m()
104    /// );
105    /// ```
106    pub fn geocentric_radius(&self, latitude: Angle) -> Length {
107        let cos_lat = latitude.as_radians().cos();
108        let sin_lat = latitude.as_radians().sin();
109        let a = self.equatorial_radius.as_metres();
110        let b = self.polar_radius.as_metres();
111        let f1 = a * a * cos_lat;
112        let f2 = b * b * sin_lat;
113        let f3 = a * cos_lat;
114        let f4 = b * sin_lat;
115        let r = (((f1 * f1) + (f2 * f2)) / ((f3 * f3) + (f4 * f4))).sqrt();
116        Length::from_metres(r)
117    }
118
119    /// Returns the radius of the parallel of the given geodetic latitude.
120    ///
121    /// See: [Radius of the Earth](https://www.oc.nps.edu/oc2902w/geodesy/radiigeo.pdf)
122    ///
123    /// # Examples
124    ///
125    /// ```
126    /// use jord::{Angle, Length};
127    /// use jord::ellipsoidal::Ellipsoid;
128    ///
129    /// assert_eq!(Ellipsoid::WGS84.equatorial_radius(), Ellipsoid::WGS84.latitude_radius(Angle::ZERO));
130    /// assert_eq!(Length::ZERO, Ellipsoid::WGS84.latitude_radius(Angle::from_degrees(90.0)));
131    /// assert_eq!(Length::ZERO, Ellipsoid::WGS84.latitude_radius(Angle::from_degrees(-90.0)));
132    /// ```
133    pub fn latitude_radius(&self, latitude: Angle) -> Length {
134        if latitude == Angle::QUARTER_CIRCLE || latitude == Angle::NEG_QUARTER_CIRCLE {
135            Length::ZERO
136        } else {
137            self.prime_vertical_radius(latitude) * latitude.as_radians().cos()
138        }
139    }
140
141    /// Returns the radius of curvature of the ellipsoid perpendicular to the meridian at the given geodetic latitude (often denoted `N`).
142    ///
143    /// See: [Radius of the Earth](https://www.oc.nps.edu/oc2902w/geodesy/radiigeo.pdf)
144    ///
145    /// # Examples
146    ///
147    /// ```
148    /// use jord::{Angle, Length};
149    /// use jord::ellipsoidal::Ellipsoid;
150    ///
151    /// assert_eq!(Ellipsoid::WGS84.equatorial_radius(), Ellipsoid::WGS84.prime_vertical_radius(Angle::ZERO));
152    /// assert_eq!(
153    ///     Length::from_metres(6388838.29),
154    ///     Ellipsoid::WGS84.prime_vertical_radius(Angle::from_degrees(45.0)).round_mm()
155    /// );
156    /// ```
157    pub fn prime_vertical_radius(&self, latitude: Angle) -> Length {
158        let e2: f64 = self.eccentricity * self.eccentricity;
159        let sin_lat = latitude.as_radians().sin();
160        let sin_lat2 = sin_lat * sin_lat;
161        let r = self.equatorial_radius.as_metres() / (1.0 - e2 * sin_lat2).sqrt();
162        Length::from_metres(r)
163    }
164
165    /// Radius of curvature in the meridian at the given geodetic latitude (often denoted `M`).
166    ///
167    /// See: [Radius of the Earth](https://www.oc.nps.edu/oc2902w/geodesy/radiigeo.pdf)
168    pub fn meridian_radius(&self, latitude: Angle) -> Length {
169        let e2: f64 = self.eccentricity * self.eccentricity;
170        let sin_lat = latitude.as_radians().sin();
171        let sin_lat2 = sin_lat * sin_lat;
172        let r =
173            self.equatorial_radius.as_metres() * (1.0 - e2) / (1.0 - e2 * sin_lat2).powf(3.0 / 2.0);
174        Length::from_metres(r)
175    }
176
177    /// Returns the mean radius (arithmetic mean) of this ellipsoid.
178    ///
179    /// # Examples
180    ///
181    /// ```
182    /// use jord::Length;
183    /// use jord::ellipsoidal::Ellipsoid;
184    ///
185    /// assert_eq!(Length::from_metres(6_371_008.8), Ellipsoid::WGS84.mean_radius().round_dm());
186    /// ```
187    pub fn mean_radius(&self) -> Length {
188        let a = self.equatorial_radius();
189        let b = self.polar_radius();
190        (2.0 * a + b) / 3.0
191    }
192
193    /// Returns the volumetric radius of this ellipsoid: the radius of sphere of same volume.
194    ///
195    /// # Examples
196    ///
197    /// ```
198    /// use jord::Length;
199    /// use jord::spherical::Sphere;
200    /// use jord::ellipsoidal::Ellipsoid;
201    ///
202    /// let r = (Ellipsoid::WGS84.volumetric_radius().as_metres() * 10.0).round() / 10.0;
203    /// assert_eq!(Sphere::EARTH.radius().as_metres(), r);
204    /// ```
205    pub fn volumetric_radius(&self) -> Length {
206        let a = self.equatorial_radius().as_metres();
207        let b = self.polar_radius().as_metres();
208        let r = (a * a * b).cbrt();
209        Length::from_metres(r)
210    }
211}
212
213impl Surface for Ellipsoid {
214    fn geodetic_to_geocentric_position(&self, pos: GeodeticPosition) -> GeocentricPosition {
215        let nv = pos.horizontal_position().as_vec3();
216        let nx = nv.x();
217        let ny = nv.y();
218        let nz = nv.z();
219        let a = self.equatorial_radius.as_metres();
220        let b = self.polar_radius.as_metres();
221        let m = (a * a) / (b * b);
222        let n = b / ((nx * nx * m) + (ny * ny * m) + (nz * nz)).sqrt();
223        let h = pos.height().as_metres();
224        let cx = n * m * nx + h * nx;
225        let cy = n * m * ny + h * ny;
226        let cz = n * nz + h * nz;
227        GeocentricPosition::from_metres(cx, cy, cz)
228    }
229
230    fn geocentric_to_geodetic_position(&self, pos: GeocentricPosition) -> GeodeticPosition {
231        let pv = pos.as_metres();
232        let px = pv.x();
233        let py = pv.y();
234        let pz = pv.z();
235        let e2 = self.eccentricity * self.eccentricity;
236        let e4 = e2 * e2;
237        let a = self.equatorial_radius.as_metres();
238        let a2 = a * a;
239        let p = (px * px + py * py) / a2;
240        let q = ((1.0 - e2) / a2) * (pz * pz);
241        let r = (p + q - e4) / 6.0;
242        let s = (e4 * p * q) / (4.0 * r * r * r);
243        let t = (1.0 + s + (s * (2.0 + s)).sqrt()).powf(1.0 / 3.0);
244        let u = r * (1.0 + t + 1.0 / t);
245        let v = (u * u + q * e4).sqrt();
246        let w = e2 * (u + v - q) / (2.0 * v);
247        let k = (u + v + w * w).sqrt() - w;
248        let d = k * (px * px + py * py).sqrt() / (k + e2);
249        let h = ((k + e2 - 1.0) / k) * (d * d + pz * pz).sqrt();
250
251        let fs = 1.0 / (d * d + pz * pz).sqrt();
252        let fa = k / (k + e2);
253        let nx = fs * fa * px;
254        let ny = fs * fa * py;
255        let nz = fs * pz;
256        GeodeticPosition::new(NVector::new(Vec3::new(nx, ny, nz)), Length::from_metres(h))
257    }
258}
259
260#[cfg(test)]
261mod tests {
262    use crate::{spherical::Sphere, Angle, Length};
263
264    use super::Ellipsoid;
265
266    #[test]
267    fn wgs84() {
268        let wgs84 = Ellipsoid::new(Length::from_metres(6_378_137.0), 298.257223563);
269        assert_eq!(
270            Ellipsoid::WGS84.equatorial_radius(),
271            wgs84.equatorial_radius()
272        );
273        assert_eq!(Ellipsoid::WGS84.polar_radius(), wgs84.polar_radius());
274        assert_eq!(Ellipsoid::WGS84.eccentricity(), wgs84.eccentricity());
275        assert_eq!(Ellipsoid::WGS84.flattening(), wgs84.flattening());
276    }
277
278    #[test]
279    fn grs80() {
280        let grs80 = Ellipsoid::new(Length::from_metres(6_378_137.0), 298.257222101);
281        assert_eq!(
282            Ellipsoid::GRS80.equatorial_radius(),
283            grs80.equatorial_radius()
284        );
285        assert_eq!(Ellipsoid::GRS80.polar_radius(), grs80.polar_radius());
286        assert_eq!(Ellipsoid::GRS80.eccentricity(), grs80.eccentricity());
287        assert_eq!(Ellipsoid::GRS80.flattening(), grs80.flattening());
288    }
289
290    #[test]
291    fn wgs72() {
292        let wgs72 = Ellipsoid::new(Length::from_metres(6_378_135.0), 298.26);
293        assert_eq!(
294            Ellipsoid::WGS72.equatorial_radius(),
295            wgs72.equatorial_radius()
296        );
297        assert_eq!(Ellipsoid::WGS72.polar_radius(), wgs72.polar_radius());
298        assert_eq!(Ellipsoid::WGS72.eccentricity(), wgs72.eccentricity());
299        assert_eq!(Ellipsoid::WGS72.flattening(), wgs72.flattening());
300    }
301
302    #[test]
303    fn mola() {
304        let mola = Ellipsoid::new(Length::from_metres(3_396_200.0), 169.8);
305        assert_eq!(
306            Ellipsoid::MOLA.equatorial_radius(),
307            mola.equatorial_radius()
308        );
309        assert_eq!(Ellipsoid::MOLA.polar_radius(), mola.polar_radius());
310        assert_eq!(Ellipsoid::MOLA.eccentricity(), mola.eccentricity());
311        assert_eq!(Ellipsoid::MOLA.flattening(), mola.flattening());
312    }
313
314    #[test]
315    fn geocentric_radius() {
316        assert_eq!(
317            Ellipsoid::WGS84.equatorial_radius(),
318            Ellipsoid::WGS84.geocentric_radius(Angle::ZERO)
319        );
320        assert_eq!(
321            Ellipsoid::WGS84.polar_radius(),
322            Ellipsoid::WGS84.geocentric_radius(Angle::from_degrees(90.0))
323        );
324        assert_eq!(
325            Length::from_metres(6_367_490.0),
326            Ellipsoid::WGS84
327                .geocentric_radius(Angle::from_degrees(45.0))
328                .round_m()
329        );
330    }
331
332    #[test]
333    fn latitude_radius() {
334        assert_eq!(
335            Ellipsoid::WGS84.equatorial_radius(),
336            Ellipsoid::WGS84.latitude_radius(Angle::ZERO)
337        );
338        assert_eq!(
339            Length::ZERO,
340            Ellipsoid::WGS84.latitude_radius(Angle::from_degrees(90.0))
341        );
342        assert_eq!(
343            Length::ZERO,
344            Ellipsoid::WGS84.latitude_radius(Angle::from_degrees(-90.0))
345        );
346    }
347
348    #[test]
349    fn prime_vertical_radius() {
350        assert_eq!(
351            Ellipsoid::WGS84.equatorial_radius(),
352            Ellipsoid::WGS84.prime_vertical_radius(Angle::ZERO)
353        );
354        assert_eq!(
355            Length::from_metres(6388838.29),
356            Ellipsoid::WGS84
357                .prime_vertical_radius(Angle::from_degrees(45.0))
358                .round_mm()
359        );
360    }
361
362    #[test]
363    fn mean_radius() {
364        assert_eq!(
365            Length::from_metres(6_371_008.8),
366            Ellipsoid::WGS84.mean_radius().round_dm()
367        );
368    }
369
370    #[test]
371    fn volumetric_radius() {
372        let r = (Ellipsoid::WGS84.volumetric_radius().as_metres() * 10.0).round() / 10.0;
373        assert_eq!(Sphere::EARTH.radius().as_metres(), r);
374    }
375}
jord/ellipsoidal/ellipsoid.rs

jord/ellipsoidal/
ellipsoid.rs
