oxiproj-core 0.1.2

//! Tissot indicatrix: local distortion ellipse at a map point.
//! Formulas from Snyder (1987) "Map Projections — A Working Manual", §4, pp. 20–25.

use crate::FactorsExact;

/// Parameters describing the Tissot indicatrix at a point on the map.
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct TissotParams {
    /// Semi-major axis of the indicatrix ellipse (dimensionless scale factor).
    pub semi_major: f64,
    /// Semi-minor axis of the indicatrix ellipse (dimensionless scale factor).
    pub semi_minor: f64,
    /// Meridian scale factor h: ratio of projected meridian element to true.
    pub meridian_scale: f64,
    /// Parallel scale factor k: ratio of projected parallel element to true.
    pub parallel_scale: f64,
    /// Maximum angular distortion ω in radians. 0 = conformal.
    pub angular_distortion: f64,
    /// Areal scale factor s = a' · b'. 1 = equal-area.
    pub areal_scale: f64,
    /// Meridian convergence γ (angle between projected meridian and true north).
    pub meridian_convergence: f64,
    /// Whether the projection is locally conformal at this point (|ω| < 1e-9).
    pub is_conformal: bool,
    /// Whether the projection is locally equal-area at this point (|s - 1| < 1e-9).
    pub is_equal_area: bool,
}

impl TissotParams {
    /// Compute Tissot parameters from exact distortion factors.
    ///
    /// Uses Snyder (1987) §4-14 formulas:
    /// - `a' = 0.5 · (√(h²+k²+2s) + √(h²+k²−2s))`
    /// - `b' = 0.5 · (√(h²+k²+2s) − √(h²+k²−2s))`
    ///
    /// where `s = h·k·sin(θ')` is the areal scale factor.
    pub fn from_exact(h: f64, k: f64, areal_s: f64, _theta_prime: f64, conv: f64) -> Self {
        // Snyder (1987) eq. 4-14:
        // a' = 0.5 · (√(h²+k²+2s) + √(h²+k²−2s))
        // b' = 0.5 · (√(h²+k²+2s) − √(h²+k²−2s))
        let hk2 = h * h + k * k;
        let two_s = 2.0 * areal_s.abs();
        let sq_plus = (hk2 + two_s).max(0.0).sqrt();
        let sq_minus = (hk2 - two_s).max(0.0).sqrt();
        let semi_major = 0.5 * (sq_plus + sq_minus);
        let semi_minor = 0.5 * (sq_plus - sq_minus).abs();

        // Max angular distortion: 2·arcsin((a'−b')/(a'+b'))
        let sum_ab = semi_major + semi_minor;
        let angular_distortion = if sum_ab > 1e-15 {
            2.0 * ((semi_major - semi_minor) / sum_ab).clamp(-1.0, 1.0).asin()
        } else {
            0.0
        };

        let tol = 1e-9;
        TissotParams {
            semi_major,
            semi_minor,
            meridian_scale: h,
            parallel_scale: k,
            angular_distortion,
            areal_scale: semi_major * semi_minor,
            meridian_convergence: conv,
            is_conformal: angular_distortion.abs() < tol,
            is_equal_area: (semi_major * semi_minor - 1.0).abs() < tol,
        }
    }

    /// Compute Tissot indicatrix from a [`FactorsExact`] struct.
    ///
    /// Convenience wrapper around [`TissotParams::from_exact`] that extracts
    /// `h`, `k`, `areal_scale`, and `meridian_convergence` from the AD result.
    ///
    /// `theta_prime` is taken as `π/2` (the intersection angle between
    /// meridians and parallels, which equals π/2 for all standard projections
    /// in the spherical Jacobian formulation).
    pub fn from_factors(f: &FactorsExact) -> Self {
        use core::f64::consts::FRAC_PI_2;
        Self::from_exact(
            f.meridian_scale,
            f.parallel_scale,
            f.areal_scale,
            FRAC_PI_2,
            f.meridian_convergence,
        )
    }

    /// Eccentricity of the indicatrix ellipse (0 = circle = conformal).
    pub fn eccentricity(&self) -> f64 {
        if self.semi_major < 1e-15 {
            return 0.0;
        }
        let ratio = self.semi_minor / self.semi_major;
        (1.0 - ratio * ratio).max(0.0).sqrt()
    }

    /// Flattening of the indicatrix ellipse.
    pub fn flattening(&self) -> f64 {
        if self.semi_major < 1e-15 {
            return 0.0;
        }
        (self.semi_major - self.semi_minor) / self.semi_major
    }
}

impl core::fmt::Display for TissotParams {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        write!(
            f,
            "Tissot(a'={:.6}, b'={:.6}, h={:.6}, k={:.6}, ω={:.4}°, s={:.6})",
            self.semi_major,
            self.semi_minor,
            self.meridian_scale,
            self.parallel_scale,
            self.angular_distortion.to_degrees(),
            self.areal_scale
        )
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    #[cfg(feature = "no_std")]
    use alloc::format;
    use core::f64::consts::PI;

    #[test]
    fn mercator_is_conformal() {
        // Mercator at φ=0: h = k = 1, s = 1, θ' = π/2
        let t = TissotParams::from_exact(1.0, 1.0, 1.0, PI / 2.0, 0.0);
        assert!((t.semi_major - 1.0).abs() < 1e-9, "a' = {}", t.semi_major);
        assert!((t.semi_minor - 1.0).abs() < 1e-9, "b' = {}", t.semi_minor);
        assert!(
            t.angular_distortion.abs() < 1e-9,
            "ω = {}",
            t.angular_distortion
        );
        assert!(t.is_conformal);
        assert!((t.areal_scale - 1.0).abs() < 1e-9);
    }

    #[test]
    fn mercator_at_45_deg() {
        // Mercator at φ=45°: h = k = sec(45°) ≈ 1.4142, s = h*k = 2.0
        let sec45 = 1.0_f64 / (PI / 4.0).cos();
        let t = TissotParams::from_exact(sec45, sec45, sec45 * sec45, PI / 2.0, 0.0);
        assert!(
            (t.semi_major - sec45).abs() < 1e-6,
            "a'={} sec45={}",
            t.semi_major,
            sec45
        );
        assert!(
            (t.semi_minor - sec45).abs() < 1e-6,
            "b'={} sec45={}",
            t.semi_minor,
            sec45
        );
        assert!(t.is_conformal, "ω={}", t.angular_distortion);
        assert!((t.areal_scale - sec45 * sec45).abs() < 1e-6);
    }

    #[test]
    fn equal_area_projection() {
        // For equal-area: a'*b' = 1
        // Example: h=2.0, k=0.5, s=h*k*sin(π/2)=1.0
        let t = TissotParams::from_exact(2.0, 0.5, 1.0, PI / 2.0, 0.0);
        assert!(
            (t.areal_scale - 1.0).abs() < 1e-6,
            "areal={}",
            t.areal_scale
        );
        assert!(t.is_equal_area);
    }

    #[test]
    fn semi_major_ge_semi_minor() {
        let t = TissotParams::from_exact(1.5, 0.8, 1.2, PI / 2.0, 0.0);
        assert!(
            t.semi_major >= t.semi_minor,
            "a'={} b'={}",
            t.semi_major,
            t.semi_minor
        );
    }

    #[test]
    fn display_format() {
        let t = TissotParams::from_exact(1.0, 1.0, 1.0, PI / 2.0, 0.0);
        let s = format!("{}", t);
        assert!(s.contains("Tissot"));
    }

    #[test]
    fn eccentricity_zero_for_conformal() {
        let t = TissotParams::from_exact(1.0, 1.0, 1.0, PI / 2.0, 0.0);
        assert!(t.eccentricity().abs() < 1e-9);
    }

    #[test]
    fn from_factors_integration() {
        use crate::FactorsExact;
        // Spherical Mercator at equator: h=k=1, s=1, conv=0
        let fe = FactorsExact::from_jacobian(1.0, 0.0, 0.0, 1.0, 0.0, 1.0);
        let t = TissotParams::from_factors(&fe);
        assert!((t.semi_major - 1.0).abs() < 1e-9, "a'={}", t.semi_major);
        assert!((t.semi_minor - 1.0).abs() < 1e-9, "b'={}", t.semi_minor);
        assert!(t.is_conformal);
    }

    #[cfg(feature = "serde")]
    #[test]
    fn tissot_params_serde_round_trip() {
        use core::f64::consts::PI;
        let original = TissotParams::from_exact(1.5, 0.8, 1.2, PI / 2.0, 0.1);
        let json = serde_json::to_string(&original).expect("serialize");
        let restored: TissotParams = serde_json::from_str(&json).expect("deserialize");
        // Use approximate comparison for f64 fields (JSON may not preserve last ULP).
        let tol = 1e-10_f64;
        assert!((original.semi_major - restored.semi_major).abs() < tol);
        assert!((original.semi_minor - restored.semi_minor).abs() < tol);
        assert!((original.meridian_scale - restored.meridian_scale).abs() < tol);
        assert!((original.parallel_scale - restored.parallel_scale).abs() < tol);
        assert!((original.angular_distortion - restored.angular_distortion).abs() < tol);
        assert!((original.areal_scale - restored.areal_scale).abs() < tol);
        assert!((original.meridian_convergence - restored.meridian_convergence).abs() < tol);
        assert_eq!(original.is_conformal, restored.is_conformal);
        assert_eq!(original.is_equal_area, restored.is_equal_area);
    }
}