oxiproj-core 0.1.2

//! Karney auxiliary-latitude series ported from PROJ `src/latitudes.cpp`.
//!
//! This module provides the uniform interface for converting between any pair
//! of auxiliary latitudes using series expansions in the third flattening `n`,
//! adapted from:
//!
//!   C. F. F. Karney, *On auxiliary latitudes*, Survey Review 56, 165-180
//!   (2024). <https://doi.org/10.1080/00396265.2023.2217604>
//!
//! The coefficient table (`COEFFS`) and pointer table (`PTRS`) are
//! machine-generated by the Maxima code `auxlat.mac` bundled with
//! GeographicLib, and are copied verbatim from `src/latitudes.cpp`. Only the
//! upper-triangular portion of the conversion matrices is included, and for
//! conversions among the "geometric" latitudes (geographic, parametric,
//! geocentric, rectifying) the Taylor series are expansions in `n^2`, so the
//! odd (zero) coefficients are excluded.

use crate::error::{ProjError, ProjResult};

/// Auxiliary latitude kinds. Ported from src/latitudes.cpp.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[repr(usize)]
#[non_exhaustive]
pub enum AuxLat {
    /// Geographic latitude, phi.
    Geographic = 0,
    /// Parametric latitude, beta.
    Parametric = 1,
    /// Geocentric latitude, theta.
    Geocentric = 2,
    /// Rectifying latitude, mu.
    Rectifying = 3,
    /// Conformal latitude, chi.
    Conformal = 4,
    /// Authalic latitude, xi.
    Authalic = 5,
}

impl AuxLat {
    /// Number of supported auxiliary latitudes. Ported from
    /// `AuxLat::NUMBER` in src/latitudes.cpp.
    pub const NUMBER: usize = 6;
    /// Order of the Taylor series in `n`. Ported from `AuxLat::ORDER` in
    /// src/latitudes.cpp.
    pub const ORDER: usize = 6;
}

/// Concatenated upper-triangular conversion-matrix coefficients.
///
/// Copied verbatim from the `coeffs[]` array in src/latitudes.cpp (generated
/// by Maxima on 2025-03-23). The starting offset of each individual matrix is
/// given by [`PTRS`].
static COEFFS: &[f64] = &[
    // C[phi,mu]; even coeffs only
    3.0 / 2.0,
    -27.0 / 32.0,
    269.0 / 512.0,
    21.0 / 16.0,
    -55.0 / 32.0,
    6759.0 / 4096.0,
    151.0 / 96.0,
    -417.0 / 128.0,
    1097.0 / 512.0,
    -15543.0 / 2560.0,
    8011.0 / 2560.0,
    293393.0 / 61440.0,
    // C[phi,chi]
    2.0,
    -2.0 / 3.0,
    -2.0,
    116.0 / 45.0,
    26.0 / 45.0,
    -2854.0 / 675.0,
    7.0 / 3.0,
    -8.0 / 5.0,
    -227.0 / 45.0,
    2704.0 / 315.0,
    2323.0 / 945.0,
    56.0 / 15.0,
    -136.0 / 35.0,
    -1262.0 / 105.0,
    73814.0 / 2835.0,
    4279.0 / 630.0,
    -332.0 / 35.0,
    -399572.0 / 14175.0,
    4174.0 / 315.0,
    -144838.0 / 6237.0,
    601676.0 / 22275.0,
    // C[phi,xi]
    4.0 / 3.0,
    4.0 / 45.0,
    -16.0 / 35.0,
    -2582.0 / 14175.0,
    60136.0 / 467775.0,
    28112932.0 / 212837625.0,
    46.0 / 45.0,
    152.0 / 945.0,
    -11966.0 / 14175.0,
    -21016.0 / 51975.0,
    251310128.0 / 638512875.0,
    3044.0 / 2835.0,
    3802.0 / 14175.0,
    -94388.0 / 66825.0,
    -8797648.0 / 10945935.0,
    6059.0 / 4725.0,
    41072.0 / 93555.0,
    -1472637812.0 / 638512875.0,
    768272.0 / 467775.0,
    455935736.0 / 638512875.0,
    4210684958.0 / 1915538625.0,
    // C[mu,phi]; even coeffs only
    -3.0 / 2.0,
    9.0 / 16.0,
    -3.0 / 32.0,
    15.0 / 16.0,
    -15.0 / 32.0,
    135.0 / 2048.0,
    -35.0 / 48.0,
    105.0 / 256.0,
    315.0 / 512.0,
    -189.0 / 512.0,
    -693.0 / 1280.0,
    1001.0 / 2048.0,
    // C[mu,chi]
    1.0 / 2.0,
    -2.0 / 3.0,
    5.0 / 16.0,
    41.0 / 180.0,
    -127.0 / 288.0,
    7891.0 / 37800.0,
    13.0 / 48.0,
    -3.0 / 5.0,
    557.0 / 1440.0,
    281.0 / 630.0,
    -1983433.0 / 1935360.0,
    61.0 / 240.0,
    -103.0 / 140.0,
    15061.0 / 26880.0,
    167603.0 / 181440.0,
    49561.0 / 161280.0,
    -179.0 / 168.0,
    6601661.0 / 7257600.0,
    34729.0 / 80640.0,
    -3418889.0 / 1995840.0,
    212378941.0 / 319334400.0,
    // C[chi,phi]
    -2.0,
    2.0 / 3.0,
    4.0 / 3.0,
    -82.0 / 45.0,
    32.0 / 45.0,
    4642.0 / 4725.0,
    5.0 / 3.0,
    -16.0 / 15.0,
    -13.0 / 9.0,
    904.0 / 315.0,
    -1522.0 / 945.0,
    -26.0 / 15.0,
    34.0 / 21.0,
    8.0 / 5.0,
    -12686.0 / 2835.0,
    1237.0 / 630.0,
    -12.0 / 5.0,
    -24832.0 / 14175.0,
    -734.0 / 315.0,
    109598.0 / 31185.0,
    444337.0 / 155925.0,
    // C[chi,mu]
    -1.0 / 2.0,
    2.0 / 3.0,
    -37.0 / 96.0,
    1.0 / 360.0,
    81.0 / 512.0,
    -96199.0 / 604800.0,
    -1.0 / 48.0,
    -1.0 / 15.0,
    437.0 / 1440.0,
    -46.0 / 105.0,
    1118711.0 / 3870720.0,
    -17.0 / 480.0,
    37.0 / 840.0,
    209.0 / 4480.0,
    -5569.0 / 90720.0,
    -4397.0 / 161280.0,
    11.0 / 504.0,
    830251.0 / 7257600.0,
    -4583.0 / 161280.0,
    108847.0 / 3991680.0,
    -20648693.0 / 638668800.0,
    // C[xi,phi]
    -4.0 / 3.0,
    -4.0 / 45.0,
    88.0 / 315.0,
    538.0 / 4725.0,
    20824.0 / 467775.0,
    -44732.0 / 2837835.0,
    34.0 / 45.0,
    8.0 / 105.0,
    -2482.0 / 14175.0,
    -37192.0 / 467775.0,
    -12467764.0 / 212837625.0,
    -1532.0 / 2835.0,
    -898.0 / 14175.0,
    54968.0 / 467775.0,
    100320856.0 / 1915538625.0,
    6007.0 / 14175.0,
    24496.0 / 467775.0,
    -5884124.0 / 70945875.0,
    -23356.0 / 66825.0,
    -839792.0 / 19348875.0,
    570284222.0 / 1915538625.0,
];

/// Starting offsets into [`COEFFS`] for each conversion matrix.
///
/// Copied verbatim from the `ptrs[]` array in src/latitudes.cpp. The length is
/// `AuxLat::NUMBER^2 + 1 == 37`. Entry `k = 6*auxout + auxin` gives the start
/// offset; `PTRS[k] == PTRS[k + 1]` means the conversion is unsupported.
static PTRS: [usize; 37] = [
    0, 0, 0, 0, 12, 33, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 66, 66, 66, 66, 87, 87,
    108, 108, 108, 129, 129, 129, 150, 150, 150, 150, 150, 150,
];

/// Ported from src/latitudes.cpp (`polyval`). Horner evaluation; `n` is the
/// highest coefficient index in `p` (so the polynomial has `n+1` terms when
/// `n >= 0`). A negative `n` yields `0.0`.
///
/// Using `.get` avoids any panic; in practice indices are in-bounds.
pub fn pj_polyval(x: f64, p: &[f64], n: isize) -> f64 {
    if n < 0 {
        return 0.0;
    }
    let mut idx = n as usize;
    let mut y = match p.get(idx) {
        Some(v) => *v,
        None => return 0.0,
    };
    while idx > 0 {
        idx -= 1;
        let c = match p.get(idx) {
            Some(v) => *v,
            None => 0.0,
        };
        y = y * x + c;
    }
    y
}

/// Ported from src/latitudes.cpp (`clenshaw` summation).
///
/// Evaluates `y = sum(f[k] * sin((2*k+2) * zeta), k, 0, K-1)` by Clenshaw
/// summation, where `zeta` is specified by its sine and cosine (`szeta`,
/// `czeta`).
pub fn pj_clenshaw(szeta: f64, czeta: f64, f: &[f64], k: usize) -> f64 {
    let mut u0 = 0.0_f64;
    let mut u1 = 0.0_f64;
    let x = 2.0 * (czeta - szeta) * (czeta + szeta);
    let mut kk = k;
    while kk > 0 {
        kk -= 1;
        let fk = match f.get(kk) {
            Some(v) => *v,
            None => 0.0,
        };
        let t = x * u0 - u1 + fk;
        u1 = u0;
        u0 = t;
    }
    2.0 * szeta * czeta * u0
}

/// Ported from src/latitudes.cpp. Converts an auxiliary latitude given as an
/// angle `zeta` to the target latitude (scalar form): computes sin/cos and adds
/// the Clenshaw series correction.
pub fn pj_auxlat_convert(zeta: f64, f: &[f64], k: usize) -> f64 {
    pj_auxlat_convert_sc(zeta, zeta.sin(), zeta.cos(), f, k)
}

/// Ported from src/latitudes.cpp. Same as `pj_auxlat_convert` but the caller
/// supplies `sin(zeta)`/`cos(zeta)`.
pub fn pj_auxlat_convert_sc(zeta: f64, szeta: f64, czeta: f64, f: &[f64], k: usize) -> f64 {
    zeta + pj_clenshaw(szeta, czeta, f, k)
}

/// Ported from src/latitudes.cpp (the `void` sin/cos-returning conversion).
/// Returns `(sin(eta), cos(eta))` where `eta = zeta + delta`, `delta` is the
/// Clenshaw series value, and the sin/cos are combined via angle-addition.
/// This provides high relative accuracy near the poles.
pub fn pj_auxlat_convert_seta_ceta(szeta: f64, czeta: f64, f: &[f64], k: usize) -> (f64, f64) {
    let delta = pj_clenshaw(szeta, czeta, f, k);
    let sdelta = delta.sin();
    let cdelta = delta.cos();
    let seta = szeta * cdelta + czeta * sdelta;
    let ceta = czeta * cdelta - szeta * sdelta;
    (seta, ceta)
}

/// Ported from src/latitudes.cpp (`pj_rectifying_radius` essence).
/// Returns the rectifying-radius scale factor as a function of the third
/// flattening `n`.
///
/// This is `(quarter meridian) / ((a+b)/2 * pi/2)` expanded as a series in
/// `n^2`; the coefficients are `((2*k - 3)!! / (2*k)!!)^2` for `k = 0..3`.
pub fn pj_rectifying_radius(n: f64) -> f64 {
    let coeff_rad = [1.0, 1.0 / 4.0, 1.0 / 64.0, 1.0 / 256.0];
    pj_polyval(n * n, &coeff_rad, 3) / (1.0 + n)
}

/// Ported from src/latitudes.cpp (`pj_auxlat_coeffs`).
///
/// Fills `f` (length 6) with the Fourier coefficients converting auxiliary
/// latitude `auxin` to `auxout` for the given third flattening `n`.
///
/// Returns `Err(ProjError::IllegalArgValue)` if the requested conversion is not
/// supported by the coefficient table (so callers never get silent wrong
/// numerics).
pub fn pj_auxlat_coeffs(n: f64, auxin: AuxLat, auxout: AuxLat, f: &mut [f64]) -> ProjResult<()> {
    let lmax = AuxLat::ORDER; // 6
    let auxin_i = auxin as usize;
    let auxout_i = auxout as usize;
    let k = AuxLat::NUMBER * auxout_i + auxin_i;
    // Bounds: k and k+1 must be valid PTRS indices.
    let o_start = match PTRS.get(k) {
        Some(v) => *v,
        None => return Err(ProjError::IllegalArgValue),
    };
    let o_next = match PTRS.get(k + 1) {
        Some(v) => *v,
        None => return Err(ProjError::IllegalArgValue),
    };
    if o_start == o_next {
        // Unsupported conversion.
        return Err(ProjError::IllegalArgValue);
    }
    let mut o = o_start;
    let mut d = n;
    let n2 = n * n;
    if auxin <= AuxLat::Rectifying && auxout <= AuxLat::Rectifying {
        // NOTE: "<=" here is by the numeric repr value: Geographic(0)..Rectifying(3).
        for l in 0..lmax {
            let m = (lmax - l - 1) / 2; // order of polynomial in n^2
            let coeff_slice = COEFFS.get(o..).unwrap_or(&[]);
            let fl = d * pj_polyval(n2, coeff_slice, m as isize);
            if let Some(slot) = f.get_mut(l) {
                *slot = fl;
            }
            o += m + 1;
            d *= n;
        }
    } else {
        for l in 0..lmax {
            let m = lmax - l - 1; // order of polynomial in n
            let coeff_slice = COEFFS.get(o..).unwrap_or(&[]);
            let fl = d * pj_polyval(n, coeff_slice, m as isize);
            if let Some(slot) = f.get_mut(l) {
                *slot = fl;
            }
            o += m + 1;
            d *= n;
        }
    }
    Ok(())
}

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

    /// WGS84 third flattening `n = f / (2 - f)` (~0.0016792203863836...).
    fn wgs84_n() -> f64 {
        let f = 1.0 / 298.257223563_f64;
        f / (2.0 - f)
    }

    #[test]
    fn phi_mu_phi_round_trip() {
        let n = wgs84_n();
        let mut c_fwd = [0.0; 6];
        let mut c_inv = [0.0; 6];
        pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Rectifying, &mut c_fwd)
            .expect("fwd coeffs");
        pj_auxlat_coeffs(n, AuxLat::Rectifying, AuxLat::Geographic, &mut c_inv)
            .expect("inv coeffs");
        for &phi in &[-1.4, -0.7, -0.1, 0.0, 0.1, 0.7, 1.4] {
            let mu = pj_auxlat_convert(phi, &c_fwd, 6);
            let back = pj_auxlat_convert(mu, &c_inv, 6);
            assert!((back - phi).abs() < 1e-12, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn phi_chi_phi_round_trip() {
        let n = wgs84_n();
        let mut c_fwd = [0.0; 6];
        let mut c_inv = [0.0; 6];
        pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Conformal, &mut c_fwd).expect("fwd coeffs");
        pj_auxlat_coeffs(n, AuxLat::Conformal, AuxLat::Geographic, &mut c_inv).expect("inv coeffs");
        for &phi in &[-1.4, -0.7, -0.1, 0.0, 0.1, 0.7, 1.4] {
            let chi = pj_auxlat_convert(phi, &c_fwd, 6);
            let back = pj_auxlat_convert(chi, &c_inv, 6);
            assert!((back - phi).abs() < 1e-12, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn unsupported_conversion_returns_err() {
        let n = wgs84_n();
        // Parametric=1, Geocentric=2; PTRS index 6*2+1=13 which has
        // PTRS[13]==54==PTRS[14] -> unsupported -> Err.
        assert!(
            pj_auxlat_coeffs(n, AuxLat::Parametric, AuxLat::Geocentric, &mut [0.0; 6]).is_err()
        );
    }

    #[test]
    fn rectifying_radius_at_zero() {
        assert_eq!(pj_rectifying_radius(0.0), 1.0);
    }

    #[test]
    fn polyval_negative_n_is_zero() {
        assert_eq!(pj_polyval(2.0, &[1.0, 2.0, 3.0], -1), 0.0);
    }

    #[test]
    fn polyval_basic_horner() {
        // p = 1 + 2x + 3x^2 at x = 2 -> 1 + 4 + 12 = 17.
        assert_eq!(pj_polyval(2.0, &[1.0, 2.0, 3.0], 2), 17.0);
    }

    #[test]
    fn clenshaw_zero_coeffs_is_zero() {
        let zeta = 0.5_f64;
        assert_eq!(pj_clenshaw(zeta.sin(), zeta.cos(), &[0.0; 6], 6), 0.0);
    }

    #[test]
    fn convert_sc_matches_scalar() {
        let n = wgs84_n();
        let mut c = [0.0; 6];
        pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Conformal, &mut c).expect("coeffs");
        let zeta = 0.7_f64;
        let a = pj_auxlat_convert(zeta, &c, 6);
        let b = pj_auxlat_convert_sc(zeta, zeta.sin(), zeta.cos(), &c, 6);
        assert_eq!(a, b);
    }

    #[test]
    fn convert_seta_ceta_consistent() {
        let n = wgs84_n();
        let mut c = [0.0; 6];
        pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Conformal, &mut c).expect("coeffs");
        let zeta = 0.7_f64;
        let eta = pj_auxlat_convert(zeta, &c, 6);
        let (seta, ceta) = pj_auxlat_convert_seta_ceta(zeta.sin(), zeta.cos(), &c, 6);
        assert!((seta - eta.sin()).abs() < 1e-12);
        assert!((ceta - eta.cos()).abs() < 1e-12);
    }

    #[test]
    fn auxlat_ordering_by_discriminant() {
        assert!(AuxLat::Geographic < AuxLat::Rectifying);
        assert!(AuxLat::Rectifying < AuxLat::Conformal);
        assert!(AuxLat::Conformal < AuxLat::Authalic);
        assert!(AuxLat::Geographic <= AuxLat::Rectifying);
    }

    #[test]
    fn phi_xi_phi_round_trip() {
        let n = wgs84_n();
        let mut c_fwd = [0.0; 6];
        let mut c_inv = [0.0; 6];
        pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Authalic, &mut c_fwd).expect("fwd coeffs");
        pj_auxlat_coeffs(n, AuxLat::Authalic, AuxLat::Geographic, &mut c_inv).expect("inv coeffs");
        for &phi in &[-1.4, -0.7, -0.1, 0.0, 0.1, 0.7, 1.4] {
            let xi = pj_auxlat_convert(phi, &c_fwd, 6);
            let back = pj_auxlat_convert(xi, &c_inv, 6);
            assert!((back - phi).abs() < 1e-12, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn chi_mu_chi_round_trip() {
        let n = wgs84_n();
        let mut c_fwd = [0.0; 6];
        let mut c_inv = [0.0; 6];
        pj_auxlat_coeffs(n, AuxLat::Conformal, AuxLat::Rectifying, &mut c_fwd).expect("fwd coeffs");
        pj_auxlat_coeffs(n, AuxLat::Rectifying, AuxLat::Conformal, &mut c_inv).expect("inv coeffs");
        for &chi in &[-1.4, -0.7, -0.1, 0.0, 0.1, 0.7, 1.4] {
            let mu = pj_auxlat_convert(chi, &c_fwd, 6);
            let back = pj_auxlat_convert(mu, &c_inv, 6);
            assert!((back - chi).abs() < 1e-12, "chi={} back={}", chi, back);
        }
    }
}