jord 0.16.0

Geographical Position Calculations (Ellipsoidal and Spherical models)
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

use crate::{numbers::eq_zero, Vec3};

/// epsilon below which expensive side is called.
const TRIAGE_SIDE_EPS: f64 = 10.0 * f64::EPSILON;

/// Computes the signed angle in radians between the given vectors.
///
/// - if vn is `None; the angle is always in [0..PI],
/// - otherwise, the angle is positive if v1 is clockzise looking along vn,
/// - and negative if anti-clockwise looking along vn
pub(crate) fn angle_radians_between(v1: Vec3, v2: Vec3, vn: Option<Vec3>) -> f64 {
    let p1xp2 = v1.cross_prod(v2);
    let norm = p1xp2.norm();
    let sin_o = match vn {
        None => norm,
        Some(v) => {
            let d = p1xp2.dot_prod(v);
            if d >= 0.0 {
                norm
            } else {
                -norm
            }
        }
    };
    let cos_o = v1.dot_prod(v2);
    sin_o.atan2(cos_o)
}

/// Easting at given *n*-vector.
pub(crate) fn easting(v: Vec3) -> Vec3 {
    // north pole = (0, 0, 1), south pole = (0, 0, -1)
    // if v.z() == 1 or -1, then v is either north or south pole, easting is therefore (0, 1, 0)
    if v.z().abs() == 1.0 {
        return Vec3::UNIT_Y;
    }
    // cross product: (0, 0, 1) x v
    Vec3::new_unit(-v.y(), v.x(), 0.0)
}

/// Determines whether v0 if right (negative integer), left (positive integer) or on the great circle (0), from v1 to v2.
pub(crate) fn side(v0: Vec3, v1: Vec3, v2: Vec3) -> i8 {
    let triage_side = v0.dot_prod(v1.cross_prod(v2));
    // The side of v0 w.r.t. (v1, v2) is given by the triple scalar product (v0 . (v1 x v2))
    // However the cross product of v1 and v2 becomes unstable if v1 and v2 are nearly parallel (coincidental or antipodal).
    // If the result if too close to 0 (using 10 * f64::EPSILON), then call the more expensive function `orthogonal_to`.
    if triage_side.abs() <= TRIAGE_SIDE_EPS {
        let side = v0.dot_prod(v1.orthogonal_to(v2));
        if eq_zero(side) {
            return 0;
        }
        if side > 0.0 {
            1
        } else {
            -1
        }
    } else {
        // either right or left.
        if triage_side > 0.0 {
            1
        } else {
            -1
        }
    }
}

#[cfg(test)]
mod tests {

    // angle_radians_between

    use std::f64::consts::PI;

    use crate::{spherical::base::angle_radians_between, Vec3};

    #[test]
    fn angle_radians_between_signed() {
        let z = Vec3::new(0.0, 0.0, -1.0);
        assert_eq!(
            angle_radians_between(Vec3::new(0.0, 1.0, 0.0), Vec3::new(1.0, 0.0, 0.0), Some(z)),
            PI / 2.0
        );
        assert_eq!(
            angle_radians_between(Vec3::new(1.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0), Some(z)),
            -PI / 2.0
        );
        assert_eq!(
            angle_radians_between(Vec3::new(1.0, 0.0, 0.0), Vec3::new(1.0, 1.0, 0.0), Some(z)),
            -PI / 4.0
        );
        assert_eq!(
            angle_radians_between(Vec3::new(1.0, 1.0, 0.0), Vec3::new(1.0, 0.0, 0.0), Some(z)),
            PI / 4.0,
        );
    }

    #[test]
    fn angle_radians_between_unsigned() {
        assert_eq!(
            angle_radians_between(Vec3::new(0.0, 1.0, 0.0), Vec3::new(1.0, 0.0, 0.0), None),
            PI / 2.0
        );
        assert_eq!(
            angle_radians_between(Vec3::new(1.0, 0.0, 0.0), Vec3::new(1.0, 1.0, 0.0), None),
            PI / 4.0
        );
    }
}