mgrs-convert 0.1.0

//! Universal Polar Stereographic, closed-form conformal mapping.
//!
//! The formulation is the polar stereographic projection from Karney 2011
//! (equivalently Snyder 1987, equations 21-33 and 21-34): the projected radius
//! is rho = 2 a k0 t' / C where t' = exp(-psi) is the conformal colatitude
//! tangent and C = (1-f) exp(e atanh e) normalizes the scale to k0 at the pole.
//
// Bounded floating-point arithmetic only; see the note in `math.rs`.

use crate::math::{es, tauf, taupf, A, F};

/// Central scale factor at the pole.
pub(crate) const K0: f64 = 0.994;

/// False easting and false northing, both 2,000,000 m.
pub(crate) const FALSE_EN: f64 = 2_000_000.0;

/// C = (1-f) exp(e atanh e), the pole-scale normalization constant.
fn c_ps() -> f64 {
    (1.0 - F) * (es() * es().atanh()).exp()
}

/// Forward projection for one polar cap. `north` picks the cap, `lat` is the
/// true latitude in degrees (negative in the south), `lon` in degrees. Returns
/// (easting, northing) with false easting and northing applied.
pub(crate) fn forward(north: bool, lat: f64, lon: f64) -> (f64, f64) {
    // Work in the north-polar aspect: flip the south cap.
    let lat_p = if north { lat } else { -lat };
    // Latitude exactly 90 must project to the pole exactly so the grid letters
    // and digits at the pole are exact; tan then loses that exactness.
    let rho = if lat_p == 90.0 {
        0.0
    } else {
        let tau = lat_p.to_radians().tan();
        let taup = taupf(tau);
        // t' = exp(-psi) = 1 / (tau' + sqrt(1 + tau'^2)).
        let tprime = 1.0 / (taup + (1.0 + taup * taup).sqrt());
        2.0 * A * K0 * tprime / c_ps()
    };
    let (slam, clam) = lon.to_radians().sin_cos();
    let x = FALSE_EN + rho * slam;
    let y = if north {
        FALSE_EN - rho * clam
    } else {
        FALSE_EN + rho * clam
    };
    (x, y)
}

/// Reverse projection for one polar cap. Returns (lat, lon) in degrees, the
/// latitude signed for the hemisphere. The pole itself reverses to longitude 0.
pub(crate) fn reverse(north: bool, easting: f64, northing: f64) -> (f64, f64) {
    let x = easting - FALSE_EN;
    let yp = if north {
        FALSE_EN - northing
    } else {
        northing - FALSE_EN
    };
    let rho = x.hypot(yp);
    let lat_p = if rho == 0.0 {
        90.0
    } else {
        let tprime = rho * c_ps() / (2.0 * A * K0);
        // tau' = (1/t' - t') / 2.
        let taup = (1.0 / tprime - tprime) / 2.0;
        tauf(taup).atan().to_degrees()
    };
    let lon = x.atan2(yp).to_degrees();
    let lat = if north { lat_p } else { -lat_p };
    (lat, lon)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn constants_pinned() {
        // Invariant 9.
        assert_eq!(K0, 0.994);
        assert_eq!(FALSE_EN, 2_000_000.0);
    }

    #[test]
    fn poles_project_to_the_exact_grid_center() {
        assert_eq!(forward(true, 90.0, 0.0), (2_000_000.0, 2_000_000.0));
        assert_eq!(forward(false, -90.0, 77.0), (2_000_000.0, 2_000_000.0));
        let (lat, lon) = reverse(true, 2_000_000.0, 2_000_000.0);
        assert_eq!((lat, lon), (90.0, 0.0));
    }

    #[test]
    fn round_trip_is_tight() {
        for &(north, lat, lon) in &[
            (true, 87.3, -122.0),
            (true, 84.001, 179.5),
            (false, -85.75, 31.0),
            (false, -80.001, -0.25),
        ] {
            let (x, y) = forward(north, lat, lon);
            let (lat2, lon2) = reverse(north, x, y);
            let (x2, y2) = forward(north, lat2, lon2);
            assert!((x - x2).abs() < 1e-7, "x drift at {lat} {lon}");
            assert!((y - y2).abs() < 1e-7, "y drift at {lat} {lon}");
        }
    }
}