mgrs-convert 0.1.0

//! Transverse Mercator, Krueger series to order n^6.
//!
//! The series coefficients are transcribed from equations 35 and 36 of Karney,
//! "Transverse Mercator with an accuracy of a few nanometers" (J. Geodesy 85(8),
//! 2011, arXiv:1002.1417), the same order GeographicLib's `TransverseMercator`
//! uses. The differential suite against `TMcoords.dat` and `GeoConvert` is what
//! proves the transcription. Truncation error inside the UTM extent is below
//! 5 nm.
//
// Bounded floating-point arithmetic only; see the note in `math.rs`.

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

/// Central scale factor on the UTM central meridian.
pub(crate) const K0: f64 = 0.9996;

/// Third flattening, n = f / (2 - f).
const N: f64 = F / (2.0 - F);
const N2: f64 = N * N;
const N3: f64 = N2 * N;
const N4: f64 = N2 * N2;
const N5: f64 = N4 * N;
const N6: f64 = N4 * N2;

/// Rectifying radius, A1 = a/(1+n) (1 + n^2/4 + n^4/64 + n^6/256).
pub(crate) const A1: f64 = A / (1.0 + N) * (1.0 + N2 / 4.0 + N4 / 64.0 + N6 / 256.0);

/// Forward series coefficients, Karney 2011 equation 35.
const ALPHA: [f64; 6] = [
    N / 2.0 - 2.0 / 3.0 * N2 + 5.0 / 16.0 * N3 + 41.0 / 180.0 * N4 - 127.0 / 288.0 * N5
        + 7891.0 / 37800.0 * N6,
    13.0 / 48.0 * N2 - 3.0 / 5.0 * N3 + 557.0 / 1440.0 * N4 + 281.0 / 630.0 * N5
        - 1_983_433.0 / 1_935_360.0 * N6,
    61.0 / 240.0 * N3 - 103.0 / 140.0 * N4 + 15061.0 / 26880.0 * N5 + 167_603.0 / 181_440.0 * N6,
    49561.0 / 161_280.0 * N4 - 179.0 / 168.0 * N5 + 6_601_661.0 / 7_257_600.0 * N6,
    34729.0 / 80640.0 * N5 - 3_418_889.0 / 1_995_840.0 * N6,
    212_378_941.0 / 319_334_400.0 * N6,
];

/// Reverse series coefficients, Karney 2011 equation 36.
const BETA: [f64; 6] = [
    N / 2.0 - 2.0 / 3.0 * N2 + 37.0 / 96.0 * N3 - 1.0 / 360.0 * N4 - 81.0 / 512.0 * N5
        + 96199.0 / 604_800.0 * N6,
    N2 / 48.0 + N3 / 15.0 - 437.0 / 1440.0 * N4 + 46.0 / 105.0 * N5
        - 1_118_711.0 / 3_870_720.0 * N6,
    17.0 / 480.0 * N3 - 37.0 / 840.0 * N4 - 209.0 / 4480.0 * N5 + 5569.0 / 90720.0 * N6,
    4397.0 / 161_280.0 * N4 - 11.0 / 504.0 * N5 - 830_251.0 / 7_257_600.0 * N6,
    4583.0 / 161_280.0 * N5 - 108_847.0 / 3_991_680.0 * N6,
    20_648_693.0 / 638_668_800.0 * N6,
];

/// Clenshaw summation of S = sum c_j sin(2j zeta) for complex zeta = `xi + i*eta`.
///
/// Returns (Re S, Im S) = (sum c_j sin 2j xi cosh 2j eta, sum c_j cos 2j xi sinh 2j eta).
fn clenshaw(c: &[f64; 6], xi: f64, eta: f64) -> (f64, f64) {
    // theta = 2 zeta; X = 2 cos theta, S = b1 sin theta with the standard
    // sine-series recurrence b_j = c_j + X b_{j+1} - b_{j+2}.
    let (s2x, c2x) = (2.0 * xi).sin_cos();
    let ch2e = (2.0 * eta).cosh();
    let sh2e = (2.0 * eta).sinh();
    let xr = 2.0 * c2x * ch2e;
    let xi_im = -2.0 * s2x * sh2e;
    let (mut b1r, mut b1i) = (0.0_f64, 0.0_f64);
    let (mut b2r, mut b2i) = (0.0_f64, 0.0_f64);
    for &cj in c.iter().rev() {
        let tr = cj + (xr * b1r - xi_im * b1i) - b2r;
        let ti = (xr * b1i + xi_im * b1r) - b2i;
        b2r = b1r;
        b2i = b1i;
        b1r = tr;
        b1i = ti;
    }
    let sr = s2x * ch2e;
    let si = c2x * sh2e;
    (b1r * sr - b1i * si, b1r * si + b1i * sr)
}

/// Forward projection. Input is latitude and the longitude offset from the
/// central meridian, both in degrees. Output is (x, y) in metres: x east of the
/// central meridian, y north of the equator, both scaled by k0, without false
/// easting or northing.
pub(crate) fn forward(lat: f64, dlon: f64) -> (f64, f64) {
    let phi = lat.to_radians();
    let lam = dlon.to_radians();
    let tau = phi.tan();
    let taup = taupf(tau);
    let (slam, clam) = lam.sin_cos();
    let xip = taup.atan2(clam);
    let etap = (slam / taup.hypot(clam)).asinh();
    let (dre, dim) = clenshaw(&ALPHA, xip, etap);
    let xi = xip + dre;
    let eta = etap + dim;
    (K0 * A1 * eta, K0 * A1 * xi)
}

/// Reverse projection. Input is (x, y) as produced by [`forward`]. Output is
/// latitude and the longitude offset from the central meridian, in degrees.
pub(crate) fn reverse(x: f64, y: f64) -> (f64, f64) {
    let xi = y / (K0 * A1);
    let eta = x / (K0 * A1);
    let (dre, dim) = clenshaw(&BETA, xi, eta);
    let xip = xi - dre;
    let etap = eta - dim;
    let (sxip, cxip) = xip.sin_cos();
    let sheta = etap.sinh();
    let taup = sxip / sheta.hypot(cxip);
    let tau = tauf(taup);
    let lat = tau.atan().to_degrees();
    let dlon = sheta.atan2(cxip).to_degrees();
    (lat, dlon)
}

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

    #[test]
    fn constants_pinned() {
        // Invariant 9. The rectifying radius for WGS84, a value every reference
        // implementation agrees on to the nanometre.
        assert!((A1 - 6_367_449.145_823_415).abs() < 1e-6);
        assert_eq!(K0, 0.9996);
    }

    #[test]
    fn equator_on_central_meridian_is_exact_zero() {
        let (x, y) = forward(0.0, 0.0);
        assert_eq!(x, 0.0);
        assert_eq!(y, 0.0);
    }

    #[test]
    fn forward_reverse_round_trip_is_tight() {
        // Invariant 7, spot check: our own round trip must close to 0.1 um.
        for &(lat, dlon) in &[
            (40.748333, -1.985278),
            (-33.856784, 2.215297),
            (83.5, 2.9),
            (-79.9, -2.9),
            (0.5, -0.5),
        ] {
            let (x, y) = forward(lat, dlon);
            let (lat2, dlon2) = reverse(x, y);
            let (x2, y2) = forward(lat2, dlon2);
            assert!((x - x2).abs() < 1e-7, "x drift {lat} {dlon}");
            assert!((y - y2).abs() < 1e-7, "y drift {lat} {dlon}");
        }
    }
}