geometry_strategy/spherical/
azimuth.rs

1//! Great-circle bearing (azimuth) on a sphere.
2//!
3//! Mirrors `boost::geometry::strategy::azimuth::spherical` from
4//! `boost/geometry/strategies/spherical/azimuth.hpp`, whose kernel is
5//! `formula::spherical_azimuth` in
6//! `boost/geometry/formulas/spherical.hpp:127-168`:
7//!
8//! ```text
9//! θ = atan2( sin(Δλ)·cos φ₂,
10//!            cos φ₁·sin φ₂ − sin φ₁·cos φ₂·cos(Δλ) )
11//! ```
12//!
13//! where φ₁, λ₁, φ₂, λ₂ are the input lat/lon in radians.
14//!
15//! # Convention
16//!
17//! `0` = north, increasing **clockwise** — the geodetic bearing
18//! convention. This is *different* from
19//! [`CartesianAzimuth`](crate::azimuth::CartesianAzimuth) (`0` = east,
20//! counter-clockwise). `azimuth(p, p)` returns `0`, mirroring the
21//! `math::equals(dlon, 0)`-style short-circuits Boost applies for
22//! coincident points.
23
24use geometry_cs::{CoordinateSystem, SphericalFamily};
25use geometry_tag::SameAs;
26use geometry_trait::Point;
27
28use crate::azimuth::AzimuthStrategy;
29
30#[cfg(feature = "std")]
31use crate::normalise::{HasAngularUnits, lonlat_radians};
32
33/// Great-circle initial bearing on a sphere, in radians (`0` = north,
34/// clockwise).
35///
36/// Mirrors `boost::geometry::strategy::azimuth::spherical`
37/// (`strategies/spherical/azimuth.hpp`). Radius-independent — a bearing
38/// is an angle, so no sphere radius is carried.
39#[derive(Debug, Default, Clone, Copy)]
40pub struct SphericalAzimuth;
41
42#[cfg(feature = "std")]
43impl<P1, P2> AzimuthStrategy<P1, P2> for SphericalAzimuth
44where
45    P1: Point<Scalar = f64>,
46    P2: Point<Scalar = f64>,
47    P1::Cs: HasAngularUnits,
48    P2::Cs: HasAngularUnits,
49    <P1::Cs as CoordinateSystem>::Family: SameAs<SphericalFamily>,
50    <P2::Cs as CoordinateSystem>::Family: SameAs<SphericalFamily>,
51{
52    type Out = f64;
53
54    #[inline]
55    fn azimuth(&self, p1: &P1, p2: &P2) -> f64 {
56        let (lon1, lat1) = lonlat_radians(p1);
57        let (lon2, lat2) = lonlat_radians(p2);
58        let dlon = lon2 - lon1;
59        // Mirrors `formula::spherical_azimuth` at
60        // `formulas/spherical.hpp:156-158`.
61        let y = dlon.sin() * lat2.cos();
62        let x = lat1.cos() * lat2.sin() - lat1.sin() * lat2.cos() * dlon.cos();
63        y.atan2(x)
64    }
65}
66
67#[cfg(all(test, feature = "std"))]
68mod tests {
69    //! Reference values from
70    //! `boost/geometry/test/algorithms/azimuth.cpp:49-56` (`test_sph`),
71    //! given in degrees; asserted here in radians.
72    #![allow(
73        clippy::float_cmp,
74        reason = "azimuths are compared with an explicit absolute tolerance, not `==`"
75    )]
76    #![allow(
77        clippy::excessive_precision,
78        reason = "reference constants are copied verbatim from Boost's azimuth.cpp; the trailing digits document the source value even past f64 precision."
79    )]
80
81    use super::SphericalAzimuth;
82    use crate::azimuth::AzimuthStrategy;
83    use geometry_adapt::{Adapt, WithCs};
84    use geometry_cs::{Degree, Spherical};
85
86    type Sp = WithCs<Adapt<[f64; 2]>, Spherical<Degree>>;
87
88    #[inline]
89    fn sp(lon: f64, lat: f64) -> Sp {
90        WithCs::new(Adapt([lon, lat]))
91    }
92
93    /// `azimuth.cpp:51` — `(0,0) → (0,0)` returns `0`.
94    #[test]
95    fn coincident_is_zero() {
96        let got = SphericalAzimuth.azimuth(&sp(0., 0.), &sp(0., 0.));
97        assert!(got.abs() < 1e-12);
98    }
99
100    /// `azimuth.cpp:52` — `(0,0) → (1,1) = 44.995636455°`.
101    #[test]
102    fn diagonal() {
103        let got = SphericalAzimuth.azimuth(&sp(0., 0.), &sp(1., 1.));
104        let expected = 44.995_636_455_344_851_f64.to_radians();
105        assert!((got - expected).abs() < 1e-9, "got {got}");
106    }
107
108    /// `azimuth.cpp:53` — `(0,0) → (100,1) = 88.984576593°`.
109    #[test]
110    fn shallow_east() {
111        let got = SphericalAzimuth.azimuth(&sp(0., 0.), &sp(100., 1.));
112        let expected = 88.984_576_593_576_293_f64.to_radians();
113        assert!((got - expected).abs() < 1e-9, "got {got}");
114    }
115
116    /// `azimuth.cpp:54` — `(0,0) → (-1,1) = -44.995636455°`.
117    #[test]
118    fn diagonal_west() {
119        let got = SphericalAzimuth.azimuth(&sp(0., 0.), &sp(-1., 1.));
120        let expected = (-44.995_636_455_344_851_f64).to_radians();
121        assert!((got - expected).abs() < 1e-9, "got {got}");
122    }
123}
geometry_strategy/spherical/azimuth.rs

geometry_strategy/spherical/
azimuth.rs
