geometry_strategy/geographic/
meridian.rs

1//! Meridian-arc formulas on a reference spheroid.
2//!
3//! Ports the coherent aggregate formed by Boost.Geometry's
4//! `meridian_direct.hpp`, `meridian_inverse.hpp`, `meridian_segment.hpp`, and
5//! `quarter_meridian.hpp`. Angles are radians and distances use the spheroid's
6//! radius unit.
7
8use geometry_cs::Spheroid;
9
10use super::direct::DirectResult;
11
12/// Classification of a spheroidal segment relative to a meridian.
13#[derive(Debug, Clone, Copy, PartialEq, Eq)]
14pub enum MeridianSegmentKind {
15    /// The endpoints do not define a meridian segment.
16    NonMeridian,
17    /// Both endpoints use the same meridian and the path avoids a pole.
18    NotCrossingPole,
19    /// The longitudes differ by π and the meridian path crosses a pole.
20    CrossingPole,
21}
22
23/// Distance result from recognizing a meridian endpoint pair.
24#[derive(Debug, Clone, Copy, PartialEq)]
25pub struct MeridianInverseResult {
26    /// Meridian distance, or zero when the pair is not meridional.
27    pub distance: f64,
28    /// Whether the endpoint pair was recognized as meridional.
29    pub meridian: bool,
30}
31
32/// Meridian formulas bound to a reference spheroid.
33#[derive(Debug, Clone, Copy, PartialEq)]
34pub struct Meridian {
35    /// Reference ellipsoid.
36    pub spheroid: Spheroid,
37}
38
39impl Meridian {
40    /// WGS84 meridian formulas.
41    pub const WGS84: Self = Self {
42        spheroid: Spheroid::WGS84,
43    };
44
45    /// Length from the equator to a pole.
46    ///
47    /// Mirrors the order-eight generalized series in
48    /// `formulas/quarter_meridian.hpp:54-76`.
49    #[must_use]
50    pub fn quarter_length(&self) -> f64 {
51        const COEFFICIENTS: [f64; 9] = [
52            1_073_741_824.0,
53            268_435_456.0,
54            16_777_216.0,
55            4_194_304.0,
56            1_638_400.0,
57            802_816.0,
58            451_584.0,
59            278_784.0,
60            184_041.0,
61        ];
62        let f = self.spheroid.flattening;
63        let n = f / (2.0 - f);
64        let ab4 = (self.spheroid.equatorial_radius + self.spheroid.polar_radius()) / 4.0;
65        let series = COEFFICIENTS[..8]
66            .iter()
67            .rev()
68            .fold(0.0, |value, coefficient| value * n * n + coefficient);
69        core::f64::consts::PI * ab4 * series / COEFFICIENTS[0]
70    }
71
72    /// Signed meridian arc from the equator to `latitude`.
73    ///
74    /// Uses Boost's maximum order-five expansion from
75    /// `formulas/meridian_inverse.hpp:116-178`.
76    #[cfg(feature = "std")]
77    #[must_use]
78    pub fn arc_length(&self, latitude: f64) -> f64 {
79        let f = self.spheroid.flattening;
80        let n = f / (2.0 - f);
81        let n2 = n * n;
82        let n3 = n2 * n;
83        let n4 = n2 * n2;
84        let n5 = n4 * n;
85        let scale = self.spheroid.equatorial_radius / (1.0 + n);
86        let c0 = 1.0 + 0.25 * n2;
87        let c2 = -1.5 * n + 0.1875 * n3;
88        let c4 = 0.9375 * n2 - 0.234_375 * n4;
89        let c6 = -0.729_166_667 * n3 + 0.227_864_583 * n5;
90        let c8 = 0.615_234_375 * n4;
91        let c10 = -0.541_406_25 * n5;
92        scale
93            * (c0 * latitude
94                + c2 * (2.0 * latitude).sin()
95                + c4 * (4.0 * latitude).sin()
96                + c6 * (6.0 * latitude).sin()
97                + c8 * (8.0 * latitude).sin()
98                + c10 * (10.0 * latitude).sin())
99    }
100
101    /// Latitude whose signed meridian arc is `distance`.
102    ///
103    /// Mirrors the order-four inverse series in
104    /// `formulas/meridian_direct.hpp:123-171`.
105    #[cfg(feature = "std")]
106    #[must_use]
107    pub fn latitude_at_arc(&self, distance: f64) -> f64 {
108        let f = self.spheroid.flattening;
109        let n = f / (2.0 - f);
110        let n2 = n * n;
111        let n3 = n2 * n;
112        let n4 = n2 * n2;
113        let mu = core::f64::consts::FRAC_PI_2 * distance / self.quarter_length();
114        let h2 = 1.5 * n - 0.843_75 * n3;
115        let h4 = 1.3125 * n2 - 1.718_75 * n4;
116        let h6 = 1.572_916_667 * n3;
117        let h8 = 2.142_578_125 * n4;
118        mu + h2 * (2.0 * mu).sin()
119            + h4 * (4.0 * mu).sin()
120            + h6 * (6.0 * mu).sin()
121            + h8 * (8.0 * mu).sin()
122    }
123
124    /// Classify the endpoint pair as a meridian segment.
125    #[cfg(feature = "std")]
126    #[must_use]
127    #[allow(
128        clippy::unused_self,
129        reason = "the instance method keeps the formula strategy API uniform with the spheroid-dependent methods"
130    )]
131    pub fn classify_segment(
132        &self,
133        longitude1: f64,
134        latitude1: f64,
135        longitude2: f64,
136        latitude2: f64,
137    ) -> MeridianSegmentKind {
138        let difference = normalize_longitude(longitude2 - longitude1);
139        if nearly_zero(difference)
140            || (nearly_equal(latitude2, core::f64::consts::FRAC_PI_2)
141                && nearly_equal(latitude1, -core::f64::consts::FRAC_PI_2))
142            || (nearly_equal(latitude1, core::f64::consts::FRAC_PI_2)
143                && nearly_equal(latitude2, -core::f64::consts::FRAC_PI_2))
144        {
145            MeridianSegmentKind::NotCrossingPole
146        } else if nearly_equal(difference.abs(), core::f64::consts::PI) {
147            MeridianSegmentKind::CrossingPole
148        } else {
149            MeridianSegmentKind::NonMeridian
150        }
151    }
152
153    /// Solve the inverse problem when the endpoints form a meridian.
154    #[cfg(feature = "std")]
155    #[must_use]
156    pub fn inverse(
157        &self,
158        longitude1: f64,
159        mut latitude1: f64,
160        longitude2: f64,
161        mut latitude2: f64,
162    ) -> MeridianInverseResult {
163        let kind = self.classify_segment(longitude1, latitude1, longitude2, latitude2);
164        if latitude1 > latitude2 {
165            core::mem::swap(&mut latitude1, &mut latitude2);
166        }
167        let distance = match kind {
168            MeridianSegmentKind::NonMeridian => 0.0,
169            MeridianSegmentKind::NotCrossingPole => {
170                (self.arc_length(latitude2) - self.arc_length(latitude1)).abs()
171            }
172            MeridianSegmentKind::CrossingPole => {
173                let latitude_sign = if latitude1 + latitude2 < 0.0 {
174                    -1.0
175                } else {
176                    1.0
177                };
178                (latitude_sign * 2.0 * self.arc_length(core::f64::consts::FRAC_PI_2)
179                    - self.arc_length(latitude1)
180                    - self.arc_length(latitude2))
181                .abs()
182            }
183        };
184        MeridianInverseResult {
185            distance,
186            meridian: kind != MeridianSegmentKind::NonMeridian,
187        }
188    }
189
190    /// Solve the direct geodesic problem along a meridian.
191    #[cfg(feature = "std")]
192    #[must_use]
193    pub fn direct(
194        &self,
195        longitude1: f64,
196        latitude1: f64,
197        distance: f64,
198        north: bool,
199    ) -> DirectResult {
200        let initial_arc = self.arc_length(latitude1);
201        let signed_distance = if north { distance } else { -distance };
202        let raw_latitude = self.latitude_at_arc(initial_arc + signed_distance);
203        let (lon2, lat2, reflected) = normalize_coordinates(longitude1, raw_latitude);
204        let final_north = north ^ reflected;
205        DirectResult::solved(
206            longitude1,
207            latitude1,
208            if north { 0.0 } else { core::f64::consts::PI },
209            self.spheroid,
210            lon2,
211            lat2,
212            if final_north {
213                0.0
214            } else {
215                core::f64::consts::PI
216            },
217        )
218    }
219}
220
221impl Default for Meridian {
222    fn default() -> Self {
223        Self::WGS84
224    }
225}
226
227#[cfg(feature = "std")]
228fn normalize_coordinates(longitude: f64, latitude: f64) -> (f64, f64, bool) {
229    let pi = core::f64::consts::PI;
230    let mut lat = (latitude + pi).rem_euclid(core::f64::consts::TAU) - pi;
231    let mut lon = longitude;
232    let reflected = if lat > core::f64::consts::FRAC_PI_2 {
233        lat = pi - lat;
234        lon += pi;
235        true
236    } else if lat < -core::f64::consts::FRAC_PI_2 {
237        lat = -pi - lat;
238        lon += pi;
239        true
240    } else {
241        false
242    };
243    (normalize_longitude(lon), lat, reflected)
244}
245
246#[cfg(feature = "std")]
247fn normalize_longitude(longitude: f64) -> f64 {
248    let pi = core::f64::consts::PI;
249    (longitude + pi).rem_euclid(core::f64::consts::TAU) - pi
250}
251
252#[cfg(feature = "std")]
253fn nearly_zero(value: f64) -> bool {
254    value.abs() <= 1e-12
255}
256
257#[cfg(feature = "std")]
258fn nearly_equal(first: f64, second: f64) -> bool {
259    (first - second).abs() <= 1e-12
260}
geometry_strategy/geographic/meridian.rs

geometry_strategy/geographic/
meridian.rs
