oxiproj-core 0.1.2

//! Small-angle helper trig functions ported from PROJ `src/aasincos.cpp` and `src/adjlon.cpp`.
//!
//! These helpers guard the standard transcendental functions against
//! out-of-domain arguments that arise from floating-point round-off in the
//! projection forward/inverse formulas. PROJ records a context errno when an
//! argument is meaningfully out of domain; this port has no projection context
//! at this layer, so it follows PROJ's numerical behaviour (clamping /
//! returning a sentinel) without the side-effecting error report.

use crate::consts::{M_HALFPI, M_PI, M_TWOPI};

/// Tolerance above unity within which an out-of-domain `asin`/`acos` argument
/// is treated as merely the result of round-off rather than a real error.
/// Mirrors `ONE_TOL` in PROJ `src/aasincos.cpp`.
const ONE_TOL: f64 = 1.00000000000001;

/// Tolerance below which both `aatan2` arguments are considered zero.
/// Mirrors `ATOL` in PROJ `src/aasincos.cpp`.
const ATOL: f64 = 1e-50;

/// Ported from src/aasincos.cpp (`aasin`).
///
/// Faithful note: PROJ records a context errno when the argument magnitude
/// exceeds `ONE_TOL`; this port has no context, so it simply clamps the
/// result to +/- pi/2 and otherwise returns `v.asin()`.
pub fn aasin(v: f64) -> f64 {
    let av = v.abs();
    if av >= 1.0 {
        // PROJ sets a context errno here when `av > ONE_TOL` (i.e. the value is
        // genuinely out of domain rather than a round-off artefact). Without a
        // context we only observe the condition and then clamp below.
        let _out_of_domain = av > ONE_TOL;
        if v < 0.0 {
            -M_HALFPI
        } else {
            M_HALFPI
        }
    } else {
        v.asin()
    }
}

/// Ported from src/aasincos.cpp (`aacos`).
pub fn aacos(v: f64) -> f64 {
    if v.abs() >= 1.0 {
        if v < 0.0 {
            M_PI
        } else {
            0.0
        }
    } else {
        v.acos()
    }
}

/// Ported from src/aasincos.cpp (`asqrt`).
pub fn asqrt(v: f64) -> f64 {
    if v <= 0.0 {
        0.0
    } else {
        v.sqrt()
    }
}

/// Ported from src/aasincos.cpp (`aatan2`).
pub fn aatan2(n: f64, d: f64) -> f64 {
    if n.abs() < ATOL && d.abs() < ATOL {
        0.0
    } else {
        n.atan2(d)
    }
}

/// Ported from src/adjlon.cpp (`adjlon`). Reduces a longitude to (-pi, pi].
pub fn adjlon(mut longitude: f64) -> f64 {
    if longitude.abs() < M_PI + 1e-12 {
        return longitude;
    }
    longitude += M_PI;
    longitude -= M_TWOPI * (longitude / M_TWOPI).floor();
    longitude -= M_PI;
    longitude
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::consts::{M_HALFPI, M_PI, M_TWOPI};

    #[test]
    fn adjlon_reduces_above_range() {
        assert!((adjlon(M_TWOPI + 1.0) - 1.0).abs() < 1e-9);
    }

    #[test]
    fn adjlon_reduces_below_range() {
        assert!((adjlon(-M_TWOPI - 1.0) - (-1.0)).abs() < 1e-9);
    }

    #[test]
    fn adjlon_passes_small_angle_unchanged() {
        assert!(adjlon(0.5) == 0.5);
    }

    #[test]
    fn aasin_clamps_positive_out_of_domain() {
        assert!((aasin(2.0) - M_HALFPI).abs() < 1e-12);
    }

    #[test]
    fn aasin_clamps_negative_out_of_domain() {
        assert!((aasin(-2.0) - (-M_HALFPI)).abs() < 1e-12);
    }

    #[test]
    fn aasin_in_domain_matches_std() {
        assert!((aasin(0.5) - 0.5_f64.asin()).abs() < 1e-12);
    }

    #[test]
    fn aacos_clamps_positive_out_of_domain() {
        assert!(aacos(2.0) == 0.0);
    }

    #[test]
    fn aacos_clamps_negative_out_of_domain() {
        assert!((aacos(-2.0) - M_PI).abs() < 1e-12);
    }

    #[test]
    fn aacos_in_domain_matches_std() {
        assert!((aacos(0.5) - 0.5_f64.acos()).abs() < 1e-12);
    }

    #[test]
    fn asqrt_negative_is_zero() {
        assert!(asqrt(-1.0) == 0.0);
    }

    #[test]
    fn asqrt_positive_matches_std() {
        assert!(asqrt(4.0) == 2.0);
    }

    #[test]
    fn aatan2_both_zero_is_zero() {
        assert!(aatan2(0.0, 0.0) == 0.0);
    }

    #[test]
    fn aatan2_nonzero_matches_std() {
        assert!((aatan2(1.0, 1.0) - 1.0_f64.atan2(1.0)).abs() < 1e-12);
    }
}