oxiproj-core 0.1.2

//! Authalic (equal-area) latitude helpers.
//! Ported from PROJ `src/latitudes.cpp` (`pj_authalic_lat*`).
//!
//! Provides Snyder's authalic-area quantity `q` (eq. 3-12) and the forward /
//! inverse authalic-latitude conversions used by the equal-area projections
//! (`aea`, `laea`, `cea`). For `|n| < 0.01` the Karney series (Clenshaw) is
//! used; otherwise the Snyder closed form (`asin(q/qp)`) with Newton refinement
//! on the inverse.

use crate::error::ProjResult;
use crate::math::aux::{pj_auxlat_coeffs, pj_auxlat_convert, pj_auxlat_convert_sc, AuxLat};
use crate::scalar::Scalar;

/// Threshold on the third flattening below which the Karney series is used.
/// Ported from `PROJ_AUTHALIC_SERIES_VALID(n)` in src/latitudes.cpp.
const PROJ_AUTHALIC_SERIES_CUTOFF: f64 = 0.01;

/// `true` when the auxiliary-latitude series is accurate enough to be preferred
/// over the analytical (closed-form) expressions. Ported from
/// `PROJ_AUTHALIC_SERIES_VALID(n) == (fabs(n) < 0.01)`.
fn series_valid(n: f64) -> bool {
    n.abs() < PROJ_AUTHALIC_SERIES_CUTOFF
}

/// Precomputed authalic-latitude coefficients (the "APA" bundle in PROJ).
///
/// Bundles the Fourier coefficients for both directions of the authalic
/// conversion together with the ellipsoid scalars needed by the closed-form
/// (`asin(q/qp)`) branch and the Newton refinement of the inverse, so callers
/// build this once per ellipsoid and reuse it per point.
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Apa {
    /// Authalic->geographic Fourier coefficients (used by the inverse). Always filled.
    pub inv: [f64; 6],
    /// Geographic->authalic Fourier coefficients (used by the forward series path).
    /// Only meaningful when [`Apa::series_valid`] is `true`.
    pub fwd: [f64; 6],
    /// `true` when `|n| < 0.01` and the series path is used.
    pub series_valid: bool,
    /// Third flattening `n` of the ellipsoid.
    pub n: f64,
    /// Eccentricity `e`.
    pub e: f64,
    /// Eccentricity squared `es = e^2`.
    pub es: f64,
    /// `one_es = 1 - e^2`.
    pub one_es: f64,
}

/// Snyder authalic `q` (eq. 3-12). Computes the authalic-area quantity for a
/// given `sin(phi)`. `e` is the eccentricity, `one_es == 1 - e^2`.
///
/// For a sphere (`e < 1e-7`) this is `2*sin(phi)`. `qp` (the value of `q` at the
/// pole) is obtained as `pj_qsfn(1.0, e, one_es)`.
pub fn pj_qsfn(sinphi: f64, e: f64, one_es: f64) -> f64 {
    const EPSILON: f64 = 1e-7;
    if e >= EPSILON {
        let e_sinphi = e * sinphi;
        let one_minus_e_sinphi_sq = 1.0 - e_sinphi * e_sinphi;
        // avoid zero division, fail gracefully
        if one_minus_e_sinphi_sq == 0.0 {
            return f64::INFINITY;
        }
        one_es * (sinphi / one_minus_e_sinphi_sq + e_sinphi.atanh() / e)
    } else {
        2.0 * sinphi
    }
}

/// Generic version of [`pj_qsfn`] for automatic-differentiation support.
///
/// Ellipsoid constants (`e`, `one_es`) stay as `f64`. `atanh` is expressed via
/// `atanh(x) = 0.5 * ln((1+x)/(1-x))` since [`Scalar`] does not provide `atanh`.
pub fn pj_qsfn_g<S: Scalar>(sinphi: S, e: f64, one_es: f64) -> S {
    const EPSILON: f64 = 1e-7;
    let one_es_s = S::from_f64(one_es);
    if e >= EPSILON {
        let e_s = S::from_f64(e);
        let e_sinphi = e_s * sinphi;
        let one = S::one();
        let one_minus = one - e_sinphi * e_sinphi;
        // atanh(e*sinphi) = 0.5 * ln((1 + e*sinphi) / (1 - e*sinphi))
        let half = S::from_f64(0.5);
        let atanh_e_sinphi = half * ((one + e_sinphi) / (one - e_sinphi)).ln();
        one_es_s * (sinphi / one_minus + atanh_e_sinphi / e_s)
    } else {
        S::from_f64(2.0) * sinphi
    }
}

/// Build the authalic coefficient bundle [`Apa`] from the ellipsoid scalars.
/// Ported from `pj_authalic_lat_compute_coeffs`.
///
/// `n` = third flattening, `e` = eccentricity, `es = e^2`, `one_es = 1 - e^2`.
/// Propagates any error from [`pj_auxlat_coeffs`].
pub fn pj_authalic_lat_compute_coeffs(n: f64, e: f64, es: f64, one_es: f64) -> ProjResult<Apa> {
    let valid = series_valid(n);
    let mut inv = [0.0_f64; 6];
    pj_auxlat_coeffs(n, AuxLat::Authalic, AuxLat::Geographic, &mut inv)?;
    let mut fwd = [0.0_f64; 6];
    if valid {
        pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Authalic, &mut fwd)?;
    }
    Ok(Apa {
        inv,
        fwd,
        series_valid: valid,
        n,
        e,
        es,
        one_es,
    })
}

/// Forward: geodetic `phi` -> authalic latitude `xi`. Ported from
/// `pj_authalic_lat`. `sinphi`/`cosphi` are `sin(phi)`/`cos(phi)`; `qp` is the
/// pole value `pj_qsfn(1.0, apa.e, apa.one_es)`.
pub fn pj_authalic_lat(phi: f64, sinphi: f64, cosphi: f64, apa: &Apa, qp: f64) -> f64 {
    if apa.series_valid {
        pj_auxlat_convert_sc(phi, sinphi, cosphi, &apa.fwd, 6)
    } else {
        let q = pj_qsfn(sinphi, apa.e, apa.one_es);
        let mut ratio = q / qp;
        if ratio.abs() > 1.0 {
            ratio = if ratio > 0.0 { 1.0 } else { -1.0 };
        }
        ratio.asin()
    }
}

/// Inverse: authalic `beta` -> geodetic `phi`. Ported from
/// `pj_authalic_lat_inverse`. For `|n| >= 0.01` performs up to 10 Newton
/// iterations. `qp` is `pj_qsfn(1.0, apa.e, apa.one_es)`.
pub fn pj_authalic_lat_inverse(beta: f64, apa: &Apa, qp: f64) -> f64 {
    let mut phi = pj_auxlat_convert(beta, &apa.inv, 6);
    if apa.series_valid {
        return phi;
    }
    let q = beta.sin() * qp;
    let q_div_one_minus_es = q / apa.one_es;
    for _ in 0..10 {
        let sinphi = phi.sin();
        let cosphi = phi.cos();
        let one_minus_es_sin2phi = 1.0 - apa.es * (sinphi * sinphi);
        let dphi = (one_minus_es_sin2phi * one_minus_es_sin2phi) / (2.0 * cosphi)
            * (q_div_one_minus_es
                - sinphi / one_minus_es_sin2phi
                - (apa.e * sinphi).atanh() / apa.e);
        if dphi.abs() >= 1e-15 {
            phi += dphi;
        } else {
            break;
        }
    }
    phi
}

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

    #[test]
    fn q_at_equator_is_zero() {
        let ell = Ellipsoid::named("WGS84").expect("wgs84");
        assert_eq!(pj_qsfn(0.0, ell.e, ell.one_es), 0.0);
    }

    #[test]
    fn qp_positive() {
        let ell = Ellipsoid::named("WGS84").expect("wgs84");
        let qp = pj_qsfn(1.0, ell.e, ell.one_es);
        assert!(qp > 0.0);
        assert!(qp.is_finite());
        assert!(qp > 1.99 && qp < 2.0);
    }

    #[test]
    fn authalic_round_trip() {
        let ell = Ellipsoid::named("WGS84").expect("wgs84");
        let apa = pj_authalic_lat_compute_coeffs(ell.n, ell.e, ell.es, ell.one_es).expect("apa");
        let qp = pj_qsfn(1.0, ell.e, ell.one_es);
        for &phi in &[0.1, 0.5, 1.0, -0.7, 0.0, 1.4, -1.3] {
            let xi = pj_authalic_lat(phi, phi.sin(), phi.cos(), &apa, qp);
            let back = pj_authalic_lat_inverse(xi, &apa, qp);
            assert!((back - phi).abs() < 1e-12, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn closed_form_round_trip_large_n() {
        let f = 1.0_f64 / 30.0;
        let es = 2.0 * f - f * f;
        let e = es.sqrt();
        let one_es = 1.0 - es;
        let n = f / (2.0 - f);
        let apa = pj_authalic_lat_compute_coeffs(n, e, es, one_es).expect("apa");
        assert!(!apa.series_valid);
        let qp = pj_qsfn(1.0, e, one_es);
        for &phi in &[0.1, 0.5, 1.0, -0.7, 1.3] {
            let xi = pj_authalic_lat(phi, phi.sin(), phi.cos(), &apa, qp);
            let back = pj_authalic_lat_inverse(xi, &apa, qp);
            assert!((back - phi).abs() < 1e-12, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn sphere_identity() {
        assert_eq!(pj_qsfn(0.3, 0.0, 1.0), 0.6);
        let apa = pj_authalic_lat_compute_coeffs(0.0, 0.0, 0.0, 1.0).expect("apa");
        let qp = pj_qsfn(1.0, 0.0, 1.0);
        for &phi in &[0.1, 0.5, -0.7, 1.0] {
            let xi = pj_authalic_lat(phi, phi.sin(), phi.cos(), &apa, qp);
            assert!((xi - phi).abs() < 1e-13, "phi={} xi={}", phi, xi);
            let back = pj_authalic_lat_inverse(xi, &apa, qp);
            assert!((back - phi).abs() < 1e-13, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn qsfn_pole_guard() {
        let ell = Ellipsoid::named("WGS84").expect("wgs84");
        let qp = pj_qsfn(1.0, ell.e, ell.one_es);
        assert!(qp.is_finite());
    }
}