oxiproj-transformations 0.1.2

//! `geoc` conversion (geographic <-> geocentric latitude).
//!
//! Ported from PROJ src/conversions/geoc.cpp.
//!
//! Converts the geographic (geodetic) latitude into the geocentric latitude
//! (forward) and back (inverse). Only the latitude component (`phi`) is
//! affected; longitude (`lam`), height (`z`) and time (`t`) pass through
//! unchanged. When the latitude is at (or extremely close to) a pole, or when
//! the ellipsoid is a sphere (`es == 0`), the latitude is left untouched.

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Operation, ProjResult, M_HALFPI};

/// The `geoc` conversion operation. Ported from PROJ src/conversions/geoc.cpp.
///
/// Implements `pj_geocentric_latitude`: forward maps geographic latitude to
/// geocentric latitude, inverse maps geocentric latitude back to geographic.
#[derive(Debug)]
struct Geoc {
    /// `1 - es` (the forward scaling factor for `tan(phi)`).
    one_es: f64,
    /// `1 / (1 - es)` (the inverse scaling factor for `tan(phi)`).
    rone_es: f64,
    /// First eccentricity squared; a sphere has `es == 0`.
    es: f64,
}

impl Operation for Geoc {
    // ported verbatim from PROJ geoc.cpp
    #[allow(clippy::float_cmp)]
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let limit = M_HALFPI - 1e-9;
        if v[1] > limit || v[1] < -limit || self.es == 0.0 {
            return Ok(c);
        }
        Ok(Coord::new(
            v[0],
            (self.one_es * v[1].tan()).atan(),
            v[2],
            v[3],
        ))
    }

    // ported verbatim from PROJ geoc.cpp
    #[allow(clippy::float_cmp)]
    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let limit = M_HALFPI - 1e-9;
        if v[1] > limit || v[1] < -limit || self.es == 0.0 {
            return Ok(c);
        }
        Ok(Coord::new(
            v[0],
            (self.rone_es * v[1].tan()).atan(),
            v[2],
            v[3],
        ))
    }

    fn has_inverse(&self) -> bool {
        true
    }
}

/// Constructor. Ported from src/conversions/geoc.cpp.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let e = p.ellipsoid;
    Ok(TransBuild::new(
        Box::new(Geoc {
            one_es: e.one_es,
            rone_es: e.rone_es,
            es: e.es,
        }),
        IoUnits::Radians,
        IoUnits::Radians,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::TransParamLookup;
    use oxiproj_core::{Ellipsoid, DEG_TO_RAD};

    struct NoParams;
    impl TransParamLookup for NoParams {
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_f64(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_str(&self, _key: &str) -> Option<&str> {
            None
        }
        fn get_bool(&self, _key: &str) -> bool {
            false
        }
        fn exists(&self, _key: &str) -> bool {
            false
        }
    }

    #[test]
    fn geoc_forward_inverse_wgs84() {
        let ell = Ellipsoid::named("WGS84").unwrap();
        let b = new(&TransParams {
            ellipsoid: &ell,
            params: &NoParams,
            registry: None,
        })
        .unwrap();
        let op = b.operation;

        let out = op
            .forward_4d(Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0))
            .unwrap();
        // Geocentric latitude for geographic latitude 55 deg on WGS84.
        assert!((out.v()[1] - 54.818973309 * DEG_TO_RAD).abs() < 1e-9);
        // Longitude passes through unchanged.
        assert!((out.v()[0] - 12.0 * DEG_TO_RAD).abs() < 1e-15);

        let back = op.inverse_4d(out).unwrap();
        // Round-trips back to the geographic latitude.
        assert!((back.v()[1] - 55.0 * DEG_TO_RAD).abs() < 1e-9);
    }
}