geometry_strategy/geographic/
differential.rs

1//! Reduced length and geodesic scale between spheroidal endpoints.
2//!
3//! Mirrors `boost::geometry::formula::differential_quantities` with the
4//! maximum third-order flattening expansion exposed by the C++ header.
5
6use geometry_cs::Spheroid;
7
8/// Differential quantities attached to a direct or inverse geodesic result.
9#[derive(Debug, Clone, Copy, PartialEq)]
10pub struct DifferentialQuantities {
11    /// Reduced geodesic length in the spheroid's radius unit.
12    pub reduced_length: f64,
13    /// Dimensionless forward geodesic scale.
14    pub geodesic_scale: f64,
15}
16
17impl Default for DifferentialQuantities {
18    fn default() -> Self {
19        Self {
20            reduced_length: 0.0,
21            geodesic_scale: 1.0,
22        }
23    }
24}
25
26/// Calculate reduced length and geodesic scale for a solved geodesic.
27#[cfg(feature = "std")]
28#[must_use]
29#[allow(
30    clippy::too_many_arguments,
31    clippy::many_single_char_names,
32    clippy::similar_names,
33    reason = "the inputs and symbols mirror Boost's differential formula"
34)]
35pub fn differential_quantities(
36    longitude1: f64,
37    latitude1: f64,
38    longitude2: f64,
39    latitude2: f64,
40    azimuth: f64,
41    reverse_azimuth: f64,
42    spheroid: Spheroid,
43) -> DifferentialQuantities {
44    let longitude_difference = longitude2 - longitude1;
45    let sin_latitude1 = latitude1.sin();
46    let cos_latitude1 = latitude1.cos();
47    let sin_latitude2 = latitude2.sin();
48    let cos_latitude2 = latitude2.cos();
49    let one_minus_f = 1.0 - spheroid.flattening;
50    let mut sin_beta1 = one_minus_f * sin_latitude1;
51    let mut sin_beta2 = one_minus_f * sin_latitude2;
52
53    if sin_beta1.abs() <= 1e-15 && sin_beta2.abs() <= 1e-15 {
54        let sigma12 = longitude_difference / one_minus_f;
55        let azimuth_sign = if azimuth >= 0.0 { 1.0 } else { -1.0 };
56        return DifferentialQuantities {
57            reduced_length: azimuth_sign * sigma12.sin() * spheroid.polar_radius(),
58            geodesic_scale: sigma12.cos(),
59        };
60    }
61
62    let f = spheroid.flattening;
63    let e2 = f * (2.0 - f);
64    let ep2 = e2 / (one_minus_f * one_minus_f);
65    let sin_alpha1 = azimuth.sin();
66    let cos_alpha1 = azimuth.cos();
67    let cos_alpha2 = reverse_azimuth.cos();
68    let mut cos_beta1 = cos_latitude1;
69    let mut cos_beta2 = cos_latitude2;
70    normalize(&mut sin_beta1, &mut cos_beta1);
71    normalize(&mut sin_beta2, &mut cos_beta2);
72    let mut sin_sigma1 = sin_beta1;
73    let mut cos_sigma1 = cos_alpha1 * cos_beta1;
74    let mut sin_sigma2 = sin_beta2;
75    let mut cos_sigma2 = cos_alpha2 * cos_beta2;
76    normalize(&mut sin_sigma1, &mut cos_sigma1);
77    normalize(&mut sin_sigma2, &mut cos_sigma2);
78    let sin_alpha0 = sin_alpha1 * cos_beta1;
79    let cos_alpha0_squared = 1.0 - sin_alpha0 * sin_alpha0;
80    let j12 = j12_flattening(
81        sin_sigma1,
82        cos_sigma1,
83        sin_sigma2,
84        cos_sigma2,
85        cos_alpha0_squared,
86        f,
87    );
88    let dn1 = (1.0 + ep2 * sin_beta1 * sin_beta1).sqrt();
89    let dn2 = (1.0 + ep2 * sin_beta2 * sin_beta2).sqrt();
90    let reduced_length = spheroid.polar_radius()
91        * (dn2 * cos_sigma1 * sin_sigma2
92            - dn1 * sin_sigma1 * cos_sigma2
93            - cos_sigma1 * cos_sigma2 * j12);
94    let cos_sigma12 = cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2;
95    let t = ep2 * (cos_beta1 - cos_beta2) * (cos_beta1 + cos_beta2) / (dn1 + dn2);
96    let geodesic_scale = cos_sigma12 + (t * sin_sigma2 - cos_sigma2 * j12) * sin_sigma1 / dn1;
97    DifferentialQuantities {
98        reduced_length,
99        geodesic_scale,
100    }
101}
102
103#[cfg(feature = "std")]
104#[allow(
105    clippy::similar_names,
106    reason = "sigma-indexed symbols mirror Boost's differential formula"
107)]
108fn j12_flattening(
109    sin_sigma1: f64,
110    cos_sigma1: f64,
111    sin_sigma2: f64,
112    cos_sigma2: f64,
113    cos_alpha0_squared: f64,
114    flattening: f64,
115) -> f64 {
116    let sigma12 = (cos_sigma1 * sin_sigma2 - sin_sigma1 * cos_sigma2)
117        .atan2(cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2);
118    let sin_2sigma1 = 2.0 * cos_sigma1 * sin_sigma1;
119    let sin_2sigma2 = 2.0 * cos_sigma2 * sin_sigma2;
120    let sin_2sigma12 = sin_2sigma2 - sin_2sigma1;
121    let l1 = sigma12 - sin_2sigma12 / 2.0;
122    let sin_4sigma1 = 2.0 * sin_2sigma1 * (cos_sigma1 * cos_sigma1 - sin_sigma1 * sin_sigma1);
123    let sin_4sigma2 = 2.0 * sin_2sigma2 * (cos_sigma2 * cos_sigma2 - sin_sigma2 * sin_sigma2);
124    let sin_4sigma12 = sin_4sigma2 - sin_4sigma1;
125    let l2 = -(cos_alpha0_squared * sin_4sigma12
126        + (-8.0 * cos_alpha0_squared + 12.0) * sin_2sigma12
127        + (12.0 * cos_alpha0_squared - 24.0) * sigma12)
128        / 16.0;
129    let cos_alpha0_fourth = cos_alpha0_squared * cos_alpha0_squared;
130    let sin_2sigma1_cubed = sin_2sigma1 * sin_2sigma1 * sin_2sigma1;
131    let sin_2sigma2_cubed = sin_2sigma2 * sin_2sigma2 * sin_2sigma2;
132    let l3 = ((9.0 * cos_alpha0_fourth - 12.0 * cos_alpha0_squared) * sin_4sigma12
133        + 4.0 * cos_alpha0_fourth * (sin_2sigma2_cubed - sin_2sigma1_cubed)
134        + (-48.0 * cos_alpha0_fourth + 96.0 * cos_alpha0_squared - 64.0) * sin_2sigma12
135        + (60.0 * cos_alpha0_fourth - 144.0 * cos_alpha0_squared + 128.0) * sigma12)
136        / 64.0;
137    cos_alpha0_squared * flattening * (l1 + flattening * (l2 + flattening * l3))
138}
139
140#[cfg(feature = "std")]
141fn normalize(x: &mut f64, y: &mut f64) {
142    let length = x.hypot(*y);
143    *x /= length;
144    *y /= length;
145}
geometry_strategy/geographic/differential.rs

geometry_strategy/geographic/
differential.rs
