geometry_strategy/geographic/
direct_vincenty.rs

1//! Vincenty solution of the direct geodesic problem on a spheroid.
2//!
3//! Mirrors `boost::geometry::formula::vincenty_direct` from
4//! `formulas/vincenty_direct.hpp:43-178`.
5
6use geometry_cs::Spheroid;
7
8use super::direct::DirectResult;
9
10#[cfg(feature = "std")]
11use super::direct::normalize_longitude;
12#[cfg(feature = "std")]
13use super::spheroid_calc::SpheroidCalc;
14
15/// Vincenty's iterative direct geodesic formula.
16///
17/// Given a longitude/latitude, distance, and initial azimuth, computes the
18/// destination and final azimuth. Inputs and output angles are radians;
19/// distance uses the spheroid radius unit. Mirrors
20/// `formula::vincenty_direct<CT>` from
21/// `formulas/vincenty_direct.hpp:43-178`.
22#[derive(Debug, Clone, Copy)]
23pub struct VincentyDirect {
24    /// Reference ellipsoid.
25    pub spheroid: Spheroid,
26    /// Iteration limit, matching Boost's default of 1000.
27    pub max_iterations: u32,
28    /// Convergence threshold for auxiliary-sphere arc length.
29    pub tolerance: f64,
30}
31
32impl VincentyDirect {
33    /// Vincenty direct on WGS84 with Boost's iteration settings.
34    pub const WGS84: Self = Self {
35        spheroid: Spheroid::WGS84,
36        max_iterations: 1000,
37        tolerance: 1e-12,
38    };
39
40    /// Solve the direct geodesic problem.
41    ///
42    /// Mirrors `vincenty_direct::apply` from
43    /// `formulas/vincenty_direct.hpp:68-176`, including the iterative
44    /// `delta_sigma` correction and final longitude normalization.
45    #[cfg(feature = "std")]
46    #[inline]
47    #[must_use]
48    #[allow(
49        clippy::many_single_char_names,
50        clippy::similar_names,
51        reason = "names follow Vincenty's equations and the cited Boost implementation"
52    )]
53    pub fn apply(&self, lon1: f64, lat1: f64, distance: f64, azimuth12: f64) -> DirectResult {
54        let calc = SpheroidCalc::from(self.spheroid);
55        let radius_a = calc.a;
56        let radius_b = calc.b;
57        let flattening = calc.f;
58
59        let sin_azimuth12 = azimuth12.sin();
60        let cos_azimuth12 = azimuth12.cos();
61        let one_min_f = 1.0 - flattening;
62        let tan_u1 = one_min_f * lat1.tan();
63        let sigma1 = tan_u1.atan2(cos_azimuth12);
64        let u1 = tan_u1.atan();
65        let sin_u1 = u1.sin();
66        let cos_u1 = u1.cos();
67
68        let sin_alpha = cos_u1 * sin_azimuth12;
69        let sin_alpha_sqr = sin_alpha * sin_alpha;
70        let cos_alpha_sqr = 1.0 - sin_alpha_sqr;
71        let b_sqr = radius_b * radius_b;
72        let u_sqr = cos_alpha_sqr * (radius_a * radius_a - b_sqr) / b_sqr;
73        let a = 1.0
74            + (u_sqr / 16_384.0) * (4096.0 + u_sqr * (-768.0 + u_sqr * (320.0 - u_sqr * 175.0)));
75        let b = (u_sqr / 1024.0) * (256.0 + u_sqr * (-128.0 + u_sqr * (74.0 - u_sqr * 47.0)));
76
77        let s_div_ba = distance / (radius_b * a);
78        let mut sigma = s_div_ba;
79        let mut sin_sigma;
80        let mut cos_sigma;
81        let mut cos_2sigma_m;
82        let mut cos_2sigma_m_sqr;
83        let mut counter = 0;
84        loop {
85            let previous_sigma = sigma;
86            let two_sigma_m = 2.0 * sigma1 + sigma;
87            sin_sigma = sigma.sin();
88            cos_sigma = sigma.cos();
89            let sin_sigma_sqr = sin_sigma * sin_sigma;
90            cos_2sigma_m = two_sigma_m.cos();
91            cos_2sigma_m_sqr = cos_2sigma_m * cos_2sigma_m;
92            let delta_sigma = b
93                * sin_sigma
94                * (cos_2sigma_m
95                    + (b / 4.0)
96                        * (cos_sigma * (-1.0 + 2.0 * cos_2sigma_m_sqr)
97                            - (b / 6.0)
98                                * cos_2sigma_m
99                                * (-3.0 + 4.0 * sin_sigma_sqr)
100                                * (-3.0 + 4.0 * cos_2sigma_m_sqr)));
101            sigma = s_div_ba + delta_sigma;
102            counter += 1;
103            if (previous_sigma - sigma).abs() <= self.tolerance || counter >= self.max_iterations {
104                break;
105            }
106        }
107
108        let lat2 = (sin_u1 * cos_sigma + cos_u1 * sin_sigma * cos_azimuth12).atan2(
109            one_min_f
110                * (sin_alpha_sqr
111                    + (sin_u1 * sin_sigma - cos_u1 * cos_sigma * cos_azimuth12).powi(2))
112                .sqrt(),
113        );
114        let lambda = (sin_sigma * sin_azimuth12)
115            .atan2(cos_u1 * cos_sigma - sin_u1 * sin_sigma * cos_azimuth12);
116        let c =
117            (flattening / 16.0) * cos_alpha_sqr * (4.0 + flattening * (4.0 - 3.0 * cos_alpha_sqr));
118        let big_l = lambda
119            - (1.0 - c)
120                * flattening
121                * sin_alpha
122                * (sigma
123                    + c * sin_sigma
124                        * (cos_2sigma_m + c * cos_sigma * (-1.0 + 2.0 * cos_2sigma_m_sqr)));
125        let reverse_azimuth =
126            sin_alpha.atan2(-sin_u1 * sin_sigma + cos_u1 * cos_sigma * cos_azimuth12);
127
128        DirectResult::solved(
129            lon1,
130            lat1,
131            azimuth12,
132            self.spheroid,
133            normalize_longitude(lon1 + big_l),
134            lat2,
135            reverse_azimuth,
136        )
137    }
138}
139
140impl Default for VincentyDirect {
141    #[inline]
142    fn default() -> Self {
143        Self::WGS84
144    }
145}
geometry_strategy/geographic/direct_vincenty.rs

geometry_strategy/geographic/
direct_vincenty.rs
