oxiproj-core 0.1.2

//! Gaussian sphere intermediary projection ported from PROJ `src/gauss.cpp`.
//!
//! The Gaussian sphere is a conformal intermediary used by several PROJ
//! projections (e.g. oblique stereographic). It maps geodetic latitude on the
//! ellipsoid to the conformal latitude on a sphere via the parameters prepared
//! in [`pj_gauss_ini`], then [`pj_gauss`] / [`pj_inv_gauss`] convert between
//! ellipsoidal and conformal-sphere coordinates.

use crate::consts::{M_FORTPI, M_HALFPI};
use crate::coord::Lp;
use crate::error::{ProjError, ProjResult};

/// Ported from src/gauss.cpp (the `GAUSS` / `pj_gauss_data` struct).
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Gauss {
    c: f64,
    k: f64,
    e: f64,
    ratexp: f64,
}

/// Ported from src/gauss.cpp (`srat`).
fn srat(esinp: f64, ratexp: f64) -> f64 {
    ((1.0 - esinp) / (1.0 + esinp)).powf(ratexp)
}

/// Ported from src/gauss.cpp (`pj_gauss_ini`).
///
/// Returns the prepared `Gauss` state plus `chi` (the conformal latitude at the
/// origin) and `rc` (the radius of conformal sphere scale), matching the two
/// out-parameters of the C function.
pub fn pj_gauss_ini(e: f64, phi0: f64) -> ProjResult<(Gauss, f64, f64)> {
    let es = e * e;
    let sphi = phi0.sin();
    let mut cphi = phi0.cos();
    cphi *= cphi;
    let rc = (1.0 - es).sqrt() / (1.0 - es * sphi * sphi);
    let c = (1.0 + es * cphi * cphi / (1.0 - es)).sqrt();
    if c == 0.0 {
        return Err(ProjError::IllegalArgValue);
    }
    let chi = (sphi / c).asin();
    let ratexp = 0.5 * c * e;
    let srat_val = srat(e * sphi, ratexp);
    if srat_val == 0.0 {
        return Err(ProjError::IllegalArgValue);
    }
    let k = if 0.5 * phi0 + M_FORTPI < 1e-10 {
        1.0 / srat_val
    } else {
        (0.5 * chi + M_FORTPI).tan() / ((0.5 * phi0 + M_FORTPI).tan().powf(c) * srat_val)
    };
    Ok((Gauss { c, k, e, ratexp }, chi, rc))
}

/// Ported from src/gauss.cpp (`pj_gauss`).
pub fn pj_gauss(elp: Lp, g: &Gauss) -> Lp {
    let phi = 2.0
        * (g.k * (0.5 * elp.phi + M_FORTPI).tan().powf(g.c) * srat(g.e * elp.phi.sin(), g.ratexp))
            .atan()
        - M_HALFPI;
    let lam = g.c * elp.lam;
    Lp::new(lam, phi)
}

/// Ported from src/gauss.cpp (`pj_inv_gauss`).
///
/// Faithful note: PROJ records a context errno if the fixed-point iteration
/// fails to converge within `MAX_ITER`; this port returns `Err(NoConvergence)`
/// in that case (PROJ still returns the last value, but we surface the failure).
pub fn pj_inv_gauss(slp: Lp, g: &Gauss) -> ProjResult<Lp> {
    const MAX_ITER: usize = 20;
    const DEL_TOL: f64 = 1e-14;
    let lam = slp.lam / g.c;
    let num = ((0.5 * slp.phi + M_FORTPI).tan() / g.k).powf(1.0 / g.c);
    let mut phi_in = slp.phi;
    let mut phi_out = phi_in;
    let mut converged = false;
    let mut iter = MAX_ITER;
    while iter > 0 {
        phi_out = 2.0 * (num * srat(g.e * phi_in.sin(), -0.5 * g.e)).atan() - M_HALFPI;
        if (phi_out - phi_in).abs() < DEL_TOL {
            converged = true;
            break;
        }
        phi_in = phi_out;
        iter -= 1;
    }
    if !converged {
        return Err(ProjError::NoConvergence);
    }
    Ok(Lp::new(lam, phi_out))
}

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

    /// WGS84 eccentricity used across the round-trip checks.
    fn wgs84_e() -> f64 {
        let es = 0.0066943799901413165_f64;
        es.sqrt()
    }

    #[test]
    fn round_trip_forward_inverse() {
        let e = wgs84_e();
        let (g, _chi, _rc) = pj_gauss_ini(e, 0.5).expect("ini");
        let cases = [
            (0.1_f64, 0.4_f64),
            (-0.2, 0.6),
            (0.0, 0.5),
            (0.3, 0.55),
            (-0.1, 0.45),
        ];
        for (lam, phi) in cases {
            let s = pj_gauss(Lp::new(lam, phi), &g);
            let back = pj_inv_gauss(s, &g).expect("inv");
            assert!(
                (back.lam - lam).abs() < 1e-12 && (back.phi - phi).abs() < 1e-12,
                "lam={} phi={} back=({},{})",
                lam,
                phi,
                back.lam,
                back.phi
            );
        }
    }

    #[test]
    fn ini_outputs_are_sane() {
        let e = wgs84_e();
        let (_g, chi, rc) = pj_gauss_ini(e, 0.5).expect("ini");
        assert!(rc > 0.0, "rc should be positive, got {}", rc);
        assert!(chi.is_finite(), "chi should be finite, got {}", chi);
    }
}