oxiproj-core 0.1.2

//! Exact distortion factors computed via automatic differentiation.
//!
//! [`FactorsExact`] is the AD-aware companion to the numeric [`crate::Factors`]
//! struct. While [`crate::Factors`] uses finite differences, [`FactorsExact`]
//! computes the exact Jacobian via forward-mode AD (`Dual1<2>`), giving
//! machine-precision scale factors and angular distortion.

/// Exact distortion factors computed via automatic differentiation.
///
/// Companion to the numeric [`crate::Factors`] struct (kept for PROJ parity).
/// Constructed via [`FactorsExact::from_jacobian`] after evaluating a
/// `ProjectGeneric` implementation with `Dual1<2>` inputs.
#[derive(Debug, Clone, Copy)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct FactorsExact {
    /// Meridional scale factor `h`. On the ellipsoid this is
    /// `‖∂(x,y)/∂φ‖ / (a·M)`, where `M = (1−es)/(1−es·sin²φ)^1.5` is the
    /// normalized meridian radius of curvature; on a sphere (`es == 0`, `M = 1`)
    /// it reduces to `‖∂(x,y)/∂φ‖ / a`.
    pub meridian_scale: f64,
    /// Parallel scale factor `k`. On the ellipsoid this is
    /// `‖∂(x,y)/∂λ‖ / (a·N·cos φ)`, where `N = 1/(1−es·sin²φ)^0.5` is the
    /// normalized prime-vertical radius of curvature; on a sphere (`es == 0`,
    /// `N = 1`) it reduces to `‖∂(x,y)/∂λ‖ / (a·cos φ)`.
    pub parallel_scale: f64,
    /// Areal scale factor `s = det J · r / (a²·cos φ)` with the PROJ ellipsoid
    /// factor `r = (1−es·sin²φ)² / (1−es)` (`r = 1` on a sphere). Matches PROJ's
    /// signed `factors.cpp` areal scale (`1` for equal-area projections).
    pub areal_scale: f64,
    /// Maximum angular distortion `ω = 2·arcsin((a′−b′)/(a′+b′))` in radians,
    /// where `a′,b′` are the Tissot semi-axes derived from `h`, `k` and `s`
    /// exactly as PROJ `factors.cpp` does. `0` for conformal projections.
    pub angular_distortion: f64,
    /// Meridian convergence γ in radians: angle from grid north to true north.
    pub meridian_convergence: f64,
    /// Partial derivative ∂x/∂λ (longitude).
    pub dx_dlam: f64,
    /// Partial derivative ∂x/∂φ (latitude).
    pub dx_dphi: f64,
    /// Partial derivative ∂y/∂λ (longitude).
    pub dy_dlam: f64,
    /// Partial derivative ∂y/∂φ (latitude).
    pub dy_dphi: f64,
    /// `true` when computed via exact AD; `false` if a finite-difference
    /// fallback was used instead.
    pub is_exact: bool,
}

impl FactorsExact {
    /// Compute distortion factors from the 2×2 Jacobian and geodetic latitude.
    ///
    /// The four partial derivatives are obtained by evaluating a
    /// `ProjectGeneric::project_fwd_generic` with `Dual1<2>` variables and
    /// reading the `.d` arrays:
    ///
    /// ```text
    /// lam_d = Dual1::<2>::variable(lam, 0);
    /// phi_d = Dual1::<2>::variable(phi, 1);
    /// let (x, y) = proj.project_fwd_generic(lam_d, phi_d)?;
    /// // dx_dlam = x.d[0], dx_dphi = x.d[1]
    /// // dy_dlam = y.d[0], dy_dphi = y.d[1]
    /// ```
    ///
    /// `a` is the semi-major axis in the projection's working units. When the
    /// `ProjectGeneric` forward runs in the engine's normalised (`a = 1`) space
    /// — as it always does — the caller **must** pass `a = 1.0`; the four
    /// partials are then already in unit-sphere units, matching PROJ's
    /// `pj_deriv`/`pj_factors` which operate before the semi-major-axis scaling
    /// applied in `fwd_finalize`.
    ///
    /// This is the **spherical** constructor (`es == 0`): the meridian and
    /// prime-vertical radii of curvature `M = N = 1`, so no ellipsoid
    /// radius-of-curvature normalization is applied. For a real ellipsoid use
    /// [`FactorsExact::from_jacobian_es`], which threads the squared
    /// eccentricity `es` and applies the same `M`/`N` correction PROJ's
    /// `factors.cpp` does.
    pub fn from_jacobian(
        dx_dlam: f64,
        dx_dphi: f64,
        dy_dlam: f64,
        dy_dphi: f64,
        phi: f64,
        a: f64,
    ) -> Self {
        Self::from_jacobian_es(dx_dlam, dx_dphi, dy_dlam, dy_dphi, phi, a, 0.0)
    }

    /// Ellipsoid-aware companion to [`FactorsExact::from_jacobian`].
    ///
    /// Extends the spherical constructor with the squared eccentricity `es`
    /// (`one_es = 1 − es` is derived internally) so the exact factors carry the
    /// same ellipsoid radius-of-curvature normalization PROJ 9.8.0
    /// `src/factors.cpp` (`pj_factors`) applies to the raw derivatives:
    ///
    /// ```text
    /// t = 1 − es·sin²φ
    /// n = √t
    /// h *= t·n / one_es      // divide meridional magnitude by M = one_es / t^1.5
    /// k *= n                 // divide parallel   magnitude by N = 1 / √t
    /// r  = t² / one_es        // areal-scale ellipsoid factor
    /// ```
    ///
    /// For a sphere (`es == 0`) `t = n = r = 1`, so this reduces **bit-for-bit**
    /// to the spherical formula: `h = ‖∂(x,y)/∂φ‖/a`,
    /// `k = ‖∂(x,y)/∂λ‖/(a·cos φ)`, `s = det J / (a²·cos φ)`.
    ///
    /// `a` obeys the same normalised-space contract as
    /// [`FactorsExact::from_jacobian`] (pass `a = 1.0` for the engine's
    /// `a = 1` `ProjectGeneric` forward).
    pub fn from_jacobian_es(
        dx_dlam: f64,
        dx_dphi: f64,
        dy_dlam: f64,
        dy_dphi: f64,
        phi: f64,
        a: f64,
        es: f64,
    ) -> Self {
        let one_es = 1.0 - es;
        let cos_phi = phi.cos().abs().max(1e-10);
        // Raw scale factors in normalised space (PROJ pj_factors, sphere first):
        //   h = ‖(dx_dphi, dy_dphi)‖ / a
        //   k = ‖(dx_dlam, dy_dlam)‖ / (a·cos φ)
        let mut h = (dx_dphi * dx_dphi + dy_dphi * dy_dphi).sqrt() / a;
        let mut k = (dx_dlam * dx_dlam + dy_dlam * dy_dlam).sqrt() / (a * cos_phi);
        // Jacobian determinant (PROJ: der.y_p·der.x_l − der.x_p·der.y_l).
        let det = dx_dlam * dy_dphi - dx_dphi * dy_dlam;
        // Ellipsoid radius-of-curvature correction (PROJ factors.cpp lines 95-103):
        // divide h by M = one_es / t^1.5 and k by N = 1/√t, and form the areal
        // factor r = t² / one_es. On a sphere (es == 0) this branch is skipped
        // and r = 1, keeping the spherical result bit-identical.
        let r = if es != 0.0 {
            let sin_phi = phi.sin();
            let t = 1.0 - es * sin_phi * sin_phi;
            let n = t.sqrt();
            h *= t * n / one_es;
            k *= n;
            t * t / one_es
        } else {
            1.0
        };
        // Areal scale (PROJ: s = (der.y_p·der.x_l − der.x_p·der.y_l)·r/cos φ),
        // kept signed to match PROJ and the numeric `Factors::compute` port.
        let s = det * r / (a * a * cos_phi);
        // Angular distortion via the Tissot semi-axes, matching PROJ's
        // factors.cpp exactly (rather than the orthogonal-graticule shortcut
        // 2·asin(|h−k|/(h+k)), which is wrong when the meridian/parallel images
        // are not perpendicular, e.g. oblique equal-area projections):
        //   t = k² + h²; a′ = ½(√(t+2s) + √(t−2s)); b′ = ½(√(t+2s) − √(t−2s))
        //   ω = 2·asin((a′−b′)/(a′+b′))
        let t_sq = k * k + h * h;
        let tmp_a = (t_sq + 2.0 * s).max(0.0).sqrt();
        let t_diff = t_sq - 2.0 * s;
        let t_diff_sqrt = if t_diff > 0.0 { t_diff.sqrt() } else { 0.0 };
        let big_b = 0.5 * (tmp_a - t_diff_sqrt);
        let big_a = 0.5 * (tmp_a + t_diff_sqrt);
        let omega = if (big_a + big_b) > 1e-14 {
            2.0 * ((big_a - big_b) / (big_a + big_b)).clamp(-1.0, 1.0).asin()
        } else {
            0.0
        };
        // Meridian convergence. PROJ `factors.cpp` computes
        // `fac->conv = -atan2(fac->der.x_p, fac->der.y_p)`; the leading minus is
        // REQUIRED for PROJ parity and to agree with the numeric
        // `Factors::compute` port (`let conv = -x_p.atan2(y_p)`). Without it the
        // exact convergence is sign-flipped — e.g. `lcc +ellps=WGS84 @ (5,39)`
        // reports −3.152° where PROJ (`proj -V`) reports +3.152°.
        let conv = -dx_dphi.atan2(dy_dphi);
        Self {
            meridian_scale: h,
            parallel_scale: k,
            areal_scale: s,
            angular_distortion: omega,
            meridian_convergence: conv,
            dx_dlam,
            dx_dphi,
            dy_dlam,
            dy_dphi,
            is_exact: true,
        }
    }
}

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

    #[test]
    fn spherical_mercator_at_equator() {
        // At phi=0: h=k=1, s=1, omega=0, conv=0
        // Spherical Mercator Jacobian at (lam, phi=0):
        //   x = lam => dx/dlam=1, dx/dphi=0
        //   y = ln(tan(pi/4 + phi/2)) => dy/dlam=0, dy/dphi=sec(0)=1
        let fe = FactorsExact::from_jacobian(1.0, 0.0, 0.0, 1.0, 0.0, 1.0);
        assert!(
            (fe.meridian_scale - 1.0).abs() < 1e-12,
            "h={}",
            fe.meridian_scale
        );
        assert!(
            (fe.parallel_scale - 1.0).abs() < 1e-12,
            "k={}",
            fe.parallel_scale
        );
        assert!((fe.areal_scale - 1.0).abs() < 1e-12, "s={}", fe.areal_scale);
        assert!(
            fe.angular_distortion.abs() < 1e-12,
            "omega={}",
            fe.angular_distortion
        );
        assert!(fe.is_exact);
    }

    #[test]
    fn spherical_mercator_at_phi30() {
        use core::f64::consts::FRAC_PI_6;
        let phi = FRAC_PI_6; // 30 deg
        let sec_phi = 1.0 / phi.cos();
        // Spherical Mercator Jacobian at phi=30deg (a=1):
        //   dx/dlam=1, dx/dphi=0, dy/dlam=0, dy/dphi=sec(phi)
        // h = sec(phi), k = 1/cos(phi) = sec(phi)  (conformal => h=k)
        // s = |det J| / (a^2 * cos phi) = sec(phi) / cos(phi) = sec^2(phi)
        //   Mercator is conformal but NOT equal-area; s != 1.
        let fe = FactorsExact::from_jacobian(1.0, 0.0, 0.0, sec_phi, phi, 1.0);
        assert!(
            (fe.meridian_scale - sec_phi).abs() < 1e-12,
            "h={}",
            fe.meridian_scale
        );
        assert!(
            (fe.parallel_scale - sec_phi).abs() < 1e-12,
            "k={}",
            fe.parallel_scale
        );
        let expected_s = sec_phi * sec_phi;
        assert!(
            (fe.areal_scale - expected_s).abs() < 1e-12,
            "s={}",
            fe.areal_scale
        );
        // ω is computed from the Tissot semi-axes exactly as PROJ's factors.cpp
        // does. For a conformal projection (h == k) that path suffers the same
        // catastrophic cancellation PROJ does — `t − 2s = h²+k²−2s` is ~0 up to
        // one ULP because `s` is formed from the determinant, not from `h·k` —
        // so ω is ~1.8e-8 rad, NOT exactly 0. PROJ itself reports this: `proj -S
        // +proj=merc` prints ω ≈ 1.05e-6°, i.e. ~1.8e-8 rad. The bound is
        // therefore PROJ-consistent (was 1e-12 under the old, non-PROJ
        // 2·asin(|h−k|/(h+k)) shortcut, which cannot render laea's ω correctly).
        assert!(
            fe.angular_distortion.abs() < 1e-6,
            "omega={}",
            fe.angular_distortion
        );
    }

    #[test]
    fn from_jacobian_es_reduces_to_sphere_when_es_zero() {
        // With es == 0 the ellipsoid-aware constructor must be bit-for-bit
        // identical to the spherical `from_jacobian`, for an arbitrary Jacobian.
        use core::f64::consts::FRAC_PI_6;
        let (dxl, dxp, dyl, dyp, phi, a) = (1.3, -0.4, 0.2, 1.7, FRAC_PI_6, 1.0);
        let sphere = FactorsExact::from_jacobian(dxl, dxp, dyl, dyp, phi, a);
        let es0 = FactorsExact::from_jacobian_es(dxl, dxp, dyl, dyp, phi, a, 0.0);
        assert_eq!(
            sphere.meridian_scale.to_bits(),
            es0.meridian_scale.to_bits()
        );
        assert_eq!(
            sphere.parallel_scale.to_bits(),
            es0.parallel_scale.to_bits()
        );
        assert_eq!(sphere.areal_scale.to_bits(), es0.areal_scale.to_bits());
        assert_eq!(
            sphere.angular_distortion.to_bits(),
            es0.angular_distortion.to_bits()
        );
        assert_eq!(
            sphere.meridian_convergence.to_bits(),
            es0.meridian_convergence.to_bits()
        );
    }

    #[test]
    fn ellipsoidal_mercator_matches_proj() {
        // Ellipsoidal Mercator (normalised, a=1, k0=1) at φ=30°:
        //   x = λ                    => dx/dλ = 1,  dx/dφ = 0
        //   y = asinh(tanφ) − e·atanh(e·sinφ)
        //       => dy/dλ = 0,  dy/dφ = one_es / (cosφ·(1 − es·sin²φ))
        // After the ellipsoid correction PROJ reports h = k = √t / cosφ with
        // t = 1 − es·sin²φ; areal scale s = h·k (conformal). Cross-checked
        // against `proj -V +proj=merc +ellps=WGS84` at (10,30): h = k = 1.15373388.
        use core::f64::consts::FRAC_PI_6;
        let phi = FRAC_PI_6; // 30°
                             // WGS84 squared eccentricity.
        let es = 0.006_694_379_990_141_316_5_f64;
        let one_es = 1.0 - es;
        let sin_phi = phi.sin();
        let cos_phi = phi.cos();
        let t = 1.0 - es * sin_phi * sin_phi;
        let dy_dphi = one_es / (cos_phi * t);
        let fe = FactorsExact::from_jacobian_es(1.0, 0.0, 0.0, dy_dphi, phi, 1.0, es);

        let expected_hk = t.sqrt() / cos_phi;
        assert!(
            (fe.meridian_scale - expected_hk).abs() < 1e-12,
            "h={} expected {}",
            fe.meridian_scale,
            expected_hk
        );
        assert!(
            (fe.parallel_scale - expected_hk).abs() < 1e-12,
            "k={} expected {}",
            fe.parallel_scale,
            expected_hk
        );
        // Anchor to the PROJ 9.x reference value.
        assert!(
            (fe.meridian_scale - 1.153_733_88).abs() < 1e-7,
            "h={} must match proj -V 1.15373388",
            fe.meridian_scale
        );
        // Conformal: h == k so ω ≈ 0 (PROJ-consistent Tissot cancellation puts
        // it at ~1.8e-8 rad, not exactly 0 — see spherical_mercator_at_phi30),
        // and areal scale s = h·k.
        assert!(
            fe.angular_distortion.abs() < 1e-6,
            "ellipsoidal Mercator is conformal: ω={}",
            fe.angular_distortion
        );
        assert!(
            (fe.areal_scale - expected_hk * expected_hk).abs() < 1e-12,
            "s={} expected h·k={}",
            fe.areal_scale,
            expected_hk * expected_hk
        );
    }
}