oxiproj-core 0.1.2

Foundation types for OxiProj: coordinates, errors, ellipsoids, datums, and units.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

//! `pj_tsfn` ported from PROJ `src/tsfn.cpp`.

use crate::scalar::Scalar;

/// Ported from src/tsfn.cpp (`pj_tsfn`).
///
/// Computes the small `t` function used by conformal projections,
/// `exp(e*atanh(e*sin(phi))) * (cos(phi)/(1+sin(phi)))`, with a numerically
/// stable branch selected on the sign of `sin(phi)`.
pub fn pj_tsfn(phi: f64, sinphi: f64, e: f64) -> f64 {
    let cosphi = phi.cos();
    (e * (e * sinphi).atanh()).exp()
        * if sinphi > 0.0 {
            cosphi / (1.0 + sinphi)
        } else {
            (1.0 - sinphi) / cosphi
        }
}

/// Generic version of [`pj_tsfn`] for automatic-differentiation support.
///
/// Branch selection (`sinphi > 0`) is made on the primal value only, so
/// differentiation through branch boundaries is handled consistently.
/// `atanh` is not in the [`Scalar`] trait; it is expressed via the identity
/// `atanh(x) = 0.5 * ln((1+x)/(1-x))`.
pub fn pj_tsfn_g<S: Scalar>(phi: S, sinphi: S, e: f64) -> S {
    let e_s = S::from_f64(e);
    let cosphi = phi.cos();
    let e_sinphi = e_s * sinphi;
    let one = S::one();
    let half = S::from_f64(0.5);
    // atanh(e*sinphi) = 0.5 * ln((1 + e*sinphi) / (1 - e*sinphi))
    let atanh_e_sinphi = half * ((one + e_sinphi) / (one - e_sinphi)).ln();
    let exp_part = (e_s * atanh_e_sinphi).exp();
    // Branch on primal value only (same numerical branch as `pj_tsfn`).
    if sinphi.to_f64() > 0.0 {
        exp_part * (cosphi / (one + sinphi))
    } else {
        exp_part * ((one - sinphi) / cosphi)
    }
}

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

    #[test]
    fn tsfn_sphere_matches_closed_form() {
        let phi = 0.5_f64;
        let t = pj_tsfn(phi, phi.sin(), 0.0);
        // exp(0) == 1, so only the trig branch remains.
        let expected = phi.cos() / (1.0 + phi.sin());
        assert!((t - expected).abs() < 1e-12);
    }

    #[test]
    fn tsfn_wgs84_is_positive() {
        let phi = 0.6_f64;
        let e = 0.0818191908426;
        let t = pj_tsfn(phi, phi.sin(), e);
        assert!(t > 0.0);
    }
}