comnum_math/
geometry.rs

1//! Vectors of the three-dimensional Euclidean space, as `[f64; 3]`: the scalar product,
2//! the difference, the norm and the direction of a vector, and the unit vector of
3//! spherical angles.
4//!
5//! The geometry of the antenna arrays of the CDL channel and of the indoor optical
6//! channel (positions, directions, orientations) is written with these.
7//!
8//! # Examples
9//! ```
10//! use comnum_math::geometry::{difference, dot, norm, unit};
11//! let (a, b) = ([1.0, 2.0, 2.0], [1.0, 0.0, 0.0]);
12//! assert_eq!(dot(a, b), 1.0);
13//! assert_eq!(norm(a), 3.0);
14//! assert_eq!(difference(a, b), [0.0, 2.0, 2.0]);
15//! assert_eq!(unit([0.0, 0.0, -2.0]), Some([0.0, 0.0, -1.0]));
16//! ```
17//!
18//! # Spherical angles
19//!
20//! [`from_spherical`] takes the zenith (polar) angle $\theta$ from $+z$ and the azimuth
21//! $\phi$ from $+x$ towards $+y$, in radians (ISO 80000-2; the angles of 3GPP TR 38.901,
22//! (7.5-23)):
23//!
24//! $$
25//! \hat{r}(\theta, \phi) = (\sin\theta \cos\phi,\ \sin\theta \sin\phi,\ \cos\theta)^{T}.
26//! $$
27//!
28//! An elevation $\epsilon$ above the $xy$ plane is the zenith $\theta = \pi/2 - \epsilon$.
29
30/// A vector of the three-dimensional space.
31pub type Vec3 = [f64; 3];
32
33/// The scalar product $a \cdot b = a_{0} b_{0} + a_{1} b_{1} + a_{2} b_{2}$.
34#[inline]
35pub fn dot(a: Vec3, b: Vec3) -> f64 {
36    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
37}
38
39/// The difference $a - b$.
40#[inline]
41pub fn difference(a: Vec3, b: Vec3) -> Vec3 {
42    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
43}
44
45/// The Euclidean norm $\lVert a \rVert = \sqrt{a \cdot a}$.
46#[inline]
47pub fn norm(a: Vec3) -> f64 {
48    dot(a, a).sqrt()
49}
50
51/// The direction $a / \lVert a \rVert$ of `a`; `None` for a zero or non-finite norm.
52pub fn unit(a: Vec3) -> Option<Vec3> {
53    let n = norm(a);
54    (n > 0.0 && n.is_finite()).then(|| [a[0] / n, a[1] / n, a[2] / n])
55}
56
57/// The unit vector $\hat{r}(\theta, \phi) = (\sin\theta \cos\phi, \sin\theta \sin\phi,
58/// \cos\theta)$ of the zenith `zenith_rad` (from $+z$) and the azimuth `azimuth_rad` (from
59/// $+x$ towards $+y$), see the module.
60///
61/// # Examples
62/// ```
63/// use comnum_math::geometry::from_spherical;
64/// use std::f64::consts::FRAC_PI_2;
65/// assert_eq!(from_spherical(0.0, 1.0), [0.0, 0.0, 1.0]);
66/// let y = from_spherical(FRAC_PI_2, FRAC_PI_2);
67/// assert!((y[1] - 1.0).abs() < 1e-15 && y[0].abs() < 1e-15 && y[2].abs() < 1e-15);
68/// ```
69#[inline]
70pub fn from_spherical(zenith_rad: f64, azimuth_rad: f64) -> Vec3 {
71    let (st, ct) = zenith_rad.sin_cos();
72    let (sp, cp) = azimuth_rad.sin_cos();
73    [st * cp, st * sp, ct]
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79
80    #[test]
81    fn a_spherical_direction_has_unit_norm() {
82        for (t, p) in [(0.3, -2.0), (1.2, 0.7), (3.0, 3.1)] {
83            let r = from_spherical(t, p);
84            assert!((norm(r) - 1.0).abs() < 1e-15);
85            assert!((r[2] - f64::cos(t)).abs() < 1e-15 && (r[1].atan2(r[0]) - p).abs() < 1e-12);
86        }
87    }
88
89    #[test]
90    fn a_direction_has_unit_norm_and_a_zero_vector_has_none() {
91        let u = unit([3.0, -4.0, 12.0]).unwrap();
92        assert!((norm(u) - 1.0).abs() < 1e-15);
93        assert_eq!(u, [3.0 / 13.0, -4.0 / 13.0, 12.0 / 13.0]);
94        assert_eq!(unit([0.0; 3]), None);
95        assert_eq!(unit([f64::INFINITY, 0.0, 0.0]), None);
96        assert_eq!(unit([f64::NAN, 0.0, 0.0]), None);
97    }
98}
comnum_math/geometry.rs

comnum_math/
geometry.rs
