geometry_strategy/geographic/
direct_karney.rs

1//! Karney eighth-order solution of the direct geodesic problem.
2//!
3//! Mirrors `boost::geometry::formula::karney_direct` from
4//! `formulas/karney_direct.hpp:52-255`, using the order-8 coefficient tables
5//! from [`geometry_coords::series_expansion`].
6
7use geometry_cs::Spheroid;
8
9use super::direct::DirectResult;
10
11#[cfg(feature = "std")]
12use geometry_coords::series_expansion::{
13    coefficients_a3, coefficients_c1, coefficients_c1p, coefficients_c3, evaluate_a1,
14    sin_cos_series,
15};
16
17#[cfg(feature = "std")]
18use super::direct::normalize_longitude;
19#[cfg(feature = "std")]
20use super::spheroid_calc::SpheroidCalc;
21
22/// Karney's eighth-order direct geodesic formula.
23///
24/// Given a longitude/latitude, distance, and initial azimuth, computes the
25/// destination and final azimuth. Inputs and outputs are radians; distance
26/// uses the spheroid radius unit. Mirrors `formula::karney_direct<CT, ..., 8>`
27/// from `formulas/karney_direct.hpp:52-255`.
28#[derive(Debug, Clone, Copy)]
29pub struct KarneyDirect {
30    /// Reference ellipsoid.
31    pub spheroid: Spheroid,
32}
33
34impl KarneyDirect {
35    /// Eighth-order Karney direct on WGS84.
36    pub const WGS84: Self = Self {
37        spheroid: Spheroid::WGS84,
38    };
39
40    /// Solve the direct geodesic problem.
41    ///
42    /// Mirrors `karney_direct::apply` at
43    /// `formulas/karney_direct.hpp:76-253`. The Rust calculation stays in
44    /// radians instead of temporarily converting angles to degrees; both
45    /// forms evaluate the same trigonometric quantities, while the final
46    /// longitude is normalized identically.
47    #[cfg(feature = "std")]
48    #[inline]
49    #[must_use]
50    #[allow(
51        clippy::many_single_char_names,
52        clippy::similar_names,
53        clippy::float_cmp,
54        reason = "names and exact zero branches follow the cited Karney formula"
55    )]
56    pub fn apply(&self, lon1: f64, lat1: f64, distance: f64, azimuth12: f64) -> DirectResult {
57        let calc = SpheroidCalc::from(self.spheroid);
58        let b = calc.b;
59        let f = calc.f;
60        let one_minus_f = 1.0 - f;
61        let two_minus_f = 2.0 - f;
62        let n = f / two_minus_f;
63        let e2 = f * two_minus_f;
64        let ep2 = e2 / (one_minus_f * one_minus_f);
65
66        let sin_alpha1 = azimuth12.sin();
67        let sin_alpha1 = if sin_alpha1.abs() <= f64::EPSILON {
68            0.0
69        } else {
70            sin_alpha1
71        };
72        let cos_alpha1 = azimuth12.cos();
73        let cos_alpha1 = if cos_alpha1.abs() <= f64::EPSILON {
74            0.0
75        } else {
76            cos_alpha1
77        };
78        let mut sin_beta1 = lat1.sin() * one_minus_f;
79        let mut cos_beta1 = lat1.cos();
80        let beta_norm = sin_beta1.hypot(cos_beta1);
81        sin_beta1 /= beta_norm;
82        cos_beta1 = (cos_beta1 / beta_norm).max(0.0);
83
84        let sin_alpha0 = sin_alpha1 * cos_beta1;
85        let cos_alpha0 = cos_alpha1.hypot(sin_alpha1 * sin_beta1);
86        let k2 = cos_alpha0 * cos_alpha0 * ep2;
87        let epsilon = k2 / (2.0 * (1.0 + (1.0 + k2).sqrt()) + k2);
88        let expansion_a1 = evaluate_a1(epsilon);
89        let coefficients_c1 = coefficients_c1(epsilon);
90        let tau12 = distance / (b * (1.0 + expansion_a1));
91        let sin_tau12 = tau12.sin();
92        let cos_tau12 = tau12.cos();
93
94        let mut sin_sigma1 = sin_beta1;
95        let sin_omega1 = sin_alpha0 * sin_beta1;
96        // Boost evaluates these terms with `sin_cos_degrees`, which returns an
97        // exact zero for the equatorial due-east/west cases. The radian
98        // trigonometric functions leave a sub-epsilon residue instead.
99        let mut cos_sigma1 = if sin_beta1.abs() > f64::EPSILON || cos_alpha1.abs() > f64::EPSILON {
100            cos_beta1 * cos_alpha1
101        } else {
102            1.0
103        };
104        let cos_omega1 = cos_sigma1;
105        let sigma_norm = sin_sigma1.hypot(cos_sigma1);
106        sin_sigma1 /= sigma_norm;
107        cos_sigma1 /= sigma_norm;
108
109        let b11 = sin_cos_series(sin_sigma1, cos_sigma1, &coefficients_c1);
110        let sin_b11 = b11.sin();
111        let cos_b11 = b11.cos();
112        let sin_tau1 = sin_sigma1 * cos_b11 + cos_sigma1 * sin_b11;
113        let cos_tau1 = cos_sigma1 * cos_b11 - sin_sigma1 * sin_b11;
114        let coefficients_c1p = coefficients_c1p(epsilon);
115        let b12 = -sin_cos_series(
116            sin_tau1 * cos_tau12 + cos_tau1 * sin_tau12,
117            cos_tau1 * cos_tau12 - sin_tau1 * sin_tau12,
118            &coefficients_c1p,
119        );
120        let sigma12 = tau12 - (b12 - b11);
121        let sin_sigma12 = sigma12.sin();
122        let cos_sigma12 = sigma12.cos();
123        let sin_sigma2 = sin_sigma1 * cos_sigma12 + cos_sigma1 * sin_sigma12;
124        let cos_sigma2 = cos_sigma1 * cos_sigma12 - sin_sigma1 * sin_sigma12;
125
126        let sin_alpha2 = sin_alpha0;
127        let cos_alpha2 = cos_alpha0 * cos_sigma2;
128        let reverse_azimuth = sin_alpha2.atan2(cos_alpha2);
129
130        let sin_beta2 = cos_alpha0 * sin_sigma2;
131        let cos_beta2 = sin_alpha0.hypot(cos_alpha0 * cos_sigma2);
132        let lat2 = sin_beta2.atan2(one_minus_f * cos_beta2);
133
134        let sin_omega2 = sin_alpha0 * sin_sigma2;
135        let cos_omega2 = cos_sigma2;
136        let omega12 = (sin_omega2 * cos_omega1 - cos_omega2 * sin_omega1)
137            .atan2(cos_omega2 * cos_omega1 + sin_omega2 * sin_omega1);
138        let coefficients_a3 = coefficients_a3(n);
139        let a3 = coefficients_a3
140            .as_slice()
141            .iter()
142            .rev()
143            .fold(0.0, |value, coefficient| value * epsilon + coefficient);
144        let a3c = -f * sin_alpha0 * a3;
145        let coefficients_c3 = coefficients_c3(n, epsilon);
146        let b31 = sin_cos_series(sin_sigma1, cos_sigma1, &coefficients_c3);
147        let b32 = sin_cos_series(sin_sigma2, cos_sigma2, &coefficients_c3);
148        let lambda12 = omega12 + a3c * (sigma12 + b32 - b31);
149
150        DirectResult::solved(
151            lon1,
152            lat1,
153            azimuth12,
154            self.spheroid,
155            normalize_longitude(lon1 + lambda12),
156            lat2,
157            reverse_azimuth,
158        )
159    }
160}
161
162impl Default for KarneyDirect {
163    #[inline]
164    fn default() -> Self {
165        Self::WGS84
166    }
167}
geometry_strategy/geographic/direct_karney.rs

geometry_strategy/geographic/
direct_karney.rs
