oxiproj-core 0.1.2

//! Map-scale factors ported from PROJ 9.8.0 src/factors.cpp and src/deriv.cpp.

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

/// Map-scale factors at a geographic point.
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Factors {
    /// Meridional scale factor h (linear scale along meridian).
    pub h: f64,
    /// Parallel scale factor k (linear scale along parallel).
    pub k: f64,
    /// Areal scale factor s. Equal to 1.0 for equal-area projections.
    pub s: f64,
    /// Angular distortion ω in radians. 0 for conformal projections.
    pub omega: f64,
    /// Tissot semi-major axis.
    pub big_a: f64,
    /// Tissot semi-minor axis.
    pub big_b: f64,
    /// Meridian-parallel angle θ' (radians).
    pub theta_prime: f64,
    /// Meridian convergence γ (radians): angle between geographic north and grid north.
    pub conv: f64,
    /// dx/dlam — partial derivative of easting w.r.t. longitude.
    pub dx_dlam: f64,
    /// dx/dphi — partial derivative of easting w.r.t. latitude.
    pub dx_dphi: f64,
    /// dy/dlam — partial derivative of northing w.r.t. longitude.
    pub dy_dlam: f64,
    /// dy/dphi — partial derivative of northing w.r.t. latitude.
    pub dy_dphi: f64,
}

impl Factors {
    /// Compute map-scale factors from the raw projection forward function.
    ///
    /// `phi` is geodetic latitude in radians, `lam` is longitude relative to the
    /// central meridian in radians. `es` is the ellipsoid's eccentricity squared
    /// (0.0 for spherical). `one_es` is 1.0 - es.
    ///
    /// `pj_fwd` is a closure calling the raw projection forward in normalized
    /// space (a=1): takes (lam, phi) and returns (x, y).
    ///
    /// Ported from PROJ 9.8.0 src/deriv.cpp (pj_deriv) and src/factors.cpp (pj_factors).
    pub fn compute<F>(phi: f64, lam: f64, es: f64, one_es: f64, pj_fwd: F) -> ProjResult<Factors>
    where
        F: Fn(f64, f64) -> ProjResult<(f64, f64)>,
    {
        const DEFAULT_H: f64 = 1e-5;
        const EPS: f64 = 1e-12;

        // Guard: check latitude not over-range
        if (phi.abs() - core::f64::consts::FRAC_PI_2) > EPS {
            return Err(ProjError::InvalidCoord);
        }

        // If lat + h would overshoot pole, clamp it back
        let phi = if phi.abs() > (core::f64::consts::FRAC_PI_2 - DEFAULT_H) {
            if phi < 0.0 {
                -(core::f64::consts::FRAC_PI_2 - DEFAULT_H)
            } else {
                core::f64::consts::FRAC_PI_2 - DEFAULT_H
            }
        } else {
            phi
        };

        let h = DEFAULT_H;

        // PROJ pj_deriv: evaluate at 4 corner points
        // Point 1: (lam+h, phi+h)
        let (x1, y1) = pj_fwd(lam + h, phi + h)?;
        // Point 2: (lam+h, phi-h)
        let (x2, y2) = pj_fwd(lam + h, phi - h)?;
        // Point 3: (lam-h, phi-h)
        let (x3, y3) = pj_fwd(lam - h, phi - h)?;
        // Point 4: (lam-h, phi+h)
        let (x4, y4) = pj_fwd(lam - h, phi + h)?;

        // Accumulation from deriv.cpp: divide by 4h
        let inv_4h = 1.0 / (4.0 * h);
        let x_l = (x1 + x2 - x3 - x4) * inv_4h; // dx/dlam
        let y_l = (y1 + y2 - y3 - y4) * inv_4h; // dy/dlam
        let x_p = (x1 - x2 - x3 + x4) * inv_4h; // dx/dphi
        let y_p = (y1 - y2 - y3 + y4) * inv_4h; // dy/dphi

        let cos_phi = phi.cos();

        // Scale factors (sphere first)
        let mut h_scale = x_p.hypot(y_p);
        let mut k_scale = x_l.hypot(y_l) / cos_phi;

        // Ellipsoid correction (factors.cpp lines 95-103)
        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_scale *= t * n / one_es;
            k_scale *= n;
            t * t / one_es
        } else {
            1.0
        };

        // Convergence (factors.cpp line 106)
        let conv = -x_p.atan2(y_p);

        // Areal scale (factors.cpp lines 109-110)
        let s = (y_p * x_l - x_p * y_l) * r / cos_phi;

        // Meridian-parallel angle theta_prime (factors.cpp line 113)
        // clamp to [-1,1] before asin (mirrors PROJ's aasin)
        let ratio = if (h_scale * k_scale).abs() < f64::EPSILON {
            0.0
        } else {
            (s / (h_scale * k_scale)).clamp(-1.0, 1.0)
        };
        let theta_prime = ratio.asin();

        // Tissot ellipse axes (factors.cpp lines 116-121)
        let t_sq = k_scale * k_scale + h_scale * h_scale;
        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);

        // Angular distortion (factors.cpp lines 124-125)
        let omega_ratio = if (big_a + big_b).abs() < f64::EPSILON {
            0.0
        } else {
            ((big_a - big_b) / (big_a + big_b)).clamp(-1.0, 1.0)
        };
        let omega = 2.0 * omega_ratio.asin();

        Ok(Factors {
            h: h_scale,
            k: k_scale,
            s,
            omega,
            big_a,
            big_b,
            theta_prime,
            conv,
            dx_dlam: x_l,
            dx_dphi: x_p,
            dy_dlam: y_l,
            dy_dphi: y_p,
        })
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use core::f64::consts::PI;

    /// Simple sphere Mercator forward: x = lam, y = ln(tan(pi/4 + phi/2))
    fn sphere_merc_fwd(lam: f64, phi: f64) -> ProjResult<(f64, f64)> {
        let y = (PI / 4.0 + phi / 2.0).tan().ln();
        Ok((lam, y))
    }

    #[test]
    fn sphere_merc_at_equator_h_eq_k_eq_1() {
        // For sphere Mercator at equator: h=k=1, s=1, omega=0
        let f = Factors::compute(0.0, 0.0, 0.0, 1.0, sphere_merc_fwd).unwrap();
        assert!(
            (f.h - 1.0).abs() < 1e-6,
            "h should be 1.0 at equator, got {}",
            f.h
        );
        assert!(
            (f.k - 1.0).abs() < 1e-6,
            "k should be 1.0 at equator, got {}",
            f.k
        );
        assert!(
            (f.s - 1.0).abs() < 1e-6,
            "s should be 1.0 (conformal), got {}",
            f.s
        );
        assert!(
            f.omega.abs() < 1e-6,
            "omega should be ~0 (conformal), got {}",
            f.omega
        );
    }

    #[test]
    fn sphere_merc_at_lat_30_conformal() {
        // For sphere Mercator: h=k=1/cos(phi), which is the secant
        let phi = 30.0_f64.to_radians();
        let f = Factors::compute(phi, 0.0, 0.0, 1.0, sphere_merc_fwd).unwrap();
        let expected_scale = 1.0 / phi.cos();
        assert!(
            (f.h - expected_scale).abs() < 1e-5,
            "h={} expected {}",
            f.h,
            expected_scale
        );
        assert!(
            (f.k - expected_scale).abs() < 1e-5,
            "k={} expected {}",
            f.k,
            expected_scale
        );
        // Conformal: omega should be ~0
        assert!(f.omega.abs() < 1e-6, "omega should be ~0 for conformal");
    }

    #[test]
    fn near_pole_clamped() {
        // Should not error out near poles (latitude clamped)
        let phi = core::f64::consts::FRAC_PI_2 - 1e-6;
        let result = Factors::compute(phi, 0.0, 0.0, 1.0, sphere_merc_fwd);
        // May succeed or fail depending on projection domain
        let _ = result;
    }

    #[test]
    fn at_pole_returns_error() {
        let phi = core::f64::consts::FRAC_PI_2 + 1e-10;
        let result = Factors::compute(phi, 0.0, 0.0, 1.0, sphere_merc_fwd);
        assert!(result.is_err(), "Latitude beyond pole should be error");
    }
}