geographiclib_rs/
geodesic_line.rs

1#![allow(non_snake_case)]
2
3use crate::geodesic::{self, GEODESIC_ORDER};
4use crate::geodesic_capability as caps;
5use crate::geomath;
6use std::collections::HashMap;
7
8/// A geodesic line.
9///
10/// Facilitates the determination of a series of points on a single geodesic.
11#[derive(Copy, Clone, PartialEq, PartialOrd, Debug)]
12pub struct GeodesicLine {
13    tiny_: f64, // This should be moved to consts
14    _A1m1: f64,
15    _A2m1: f64,
16    _A3c: f64,
17    _A4: f64,
18    _B11: f64,
19    _B21: f64,
20    _B31: f64,
21    _B41: f64,
22    _C1a: [f64; GEODESIC_ORDER + 1],
23    _C1pa: [f64; GEODESIC_ORDER + 1],
24    _C2a: [f64; GEODESIC_ORDER + 1],
25    _C3a: [f64; GEODESIC_ORDER],
26    _C4a: [f64; GEODESIC_ORDER],
27    _b: f64,
28    _c2: f64,
29    _calp0: f64,
30    _csig1: f64,
31    _comg1: f64,
32    _ctau1: f64,
33    _dn1: f64,
34    _f1: f64,
35    _k2: f64,
36    _salp0: f64,
37    _somg1: f64,
38    _ssig1: f64,
39    _stau1: f64,
40    _a13: f64,
41    _a: f64,
42    azi1: f64,
43    calp1: f64,
44    caps: u64,
45    f: f64,
46    lat1: f64,
47    lon1: f64,
48    _s13: f64,
49    salp1: f64,
50}
51
52impl GeodesicLine {
53    /// Create a new geodesic line from an initial point and azimuth.
54    ///
55    /// # Arguments
56    ///   - geod - geodesic parameters
57    ///   - lat1 - initial latitude (degrees)
58    ///   - lon1 - initial longitude (degrees)
59    ///   - azi1 - initial azimuth (degrees)
60    ///   - caps - bitor'ed combination of [capabilities](caps); defaults to `STANDARD | DISTANCE_IN`
61    ///     - `caps |= LATITUDE` for the latitude lat2; this is added automatically
62    ///     - `caps |= LONGITUDE` for the longitude lon2
63    ///     - `caps |= AZIMUTH` for the azimuth azi2; this is added automatically
64    ///     - `caps |= DISTANCE` for the distance s12
65    ///     - `caps |= REDUCEDLENGTH` for the reduced length m12
66    ///     - `caps |= GEODESICSCALE` for the geodesic scales M12 and M21
67    ///     - `caps |= AREA` for the area S12
68    ///     - `caps |= DISTANCE_IN` permits the length of the geodesic to be given in terms of s12;
69    ///       without this capability the length can only be specified in terms of arc length
70    ///     - `caps |= ALL` for all of the above
71    ///   - salp1 - sine of initial azimuth
72    ///   - calp1 - cosine of initial azimuth
73    pub fn new(
74        geod: &geodesic::Geodesic,
75        lat1: f64,
76        lon1: f64,
77        azi1: f64,
78        caps: Option<u64>,
79        salp1: Option<f64>,
80        calp1: Option<f64>,
81    ) -> Self {
82        let caps = match caps {
83            None => caps::STANDARD | caps::DISTANCE_IN,
84            Some(caps) => caps,
85        };
86        let salp1 = match salp1 {
87            None => f64::NAN,
88            Some(salp1) => salp1,
89        };
90        let calp1 = match calp1 {
91            None => f64::NAN,
92            Some(calp1) => calp1,
93        };
94
95        // This was taken from geodesic, putting it here for convenience
96        let tiny_ = geomath::get_min_val().sqrt();
97
98        let _a = geod.a;
99        let f = geod.f;
100        let _b = geod._b;
101        let _c2 = geod._c2;
102        let _f1 = geod._f1;
103        let caps = caps | caps::LATITUDE | caps::AZIMUTH | caps::LONG_UNROLL;
104        let (azi1, salp1, calp1) = if salp1.is_nan() || calp1.is_nan() {
105            let azi1 = geomath::ang_normalize(azi1);
106            let (salp1, calp1) = geomath::sincosd(geomath::ang_round(azi1));
107            (azi1, salp1, calp1)
108        } else {
109            (azi1, salp1, calp1)
110        };
111        let lat1 = geomath::lat_fix(lat1);
112
113        let (mut sbet1, mut cbet1) = geomath::sincosd(geomath::ang_round(lat1));
114        sbet1 *= _f1;
115        geomath::norm(&mut sbet1, &mut cbet1);
116        cbet1 = tiny_.max(cbet1);
117        let _dn1 = (1.0 + geod._ep2 * geomath::sq(sbet1)).sqrt();
118        let _salp0 = salp1 * cbet1;
119        let _calp0 = calp1.hypot(salp1 * sbet1);
120        let mut _ssig1 = sbet1;
121        let _somg1 = _salp0 * sbet1;
122        let mut _csig1 = if sbet1 != 0.0 || calp1 != 0.0 {
123            cbet1 * calp1
124        } else {
125            1.0
126        };
127        let _comg1 = _csig1;
128        geomath::norm(&mut _ssig1, &mut _csig1);
129        let _k2 = geomath::sq(_calp0) * geod._ep2;
130        let eps = _k2 / (2.0 * (1.0 + (1.0 + _k2).sqrt()) + _k2);
131
132        let mut _A1m1 = 0.0;
133        let mut _C1a: [f64; GEODESIC_ORDER + 1] = [0.0; GEODESIC_ORDER + 1];
134        let mut _B11 = 0.0;
135        let mut _stau1 = 0.0;
136        let mut _ctau1 = 0.0;
137        if caps & caps::CAP_C1 != 0 {
138            _A1m1 = geomath::_A1m1f(eps, geod.GEODESIC_ORDER);
139            geomath::_C1f(eps, &mut _C1a, geod.GEODESIC_ORDER);
140            _B11 = geomath::sin_cos_series(true, _ssig1, _csig1, &_C1a);
141            let s = _B11.sin();
142            let c = _B11.cos();
143            _stau1 = _ssig1 * c + _csig1 * s;
144            _ctau1 = _csig1 * c - _ssig1 * s;
145        }
146
147        let mut _C1pa: [f64; GEODESIC_ORDER + 1] = [0.0; GEODESIC_ORDER + 1];
148        if caps & caps::CAP_C1p != 0 {
149            geomath::_C1pf(eps, &mut _C1pa, geod.GEODESIC_ORDER);
150        }
151
152        let mut _A2m1 = 0.0;
153        let mut _C2a: [f64; GEODESIC_ORDER + 1] = [0.0; GEODESIC_ORDER + 1];
154        let mut _B21 = 0.0;
155        if caps & caps::CAP_C2 != 0 {
156            _A2m1 = geomath::_A2m1f(eps, geod.GEODESIC_ORDER);
157            geomath::_C2f(eps, &mut _C2a, geod.GEODESIC_ORDER);
158            _B21 = geomath::sin_cos_series(true, _ssig1, _csig1, &_C2a);
159        }
160
161        let mut _C3a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER];
162        let mut _A3c = 0.0;
163        let mut _B31 = 0.0;
164        if caps & caps::CAP_C3 != 0 {
165            geod._C3f(eps, &mut _C3a);
166            _A3c = -f * _salp0 * geod._A3f(eps);
167            _B31 = geomath::sin_cos_series(true, _ssig1, _csig1, &_C3a);
168        }
169
170        let mut _C4a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER];
171        let mut _A4 = 0.0;
172        let mut _B41 = 0.0;
173        if caps & caps::CAP_C4 != 0 {
174            geod._C4f(eps, &mut _C4a);
175            _A4 = geomath::sq(_a) * _calp0 * _salp0 * geod._e2;
176            _B41 = geomath::sin_cos_series(false, _ssig1, _csig1, &_C4a);
177        }
178
179        let _s13 = f64::NAN;
180        let _a13 = f64::NAN;
181
182        GeodesicLine {
183            tiny_,
184            _A1m1,
185            _A2m1,
186            _A3c,
187            _A4,
188            _B11,
189            _B21,
190            _B31,
191            _B41,
192            _C1a,
193            _C1pa,
194            _comg1,
195            _C2a,
196            _C3a,
197            _C4a,
198            _b,
199            _c2,
200            _calp0,
201            _csig1,
202            _ctau1,
203            _dn1,
204            _f1,
205            _k2,
206            _salp0,
207            _somg1,
208            _ssig1,
209            _stau1,
210            _a,
211            _a13,
212            azi1,
213            calp1,
214            caps,
215            f,
216            lat1,
217            lon1,
218            _s13,
219            salp1,
220        }
221    }
222
223    /// Place a point at a given distance along this line.
224    ///
225    /// # Arguments
226    ///   - arcmode - Indicates if s12_a12 is an arc length (true) or distance (false)
227    ///   - s12_a12 - Distance to point (meters) or arc length (degrees); can be negative
228    ///   - outmask - bitor'ed combination of [capabilities](caps) defining which values to return
229    ///     - `outmask |= LATITUDE` for the latitude lat2
230    ///     - `outmask |= LONGITUDE` for the longitude lon2
231    ///     - `outmask |= AZIMUTH` for the azimuth azi2
232    ///     - `outmask |= DISTANCE` for the distance s12
233    ///     - `outmask |= REDUCEDLENGTH` for the reduced length m12
234    ///     - `outmask |= GEODESICSCALE` for the geodesic scales M12 and M21
235    ///     - `outmask |= ALL` for all of the above
236    ///     - `outmask |= LONG_UNROLL` to unroll lon2 (instead of reducing it to the range [−180°, 180°])
237    ///
238    /// # Returns
239    /// Values are [f64::NAN] if not supported by the specified capabilities.
240    ///  - a12 arc length between point 1 and point 2 (degrees)
241    ///  - lat2 latitude of point 2 (degrees)
242    ///  - lon2 longitude of point 2 (degrees)
243    ///  - azi2 (forward) azimuth at point 2 (degrees)
244    ///  - s12 distance between point 1 and point 2 (meters)
245    ///  - m12 reduced length of geodesic (meters)
246    ///  - M12 geodesic scale of point 2 relative to point 1 (dimensionless)
247    ///  - M21 geodesic scale of point 1 relative to point 2 (dimensionless)
248    ///  - S12 area under the geodesic (meters<sup>2</sup>)
249    ///
250    /// # Example
251    /// ```rust
252    /// use geographiclib_rs::{Geodesic, GeodesicLine, geodesic_capability as caps};
253    /// let g = Geodesic::wgs84();
254    /// let outmask = caps::LATITUDE | caps::LONGITUDE | caps::AZIMUTH | caps::DISTANCE_IN;
255    /// let line = GeodesicLine::new(&g, -11.95887, -116.94513, 92.712619830452549, Some(outmask), None, None);
256    /// let (_, lat2, lon2, azi2, _, _, _, _, _) = line._gen_position(false, 13834722.5801401374, outmask);
257    ///
258    /// use approx::assert_relative_eq;
259    /// assert_relative_eq!(lat2, 4.57352, epsilon=1e-13);
260    /// assert_relative_eq!(lon2, 7.16501, epsilon=1e-13);
261    /// assert_relative_eq!(azi2, 78.64960934409585, epsilon=1e-13);
262    /// ```
263    pub fn _gen_position(
264        &self,
265        arcmode: bool,
266        s12_a12: f64,
267        outmask: u64,
268    ) -> (f64, f64, f64, f64, f64, f64, f64, f64, f64) {
269        let mut a12 = f64::NAN;
270        let mut lat2 = f64::NAN;
271        let mut lon2 = f64::NAN;
272        let mut azi2 = f64::NAN;
273        let mut s12 = f64::NAN;
274        let mut m12 = f64::NAN;
275        let mut M12 = f64::NAN;
276        let mut M21 = f64::NAN;
277        let mut S12 = f64::NAN;
278        let outmask = outmask & (self.caps & caps::OUT_MASK);
279        if !(arcmode || (self.caps & (caps::OUT_MASK & caps::DISTANCE_IN) != 0)) {
280            return (a12, lat2, lon2, azi2, s12, m12, M12, M21, S12);
281        }
282
283        let mut B12 = 0.0;
284        let mut AB1 = 0.0;
285        let mut sig12: f64;
286        let mut ssig12: f64;
287        let mut csig12: f64;
288        let mut ssig2: f64;
289        let mut csig2: f64;
290        if arcmode {
291            sig12 = s12_a12.to_radians();
292            let res = geomath::sincosd(s12_a12);
293            ssig12 = res.0;
294            csig12 = res.1;
295        } else {
296            // tau12 = s12_a12 / (self._b * (1 + self._A1m1))
297            let tau12 = s12_a12 / (self._b * (1.0 + self._A1m1));
298
299            let s = tau12.sin();
300            let c = tau12.cos();
301
302            B12 = -geomath::sin_cos_series(
303                true,
304                self._stau1 * c + self._ctau1 * s,
305                self._ctau1 * c - self._stau1 * s,
306                &self._C1pa,
307            );
308            sig12 = tau12 - (B12 - self._B11);
309            ssig12 = sig12.sin();
310            csig12 = sig12.cos();
311            if self.f.abs() > 0.01 {
312                ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12;
313                csig2 = self._csig1 * csig12 - self._ssig1 * ssig12;
314                B12 = geomath::sin_cos_series(true, ssig2, csig2, &self._C1a);
315                let serr = (1.0 + self._A1m1) * (sig12 + (B12 - self._B11)) - s12_a12 / self._b;
316                sig12 -= serr / (1.0 + self._k2 * geomath::sq(ssig2)).sqrt();
317                ssig12 = sig12.sin();
318                csig12 = sig12.cos();
319            }
320        };
321        ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12;
322        csig2 = self._csig1 * csig12 - self._ssig1 * ssig12;
323        let dn2 = (1.0 + self._k2 * geomath::sq(ssig2)).sqrt();
324        if outmask & (caps::DISTANCE | caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 {
325            if arcmode || self.f.abs() > 0.01 {
326                B12 = geomath::sin_cos_series(true, ssig2, csig2, &self._C1a);
327            }
328            AB1 = (1.0 + self._A1m1) * (B12 - self._B11);
329        }
330
331        let sbet2 = self._calp0 * ssig2;
332        let mut cbet2 = self._salp0.hypot(self._calp0 * csig2);
333        if cbet2 == 0.0 {
334            cbet2 = self.tiny_;
335            csig2 = self.tiny_;
336        }
337        let salp2 = self._salp0;
338        let calp2 = self._calp0 * csig2;
339        if outmask & caps::DISTANCE != 0 {
340            s12 = if arcmode {
341                self._b * ((1.0 + self._A1m1) * sig12 + AB1)
342            } else {
343                s12_a12
344            }
345        }
346        if outmask & caps::LONGITUDE != 0 {
347            let somg2 = self._salp0 * ssig2;
348            let comg2 = csig2;
349            let E = 1.0_f64.copysign(self._salp0);
350            let omg12 = if outmask & caps::LONG_UNROLL != 0 {
351                E * (sig12 - (ssig2.atan2(csig2) - self._ssig1.atan2(self._csig1))
352                    + ((E * somg2).atan2(comg2) - (E * self._somg1).atan2(self._comg1)))
353            } else {
354                (somg2 * self._comg1 - comg2 * self._somg1)
355                    .atan2(comg2 * self._comg1 + somg2 * self._somg1)
356            };
357            let lam12 = omg12
358                + self._A3c
359                    * (sig12
360                        + (geomath::sin_cos_series(true, ssig2, csig2, &self._C3a) - self._B31));
361            let lon12 = lam12.to_degrees();
362            lon2 = if outmask & caps::LONG_UNROLL != 0 {
363                self.lon1 + lon12
364            } else {
365                geomath::ang_normalize(
366                    geomath::ang_normalize(self.lon1) + geomath::ang_normalize(lon12),
367                )
368            };
369        };
370
371        if outmask & caps::LATITUDE != 0 {
372            lat2 = geomath::atan2d(sbet2, self._f1 * cbet2);
373        }
374        if outmask & caps::AZIMUTH != 0 {
375            azi2 = geomath::atan2d(salp2, calp2);
376        }
377        if outmask & (caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 {
378            let B22 = geomath::sin_cos_series(true, ssig2, csig2, &self._C2a);
379            let AB2 = (1.0 + self._A2m1) * (B22 - self._B21);
380            let J12 = (self._A1m1 - self._A2m1) * sig12 + (AB1 - AB2);
381            if outmask & caps::REDUCEDLENGTH != 0 {
382                m12 = self._b
383                    * ((dn2 * (self._csig1 * ssig2) - self._dn1 * (self._ssig1 * csig2))
384                        - self._csig1 * csig2 * J12);
385            }
386            if outmask & caps::GEODESICSCALE != 0 {
387                let t =
388                    self._k2 * (ssig2 - self._ssig1) * (ssig2 + self._ssig1) / (self._dn1 + dn2);
389                M12 = csig12 + (t * ssig2 - csig2 * J12) * self._ssig1 / self._dn1;
390                M21 = csig12 - (t * self._ssig1 - self._csig1 * J12) * ssig2 / dn2;
391            }
392        }
393        if outmask & caps::AREA != 0 {
394            let B42 = geomath::sin_cos_series(false, ssig2, csig2, &self._C4a);
395            let salp12: f64;
396            let calp12: f64;
397            if self._calp0 == 0.0 || self._salp0 == 0.0 {
398                salp12 = salp2 * self.calp1 - calp2 * self.salp1;
399                calp12 = calp2 * self.calp1 + salp2 * self.salp1;
400            } else {
401                salp12 = self._calp0
402                    * self._salp0
403                    * (if csig12 <= 0.0 {
404                        self._csig1 * (1.0 - csig12) + ssig12 * self._ssig1
405                    } else {
406                        ssig12 * (self._csig1 * ssig12 / (1.0 + csig12) + self._ssig1)
407                    });
408                calp12 = geomath::sq(self._salp0) + geomath::sq(self._calp0) * self._csig1 * csig2;
409            }
410            S12 = self._c2 * salp12.atan2(calp12) + self._A4 * (B42 - self._B41);
411        }
412        a12 = if arcmode { s12_a12 } else { sig12.to_degrees() };
413        (a12, lat2, lon2, azi2, s12, m12, M12, M21, S12)
414    }
415
416    // not currently used, but maybe some day
417    #[allow(dead_code)]
418    pub fn Position(&self, s12: f64, outmask: Option<u64>) -> HashMap<String, f64> {
419        let outmask = match outmask {
420            Some(outmask) => outmask,
421            None => caps::STANDARD,
422        };
423        let mut result: HashMap<String, f64> = HashMap::new();
424        result.insert("lat1".to_string(), self.lat1);
425        result.insert("azi1".to_string(), self.azi1);
426        result.insert("s12".to_string(), s12);
427        let lon1 = if outmask & caps::LONG_UNROLL != 0 {
428            self.lon1
429        } else {
430            geomath::ang_normalize(self.lon1)
431        };
432        result.insert("lon1".to_string(), lon1);
433
434        let (a12, lat2, lon2, azi2, _s12, m12, M12, M21, S12) =
435            self._gen_position(false, s12, outmask);
436        let outmask = outmask & caps::OUT_MASK;
437        result.insert("a12".to_string(), a12);
438        if outmask & caps::LATITUDE != 0 {
439            result.insert("lat2".to_string(), lat2);
440        }
441        if outmask & caps::LONGITUDE != 0 {
442            result.insert("lon2".to_string(), lon2);
443        }
444        if outmask & caps::AZIMUTH != 0 {
445            result.insert("azi2".to_string(), azi2);
446        }
447        if outmask & caps::REDUCEDLENGTH != 0 {
448            result.insert("m12".to_string(), m12);
449        }
450        if outmask & caps::GEODESICSCALE != 0 {
451            result.insert("M12".to_string(), M12);
452            result.insert("M21".to_string(), M21);
453        }
454        if outmask & caps::AREA != 0 {
455            result.insert("S12".to_string(), S12);
456        }
457        result
458    }
459}
460
461#[cfg(test)]
462mod tests {
463    use super::*;
464    use geodesic::Geodesic;
465
466    #[test]
467    fn test_gen_position() {
468        let geod = Geodesic::wgs84();
469        let gl = GeodesicLine::new(&geod, 0.0, 0.0, 10.0, None, None, None);
470        let res = gl._gen_position(false, 150.0, 3979);
471        assert_eq!(res.0, 0.0013520059461334633);
472        assert_eq!(res.1, 0.0013359451088740494);
473        assert_eq!(res.2, 0.00023398621812867812);
474        assert_eq!(res.3, 10.000000002727887);
475        assert_eq!(res.4, 150.0);
476        assert!(res.5.is_nan());
477        assert!(res.6.is_nan());
478        assert!(res.7.is_nan());
479        assert!(res.8.is_nan());
480    }
481
482    #[test]
483    fn test_init() {
484        let geod = Geodesic::wgs84();
485        let gl = GeodesicLine::new(&geod, 0.0, 0.0, 0.0, None, None, None);
486        assert_eq!(gl._a, 6378137.0);
487        assert_eq!(gl.f, 0.0033528106647474805);
488        assert_eq!(gl._b, 6356752.314245179);
489        assert_eq!(gl._c2, 40589732499314.76);
490        assert_eq!(gl._f1, 0.9966471893352525);
491        assert_eq!(gl.caps, 36747);
492        assert_eq!(gl.lat1, 0.0);
493        assert_eq!(gl.lon1, 0.0);
494        assert_eq!(gl.azi1, 0.0);
495        assert_eq!(gl.salp1, 0.0);
496        assert_eq!(gl.calp1, 1.0);
497        assert_eq!(gl._dn1, 1.0);
498        assert_eq!(gl._salp0, 0.0);
499        assert_eq!(gl._calp0, 1.0);
500        assert_eq!(gl._ssig1, 0.0);
501        assert_eq!(gl._somg1, 0.0);
502        assert_eq!(gl._csig1, 1.0);
503        assert_eq!(gl._comg1, 1.0);
504        assert_eq!(gl._k2, geod._ep2);
505        assert!(gl._s13.is_nan());
506        assert!(gl._a13.is_nan());
507    }
508}
geographiclib_rs/geodesic_line.rs

geographiclib_rs/
geodesic_line.rs
