Trait based refactor of Ellipsoid

The Ellipsoid impl was huge and unmaintainable.

Now it has been split into an EllipsoidBase trait,
defining the fundamental shape and size parameters,
and a number of more specialized traits for latitudes,
meridian-geometry, geodesics, cartesians (geocart),
and gravity.

Also lib.rs was cleaned up to reflect the new structure,
and the prelude has been rebuilt using thematical modules:
coord, ctx, and ellps, representing coordinate, context,
and ellipsoidal material, respectively.

Likewise, the "authoring" extended prelude has been
restructured through the introduction of the grd,
ops, and parse modules.

Finally, a number of documentation improvements were
introduced while refactoring the existing material to
reflect the new structure.

For usage based on the "use geodesy::prelude::*" idiom,
no user visible changes are expected.

Author: busstoptaktik

CHANGELOG.md

@@ -9,14 +9,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
-- CoordinateSet: xyz(), set_xyz(), xyzt(), set_xyzt() methods
+- `CoordinateSet`: xyz(), set_xyz(), xyzt(), set_xyzt() methods
 - Vector space operators (Add, Sub, Mul, Div) for all built
   in coordinate tuple types (Coor4D, Coor3D, Coor2D, Coor32)
+- `TriaxialEllpisoid`, mostly as a placeholder
 
 ### Changed
 
 - `CoordinateTuple` trait now requires implementation of the constructor
   method `new(fill: f64)`, returning an object of `dim()` copies of `fill`.
+- The huge `Ellipsoid`-implementation switched to a new trait `EllipsoidBase`,
+  and a number of associated traits requiring `EllipsoidBase`.
+- Meaning of `EllipsoidBase::aspect_ratio()` switched from *b / a* to the
+  apparently more common *a / b*
+- Major restructuring and clean up of `lib.rs`. Only marignally visible externally,
+  if using `use geodesy::prelude::*`
 
 ### Removed
 

examples/03-user_defined_operators.rs

@@ -51,7 +51,7 @@ pub fn add42_constructor(parameters: &RawParameters, ctx: &dyn Context) -> Resul
 }
 
 fn main() -> anyhow::Result<()> {
-    let mut prv = geodesy::Minimal::new();
+    let mut prv = geodesy::prelude::Minimal::new();
     prv.register_op("add42", OpConstructor(add42_constructor));
     let add42 = prv.op("add42")?;
 

examples/05-pq.rs

@@ -66,7 +66,7 @@ fn main() -> Result<(), anyhow::Error> {
 
     // We use ::new() instead of ::default() in order to gain access to the
     // BUILTIN_ADAPTORS
-    let mut ctx = geodesy::Minimal::new();
+    let mut ctx = geodesy::prelude::Minimal::new();
     trace!("trace message 2");
     debug!("debug message 2");
 

examples/06-user_defined_coordinate_types_and_containers.rs

@@ -87,7 +87,7 @@ fn main() -> Result<(), anyhow::Error> {
 
     // We use ::new() instead of ::default() in order to gain access to the
     // BUILTIN_ADAPTORS (geo:in, geo:out etc.)
-    let mut ctx = geodesy::Minimal::new();
+    let mut ctx = geodesy::prelude::Minimal::new();
     trace!("have context");
 
     let copenhagen = Coor4D([55., 12., 0., 0.]);

examples/08-user_defined_operators_using_proj.rs

@@ -171,7 +171,7 @@ pub fn proj_constructor(parameters: &RawParameters, _ctx: &dyn Context) -> Resul
 }
 
 fn main() -> anyhow::Result<()> {
-    let mut prv = geodesy::Minimal::new();
+    let mut prv = geodesy::prelude::Minimal::new();
     prv.register_op("proj", OpConstructor(proj_constructor));
     let e = Ellipsoid::default();
 

src/coordinate/coor32.rs

@@ -1,7 +1,7 @@
 use super::*;
 
 /// Generic 2D Coordinate tuple, with no fixed interpretation of the elements.
-/// A tiny coordinate type: Just one fourth the weight of a [`Coor4D`](crate::Coor4D).
+/// A tiny coordinate type: Just one fourth the weight of a [`Coor4D`](super::Coor4D).
 /// Probably only useful for small scale world maps, without too much zoom.
 #[derive(Debug, Default, PartialEq, Copy, Clone)]
 pub struct Coor32(pub [f32; 2]);

src/coordinate/mod.rs

@@ -63,7 +63,7 @@ where
 }
 
 /// For Rust Geodesy, the ISO-19111 concept of `DirectPosition` is represented
-/// as a `geodesy::Coo4D`.
+/// as a `geodesy::Coor4D`.
 ///
 /// The strict connection between the ISO19107 "DirectPosition" datatype
 /// and the ISO19111/OGC Topic 2 "CoordinateSet" interface (i.e. trait)

src/coordinate/tuple.rs

@@ -109,7 +109,7 @@ all_coord_operators!(Coor32, Coor32, coor32);
 
 /// CoordinateTuple is the ISO-19111 atomic spatial/spatiotemporal
 /// referencing element. So loosely speaking, a CoordinateSet is a
-///  collection of CoordinateTuples.
+/// collection of CoordinateTuples.
 ///
 /// Note that (despite the formal name) the underlying data structure
 /// need not be a tuple: It can be any item, for which it makes sense
@@ -357,7 +357,7 @@ pub trait CoordinateTuple {
     /// # See also:
     ///
     /// [`hypot3`](Self::hypot3),
-    /// [`distance`](crate::ellipsoid::Ellipsoid::distance)
+    /// [`distance`](crate::ellps::Geodesics::distance)
     ///
     /// # Examples
     ///
@@ -391,7 +391,7 @@ pub trait CoordinateTuple {
     /// # See also:
     ///
     /// [`hypot2()`](Self::hypot2),
-    /// [`distance`](crate::ellipsoid::Ellipsoid::distance)
+    /// [`distance`](crate::ellps::Geodesics::distance)
     ///
     /// # Examples
     ///

src/ellipsoid/biaxial.rs

@@ -0,0 +1,99 @@
+use crate::prelude::*;
+
+/// An ellipsoid of revolution.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Ellipsoid {
+    a: f64,
+    f: f64,
+}
+
+/// GRS80 is the default ellipsoid.
+impl Default for Ellipsoid {
+    fn default() -> Ellipsoid {
+        Ellipsoid::new(6_378_137.0, 1. / 298.257_222_100_882_7)
+    }
+}
+
+impl EllipsoidBase for Ellipsoid {
+    fn semimajor_axis(&self) -> f64 {
+        self.a
+    }
+
+    fn flattening(&self) -> f64 {
+        self.f
+    }
+}
+
+/// Constructors for `Ellipsoid`
+impl Ellipsoid {
+    /// User defined ellipsoid
+    #[must_use]
+    pub fn new(semimajor_axis: f64, flattening: f64) -> Ellipsoid {
+        Ellipsoid {
+            a: semimajor_axis,
+            f: flattening,
+        }
+    }
+
+    /// Predefined ellipsoid; built-in, defined in asset collections, or given as a
+    /// string formatted (a, rf) tuple, e.g. "6378137, 298.25"
+    pub fn named(name: &str) -> Result<Ellipsoid, Error> {
+        // Is it one of the few builtins?
+        if let Some(index) = super::constants::ELLIPSOID_LIST
+            .iter()
+            .position(|&ellps| ellps.0 == name)
+        {
+            let e = super::constants::ELLIPSOID_LIST[index];
+            let ax: f64 = e.1.parse().unwrap();
+            let rf: f64 = e.3.parse().unwrap();
+            // EPSG convention: zero reciproque flattening indicates zero flattening
+            let f = if rf != 0.0 { 1.0 / rf } else { rf };
+            return Ok(Ellipsoid::new(ax, f));
+        }
+
+        // The "semimajor, reciproque-flattening" form, e.g. "6378137, 298.3"
+        let a_and_rf = name.split(',').collect::<Vec<_>>();
+        if a_and_rf.len() == 2_usize {
+            if let Ok(a) = a_and_rf[0].trim().parse::<f64>() {
+                if let Ok(rf) = a_and_rf[1].trim().parse::<f64>() {
+                    return Ok(Ellipsoid::new(a, 1. / rf));
+                }
+            }
+        }
+
+        // TODO: Search asset collection
+        Err(Error::NotFound(
+            String::from(name),
+            String::from("Ellipsoid::named()"),
+        ))
+    }
+}
+
+// ----- Tests ---------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use crate::ellps::Ellipsoid;
+    use crate::ellps::EllipsoidBase;
+    use crate::ellps::Meridians;
+    use crate::Error;
+
+    #[test]
+    fn test_ellipsoid() -> Result<(), Error> {
+        // Constructors
+        let ellps = Ellipsoid::named("intl")?;
+        assert_eq!(ellps.flattening(), 1. / 297.);
+
+        let ellps = Ellipsoid::named("6378137, 298.25")?;
+        assert_eq!(ellps.semimajor_axis(), 6378137.0);
+        assert_eq!(ellps.flattening(), 1. / 298.25);
+
+        let ellps = Ellipsoid::named("GRS80")?;
+        assert_eq!(ellps.semimajor_axis(), 6378137.0);
+        assert_eq!(ellps.flattening(), 1. / 298.257_222_100_882_7);
+
+        assert!((ellps.normalized_meridian_arc_unit() - 0.998_324_298_423_041_5).abs() < 1e-13);
+        assert!((4.0 * ellps.meridian_quadrant() - 40_007_862.916_921_8).abs() < 1e-7);
+        Ok(())
+    }
+}

src/ellipsoid/cartesians.rs src/ellipsoid/geocart.rs

@@ -2,9 +2,8 @@ use super::*;
 
 use std::f64::consts::FRAC_PI_2;
 
-impl Ellipsoid {
-    // ----- Cartesian <--> Geographic conversion ----------------------------------
-
+/// Geographic <--> Cartesian conversion
+pub trait GeoCart: EllipsoidBase {
     /// Geographic to cartesian conversion.
     ///
     /// Follows the the derivation given by
@@ -13,11 +12,8 @@ impl Ellipsoid {
     #[must_use]
     #[allow(non_snake_case)] // make it possible to mimic math notation from original paper
     #[allow(clippy::many_single_char_names)] // ditto
-    pub fn cartesian(&self, geographic: &Coor4D) -> Coor4D {
-        let lam = geographic[0];
-        let phi = geographic[1];
-        let h = geographic[2];
-        let t = geographic[3];
+    fn cartesian<C: CoordinateTuple>(&self, geographic: &C) -> Coor4D {
+        let (lam, phi, h, t) = geographic.xyzt();
 
         let N = self.prime_vertical_radius_of_curvature(phi);
         let cosphi = phi.cos();
@@ -32,26 +28,23 @@ impl Ellipsoid {
         Coor4D::raw(X, Y, Z, t)
     }
 
-    /// Cartesian to geogaphic conversion.
+    /// Cartesian to geographic conversion.
     ///
     /// Follows the the derivation given by
     /// Bowring ([1976](crate::Bibliography::Bow76) and
     /// [1985](crate::Bibliography::Bow85))
     #[must_use]
     #[allow(non_snake_case)] // make it possible to mimic math notation from original paper
     #[allow(clippy::many_single_char_names)] // ditto
-    pub fn geographic(&self, cartesian: &Coor4D) -> Coor4D {
-        let X = cartesian[0];
-        let Y = cartesian[1];
-        let Z = cartesian[2];
-        let t = cartesian[3];
+    fn geographic<C: CoordinateTuple>(&self, cartesian: &C) -> Coor4D {
+        let (X, Y, Z, t) = cartesian.xyzt();
 
         // We need a few additional ellipsoidal parameters
         let b = self.semiminor_axis();
         let eps = self.second_eccentricity_squared();
         let es = self.eccentricity_squared();
 
-        // The longitude is straightforward
+        // The longitude is straightforward: Plain geometry in the equatoreal plane
         let lam = Y.atan2(X);
 
         // The perpendicular distance from the point coordinate to the Z-axis
@@ -79,26 +72,29 @@ impl Ellipsoid {
         // Fukushima (1999), Appendix B: Equivalent to Even Rouault's, implementation,
         // but not as clear - although a bit faster due to the substitution of sqrt
         // for hypot.
-        let T = (Z * self.a) / (p * b);
+        let a = self.semimajor_axis();
+        let T = (Z * a) / (p * b);
         let c = 1.0 / (1.0 + T * T).sqrt();
         let s = c * T;
 
         let phi_num = Z + eps * b * s.powi(3);
-        let phi_denom = p - es * self.a * c.powi(3);
+        let phi_denom = p - es * a * c.powi(3);
         let phi = phi_num.atan2(phi_denom);
 
         let lenphi = phi_num.hypot(phi_denom);
         let sinphi = phi_num / lenphi;
         let cosphi = phi_denom / lenphi;
 
+        let a = self.semimajor_axis();
+
         // We already have sinphi and es, so we can compute the radius
         // of curvature faster by inlining, rather than calling the
         // prime_vertical_radius_of_curvature() method.
-        let N = self.a / (1.0 - sinphi.powi(2) * es).sqrt();
+        let N = a / (1.0 - sinphi.powi(2) * es).sqrt();
 
         // Bowring (1985), as quoted by Burtch (2006), suggests this expression
         // as more accurate than the commonly used h = p / cosphi - N;
-        let h = p * cosphi + Z * sinphi - self.a * self.a / N;
+        let h = p * cosphi + Z * sinphi - a * a / N;
 
         Coor4D::raw(lam, phi, h, t)
     }