geographiclib_rs/
geodesic_line.rs

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

geographiclib_rs/
geodesic_line.rs
