/rust/registry/src/index.crates.io-1949cf8c6b5b557f/geographiclib-rs-0.2.7/src/geodesic.rs
Line | Count | Source |
1 | | #![allow(non_snake_case)] |
2 | | #![allow(clippy::excessive_precision)] |
3 | | |
4 | | use crate::geodesic_capability as caps; |
5 | | use crate::geodesic_line; |
6 | | use crate::geomath; |
7 | | use std::sync; |
8 | | |
9 | | use std::f64::consts::{FRAC_1_SQRT_2, PI}; |
10 | | |
11 | | pub const WGS84_A: f64 = 6378137.0; |
12 | | // Evaluating this as 1000000000.0 / (298257223563f64) reduces the |
13 | | // round-off error by about 10%. However, expressing the flattening as |
14 | | // 1/298.257223563 is well ingrained. |
15 | | pub const WGS84_F: f64 = 1.0 / ((298257223563f64) / 1000000000.0); |
16 | | |
17 | | #[derive(Copy, Clone, PartialEq, PartialOrd, Debug)] |
18 | | pub struct Geodesic { |
19 | | pub a: f64, |
20 | | pub f: f64, |
21 | | pub _f1: f64, |
22 | | pub _e2: f64, |
23 | | pub _ep2: f64, |
24 | | _n: f64, |
25 | | pub _b: f64, |
26 | | pub _c2: f64, |
27 | | _etol2: f64, |
28 | | _A3x: [f64; GEODESIC_ORDER], |
29 | | _C3x: [f64; _nC3x_], |
30 | | _C4x: [f64; _nC4x_], |
31 | | |
32 | | _nC3x_: usize, |
33 | | _nC4x_: usize, |
34 | | maxit1_: u64, |
35 | | maxit2_: u64, |
36 | | |
37 | | pub tiny_: f64, |
38 | | tol0_: f64, |
39 | | tol1_: f64, |
40 | | _tol2_: f64, |
41 | | tolb_: f64, |
42 | | xthresh_: f64, |
43 | | } |
44 | | |
45 | | static WGS84_GEOD: sync::OnceLock<Geodesic> = sync::OnceLock::new(); |
46 | | |
47 | | impl Geodesic { |
48 | 0 | pub fn wgs84() -> Self { |
49 | 0 | *WGS84_GEOD.get_or_init(|| Geodesic::new(WGS84_A, WGS84_F)) |
50 | 0 | } |
51 | | |
52 | 0 | pub fn equatorial_radius(&self) -> f64 { |
53 | 0 | self.a |
54 | 0 | } |
55 | | |
56 | 0 | pub fn flattening(&self) -> f64 { |
57 | 0 | self.f |
58 | 0 | } |
59 | | } |
60 | | |
61 | | const COEFF_A3: [f64; 18] = [ |
62 | | -3.0, 128.0, -2.0, -3.0, 64.0, -1.0, -3.0, -1.0, 16.0, 3.0, -1.0, -2.0, 8.0, 1.0, -1.0, 2.0, |
63 | | 1.0, 1.0, |
64 | | ]; |
65 | | |
66 | | const COEFF_C3: [f64; 45] = [ |
67 | | 3.0, 128.0, 2.0, 5.0, 128.0, -1.0, 3.0, 3.0, 64.0, -1.0, 0.0, 1.0, 8.0, -1.0, 1.0, 4.0, 5.0, |
68 | | 256.0, 1.0, 3.0, 128.0, -3.0, -2.0, 3.0, 64.0, 1.0, -3.0, 2.0, 32.0, 7.0, 512.0, -10.0, 9.0, |
69 | | 384.0, 5.0, -9.0, 5.0, 192.0, 7.0, 512.0, -14.0, 7.0, 512.0, 21.0, 2560.0, |
70 | | ]; |
71 | | |
72 | | const COEFF_C4: [f64; 77] = [ |
73 | | 97.0, 15015.0, 1088.0, 156.0, 45045.0, -224.0, -4784.0, 1573.0, 45045.0, -10656.0, 14144.0, |
74 | | -4576.0, -858.0, 45045.0, 64.0, 624.0, -4576.0, 6864.0, -3003.0, 15015.0, 100.0, 208.0, 572.0, |
75 | | 3432.0, -12012.0, 30030.0, 45045.0, 1.0, 9009.0, -2944.0, 468.0, 135135.0, 5792.0, 1040.0, |
76 | | -1287.0, 135135.0, 5952.0, -11648.0, 9152.0, -2574.0, 135135.0, -64.0, -624.0, 4576.0, -6864.0, |
77 | | 3003.0, 135135.0, 8.0, 10725.0, 1856.0, -936.0, 225225.0, -8448.0, 4992.0, -1144.0, 225225.0, |
78 | | -1440.0, 4160.0, -4576.0, 1716.0, 225225.0, -136.0, 63063.0, 1024.0, -208.0, 105105.0, 3584.0, |
79 | | -3328.0, 1144.0, 315315.0, -128.0, 135135.0, -2560.0, 832.0, 405405.0, 128.0, 99099.0, |
80 | | ]; |
81 | | |
82 | | pub const GEODESIC_ORDER: usize = 6; |
83 | | pub(crate) const CARR_SIZE: usize = GEODESIC_ORDER + 1; |
84 | | |
85 | | #[allow(non_upper_case_globals)] |
86 | | const _nC3x_: usize = 15; |
87 | | #[allow(non_upper_case_globals)] |
88 | | const _nC4x_: usize = 21; |
89 | | |
90 | | impl Geodesic { |
91 | 0 | pub fn new(a: f64, f: f64) -> Self { |
92 | 0 | let maxit1_ = 20; |
93 | 0 | let maxit2_ = maxit1_ + geomath::DIGITS + 10; |
94 | 0 | let tiny_ = f64::MIN_POSITIVE.sqrt(); |
95 | 0 | let tol0_ = f64::EPSILON; |
96 | 0 | let tol1_ = 200.0 * tol0_; |
97 | 0 | let _tol2_ = tol0_.sqrt(); |
98 | 0 | let tolb_ = tol0_ * _tol2_; |
99 | 0 | let xthresh_ = 1000.0 * _tol2_; |
100 | | |
101 | 0 | let _f1 = 1.0 - f; |
102 | 0 | let _e2 = f * (2.0 - f); |
103 | 0 | let _ep2 = _e2 / geomath::sq(_f1); |
104 | 0 | let _n = f / (2.0 - f); |
105 | 0 | let _b = a * _f1; |
106 | 0 | let _c2 = (geomath::sq(a) |
107 | 0 | + geomath::sq(_b) |
108 | 0 | * (if _e2 == 0.0 { |
109 | 0 | 1.0 |
110 | | } else { |
111 | 0 | geomath::eatanhe(1.0, (if f < 0.0 { -1.0 } else { 1.0 }) * _e2.abs().sqrt()) |
112 | 0 | / _e2 |
113 | | })) |
114 | | / 2.0; |
115 | 0 | let _etol2 = 0.1 * _tol2_ / (f.abs().max(0.001) * (1.0 - f / 2.0).min(1.0) / 2.0).sqrt(); |
116 | | |
117 | 0 | let mut _A3x: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER]; |
118 | 0 | let mut _C3x: [f64; _nC3x_] = [0.0; _nC3x_]; |
119 | 0 | let mut _C4x: [f64; _nC4x_] = [0.0; _nC4x_]; |
120 | | |
121 | | // Call a3coeff |
122 | 0 | let mut o: usize = 0; |
123 | 0 | for (k, j) in (0..GEODESIC_ORDER).rev().enumerate() { |
124 | 0 | let m = j.min(GEODESIC_ORDER - j - 1); |
125 | 0 | _A3x[k] = geomath::polyval(m, &COEFF_A3[o..], _n) / COEFF_A3[o + m + 1]; |
126 | 0 | o += m + 2; |
127 | 0 | } |
128 | | |
129 | | // c3coeff |
130 | 0 | let mut o = 0; |
131 | 0 | let mut k = 0; |
132 | 0 | for l in 1..GEODESIC_ORDER { |
133 | 0 | for j in (l..GEODESIC_ORDER).rev() { |
134 | 0 | let m = j.min(GEODESIC_ORDER - j - 1); |
135 | 0 | _C3x[k] = geomath::polyval(m, &COEFF_C3[o..], _n) / COEFF_C3[o + m + 1]; |
136 | 0 | k += 1; |
137 | 0 | o += m + 2; |
138 | 0 | } |
139 | | } |
140 | | |
141 | | // c4coeff |
142 | 0 | let mut o = 0; |
143 | 0 | let mut k = 0; |
144 | 0 | for l in 0..GEODESIC_ORDER { |
145 | 0 | for j in (l..GEODESIC_ORDER).rev() { |
146 | 0 | let m = GEODESIC_ORDER - j - 1; |
147 | 0 | _C4x[k] = geomath::polyval(m, &COEFF_C4[o..], _n) / COEFF_C4[o + m + 1]; |
148 | 0 | k += 1; |
149 | 0 | o += m + 2; |
150 | 0 | } |
151 | | } |
152 | | |
153 | 0 | Geodesic { |
154 | 0 | a, |
155 | 0 | f, |
156 | 0 | _f1, |
157 | 0 | _e2, |
158 | 0 | _ep2, |
159 | 0 | _n, |
160 | 0 | _b, |
161 | 0 | _c2, |
162 | 0 | _etol2, |
163 | 0 | _A3x, |
164 | 0 | _C3x, |
165 | 0 | _C4x, |
166 | 0 |
|
167 | 0 | _nC3x_, |
168 | 0 | _nC4x_, |
169 | 0 | maxit1_, |
170 | 0 | maxit2_, |
171 | 0 |
|
172 | 0 | tiny_, |
173 | 0 | tol0_, |
174 | 0 | tol1_, |
175 | 0 | _tol2_, |
176 | 0 | tolb_, |
177 | 0 | xthresh_, |
178 | 0 | } |
179 | 0 | } |
180 | | |
181 | 0 | pub fn _A3f(&self, eps: f64) -> f64 { |
182 | 0 | geomath::polyval(GEODESIC_ORDER - 1, &self._A3x, eps) |
183 | 0 | } |
184 | | |
185 | 0 | pub fn _C3f(&self, eps: f64, c: &mut [f64; GEODESIC_ORDER]) { |
186 | 0 | let mut mult = 1.0; |
187 | 0 | let mut o = 0; |
188 | | // Clippy wants us to turn this into `c.iter_mut().enumerate().take(geodesic_order + 1).skip(1)` |
189 | | // but benching (rust-1.75) shows that it would be slower. |
190 | | #[allow(clippy::needless_range_loop)] |
191 | 0 | for l in 1..GEODESIC_ORDER { |
192 | 0 | let m = GEODESIC_ORDER - l - 1; |
193 | 0 | mult *= eps; |
194 | 0 | c[l] = mult * geomath::polyval(m, &self._C3x[o..], eps); |
195 | 0 | o += m + 1; |
196 | 0 | } |
197 | 0 | } |
198 | | |
199 | 0 | pub fn _C4f(&self, eps: f64, c: &mut [f64; GEODESIC_ORDER]) { |
200 | 0 | let mut mult = 1.0; |
201 | 0 | let mut o = 0; |
202 | | // Clippy wants us to turn this into `c.iter_mut().enumerate().take(geodesic_order + 1).skip(1)` |
203 | | // but benching (rust-1.75) shows that it would be slower. |
204 | | #[allow(clippy::needless_range_loop)] |
205 | 0 | for l in 0..GEODESIC_ORDER { |
206 | 0 | let m = GEODESIC_ORDER - l - 1; |
207 | 0 | c[l] = mult * geomath::polyval(m, &self._C4x[o..], eps); |
208 | 0 | o += m + 1; |
209 | 0 | mult *= eps; |
210 | 0 | } |
211 | 0 | } |
212 | | |
213 | | #[allow(clippy::too_many_arguments)] |
214 | 0 | pub fn _Lengths( |
215 | 0 | &self, |
216 | 0 | eps: f64, |
217 | 0 | sig12: f64, |
218 | 0 | ssig1: f64, |
219 | 0 | csig1: f64, |
220 | 0 | dn1: f64, |
221 | 0 | ssig2: f64, |
222 | 0 | csig2: f64, |
223 | 0 | dn2: f64, |
224 | 0 | cbet1: f64, |
225 | 0 | cbet2: f64, |
226 | 0 | outmask: u64, |
227 | 0 | C1a: &mut [f64; CARR_SIZE], |
228 | 0 | C2a: &mut [f64; CARR_SIZE], |
229 | 0 | ) -> (f64, f64, f64, f64, f64) { |
230 | 0 | let outmask = outmask & caps::OUT_MASK; |
231 | 0 | let mut s12b = f64::NAN; |
232 | 0 | let mut m12b = f64::NAN; |
233 | 0 | let mut m0 = f64::NAN; |
234 | 0 | let mut M12 = f64::NAN; |
235 | 0 | let mut M21 = f64::NAN; |
236 | | |
237 | 0 | let mut A1 = 0.0; |
238 | 0 | let mut A2 = 0.0; |
239 | 0 | let mut m0x = 0.0; |
240 | 0 | let mut J12 = 0.0; |
241 | | |
242 | 0 | if outmask & (caps::DISTANCE | caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 { |
243 | 0 | A1 = geomath::_A1m1f(eps); |
244 | 0 | geomath::_C1f(eps, C1a); |
245 | 0 | if outmask & (caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 { |
246 | 0 | A2 = geomath::_A2m1f(eps); |
247 | 0 | geomath::_C2f(eps, C2a); |
248 | 0 | m0x = A1 - A2; |
249 | 0 | A2 += 1.0; |
250 | 0 | } |
251 | 0 | A1 += 1.0; |
252 | 0 | } |
253 | 0 | if outmask & caps::DISTANCE != 0 { |
254 | 0 | let B1 = geomath::sin_cos_series(true, ssig2, csig2, C1a) |
255 | 0 | - geomath::sin_cos_series(true, ssig1, csig1, C1a); |
256 | 0 | s12b = A1 * (sig12 + B1); |
257 | 0 | if outmask & (caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 { |
258 | 0 | let B2 = geomath::sin_cos_series(true, ssig2, csig2, C2a) |
259 | 0 | - geomath::sin_cos_series(true, ssig1, csig1, C2a); |
260 | 0 | J12 = m0x * sig12 + (A1 * B1 - A2 * B2); |
261 | 0 | } |
262 | 0 | } else if outmask & (caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 { |
263 | 0 | for l in 1..=GEODESIC_ORDER { |
264 | 0 | C2a[l] = A1 * C1a[l] - A2 * C2a[l]; |
265 | 0 | } |
266 | 0 | J12 = m0x * sig12 |
267 | 0 | + (geomath::sin_cos_series(true, ssig2, csig2, C2a) |
268 | 0 | - geomath::sin_cos_series(true, ssig1, csig1, C2a)); |
269 | 0 | } |
270 | 0 | if outmask & caps::REDUCEDLENGTH != 0 { |
271 | 0 | m0 = m0x; |
272 | 0 | // J12 is wrong |
273 | 0 | m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12; |
274 | 0 | } |
275 | 0 | if outmask & caps::GEODESICSCALE != 0 { |
276 | 0 | let csig12 = csig1 * csig2 + ssig1 * ssig2; |
277 | 0 | let t = self._ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2); |
278 | 0 | M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1; |
279 | 0 | M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2; |
280 | 0 | } |
281 | 0 | (s12b, m12b, m0, M12, M21) |
282 | 0 | } |
283 | | |
284 | | #[allow(clippy::too_many_arguments)] |
285 | 0 | pub fn _InverseStart( |
286 | 0 | &self, |
287 | 0 | sbet1: f64, |
288 | 0 | cbet1: f64, |
289 | 0 | dn1: f64, |
290 | 0 | sbet2: f64, |
291 | 0 | cbet2: f64, |
292 | 0 | dn2: f64, |
293 | 0 | lam12: f64, |
294 | 0 | slam12: f64, |
295 | 0 | clam12: f64, |
296 | 0 | C1a: &mut [f64; CARR_SIZE], |
297 | 0 | C2a: &mut [f64; CARR_SIZE], |
298 | 0 | ) -> (f64, f64, f64, f64, f64, f64) { |
299 | 0 | let mut sig12 = -1.0; |
300 | 0 | let mut salp2 = f64::NAN; |
301 | 0 | let mut calp2 = f64::NAN; |
302 | 0 | let mut dnm = f64::NAN; |
303 | | |
304 | | let mut somg12: f64; |
305 | | let mut comg12: f64; |
306 | | |
307 | 0 | let sbet12 = sbet2 * cbet1 - cbet2 * sbet1; |
308 | 0 | let cbet12 = cbet2 * cbet1 + sbet2 * sbet1; |
309 | | |
310 | 0 | let mut sbet12a = sbet2 * cbet1; |
311 | 0 | sbet12a += cbet2 * sbet1; |
312 | | |
313 | 0 | let shortline = cbet12 >= 0.0 && sbet12 < 0.5 && cbet2 * lam12 < 0.5; |
314 | 0 | if shortline { |
315 | 0 | let mut sbetm2 = geomath::sq(sbet1 + sbet2); |
316 | 0 | sbetm2 /= sbetm2 + geomath::sq(cbet1 + cbet2); |
317 | 0 | dnm = (1.0 + self._ep2 * sbetm2).sqrt(); |
318 | 0 | let omg12 = lam12 / (self._f1 * dnm); |
319 | 0 | somg12 = omg12.sin(); |
320 | 0 | comg12 = omg12.cos(); |
321 | 0 | } else { |
322 | 0 | somg12 = slam12; |
323 | 0 | comg12 = clam12; |
324 | 0 | } |
325 | | |
326 | 0 | let mut salp1 = cbet2 * somg12; |
327 | | |
328 | 0 | let mut calp1 = if comg12 >= 0.0 { |
329 | 0 | sbet12 + cbet2 * sbet1 * geomath::sq(somg12) / (1.0 + comg12) |
330 | | } else { |
331 | 0 | sbet12a - cbet2 * sbet1 * geomath::sq(somg12) / (1.0 - comg12) |
332 | | }; |
333 | | |
334 | 0 | let ssig12 = salp1.hypot(calp1); |
335 | 0 | let csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12; |
336 | | |
337 | 0 | if shortline && ssig12 < self._etol2 { |
338 | 0 | salp2 = cbet1 * somg12; |
339 | 0 | calp2 = sbet12 |
340 | 0 | - cbet1 |
341 | 0 | * sbet2 |
342 | 0 | * (if comg12 >= 0.0 { |
343 | 0 | geomath::sq(somg12) / (1.0 + comg12) |
344 | | } else { |
345 | 0 | 1.0 - comg12 |
346 | | }); |
347 | 0 | geomath::norm(&mut salp2, &mut calp2); |
348 | 0 | sig12 = ssig12.atan2(csig12); |
349 | 0 | } else if self._n.abs() > 0.1 |
350 | 0 | || csig12 >= 0.0 |
351 | 0 | || ssig12 >= 6.0 * self._n.abs() * PI * geomath::sq(cbet1) |
352 | 0 | { |
353 | 0 | } else { |
354 | | let x: f64; |
355 | | let y: f64; |
356 | | let betscale: f64; |
357 | | let lamscale: f64; |
358 | 0 | let lam12x = (-slam12).atan2(-clam12); |
359 | 0 | if self.f >= 0.0 { |
360 | 0 | let k2 = geomath::sq(sbet1) * self._ep2; |
361 | 0 | let eps = k2 / (2.0 * (1.0 + (1.0 + k2).sqrt()) + k2); |
362 | 0 | lamscale = self.f * cbet1 * self._A3f(eps) * PI; |
363 | 0 | betscale = lamscale * cbet1; |
364 | 0 | x = lam12x / lamscale; |
365 | 0 | y = sbet12a / betscale; |
366 | 0 | } else { |
367 | 0 | let cbet12a = cbet2 * cbet1 - sbet2 * sbet1; |
368 | 0 | let bet12a = sbet12a.atan2(cbet12a); |
369 | 0 | let (_, m12b, m0, _, _) = self._Lengths( |
370 | 0 | self._n, |
371 | 0 | PI + bet12a, |
372 | 0 | sbet1, |
373 | 0 | -cbet1, |
374 | 0 | dn1, |
375 | 0 | sbet2, |
376 | 0 | cbet2, |
377 | 0 | dn2, |
378 | 0 | cbet1, |
379 | 0 | cbet2, |
380 | 0 | caps::REDUCEDLENGTH, |
381 | 0 | C1a, |
382 | 0 | C2a, |
383 | 0 | ); |
384 | 0 | x = -1.0 + m12b / (cbet1 * cbet2 * m0 * PI); |
385 | 0 | betscale = if x < -0.01 { |
386 | 0 | sbet12a / x |
387 | | } else { |
388 | 0 | -self.f * geomath::sq(cbet1) * PI |
389 | | }; |
390 | 0 | lamscale = betscale / cbet1; |
391 | 0 | y = lam12x / lamscale; |
392 | | } |
393 | 0 | if y > -self.tol1_ && x > -1.0 - self.xthresh_ { |
394 | 0 | if self.f >= 0.0 { |
395 | 0 | salp1 = (-x).min(1.0); |
396 | 0 | calp1 = -(1.0 - geomath::sq(salp1)).sqrt() |
397 | | } else { |
398 | 0 | calp1 = x.max(if x > -self.tol1_ { 0.0 } else { -1.0 }); |
399 | 0 | salp1 = (1.0 - geomath::sq(calp1)).sqrt(); |
400 | | } |
401 | | } else { |
402 | 0 | let k = geomath::astroid(x, y); |
403 | 0 | let omg12a = lamscale |
404 | 0 | * if self.f >= 0.0 { |
405 | 0 | -x * k / (1.0 + k) |
406 | | } else { |
407 | 0 | -y * (1.0 + k) / k |
408 | | }; |
409 | 0 | somg12 = omg12a.sin(); |
410 | 0 | comg12 = -(omg12a.cos()); |
411 | 0 | salp1 = cbet2 * somg12; |
412 | 0 | calp1 = sbet12a - cbet2 * sbet1 * geomath::sq(somg12) / (1.0 - comg12); |
413 | | } |
414 | | } |
415 | | |
416 | 0 | if salp1 > 0.0 || salp1.is_nan() { |
417 | 0 | geomath::norm(&mut salp1, &mut calp1); |
418 | 0 | } else { |
419 | 0 | salp1 = 1.0; |
420 | 0 | calp1 = 0.0; |
421 | 0 | }; |
422 | 0 | (sig12, salp1, calp1, salp2, calp2, dnm) |
423 | 0 | } |
424 | | |
425 | | #[allow(clippy::too_many_arguments)] |
426 | 0 | pub fn _Lambda12( |
427 | 0 | &self, |
428 | 0 | sbet1: f64, |
429 | 0 | cbet1: f64, |
430 | 0 | dn1: f64, |
431 | 0 | sbet2: f64, |
432 | 0 | cbet2: f64, |
433 | 0 | dn2: f64, |
434 | 0 | salp1: f64, |
435 | 0 | mut calp1: f64, |
436 | 0 | slam120: f64, |
437 | 0 | clam120: f64, |
438 | 0 | diffp: bool, |
439 | 0 | C1a: &mut [f64; CARR_SIZE], |
440 | 0 | C2a: &mut [f64; CARR_SIZE], |
441 | 0 | C3a: &mut [f64; GEODESIC_ORDER], |
442 | 0 | ) -> (f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64) { |
443 | 0 | if sbet1 == 0.0 && calp1 == 0.0 { |
444 | 0 | calp1 = -self.tiny_; |
445 | 0 | } |
446 | 0 | let salp0 = salp1 * cbet1; |
447 | 0 | let calp0 = calp1.hypot(salp1 * sbet1); |
448 | | |
449 | 0 | let mut ssig1 = sbet1; |
450 | 0 | let somg1 = salp0 * sbet1; |
451 | 0 | let mut csig1 = calp1 * cbet1; |
452 | 0 | let comg1 = calp1 * cbet1; |
453 | 0 | geomath::norm(&mut ssig1, &mut csig1); |
454 | | |
455 | 0 | let salp2 = if cbet2 != cbet1 { salp0 / cbet2 } else { salp1 }; |
456 | 0 | let calp2 = if cbet2 != cbet1 || sbet2.abs() != -sbet1 { |
457 | 0 | (geomath::sq(calp1 * cbet1) |
458 | 0 | + if cbet1 < -sbet1 { |
459 | 0 | (cbet2 - cbet1) * (cbet1 + cbet2) |
460 | | } else { |
461 | 0 | (sbet1 - sbet2) * (sbet1 + sbet2) |
462 | | }) |
463 | 0 | .sqrt() |
464 | 0 | / cbet2 |
465 | | } else { |
466 | 0 | calp1.abs() |
467 | | }; |
468 | 0 | let mut ssig2 = sbet2; |
469 | 0 | let somg2 = salp0 * sbet2; |
470 | 0 | let mut csig2 = calp2 * cbet2; |
471 | 0 | let comg2 = calp2 * cbet2; |
472 | 0 | geomath::norm(&mut ssig2, &mut csig2); |
473 | | |
474 | 0 | let sig12 = ((csig1 * ssig2 - ssig1 * csig2).max(0.0)).atan2(csig1 * csig2 + ssig1 * ssig2); |
475 | 0 | let somg12 = (comg1 * somg2 - somg1 * comg2).max(0.0); |
476 | 0 | let comg12 = comg1 * comg2 + somg1 * somg2; |
477 | 0 | let eta = (somg12 * clam120 - comg12 * slam120).atan2(comg12 * clam120 + somg12 * slam120); |
478 | | |
479 | 0 | let k2 = geomath::sq(calp0) * self._ep2; |
480 | 0 | let eps = k2 / (2.0 * (1.0 + (1.0 + k2).sqrt()) + k2); |
481 | 0 | self._C3f(eps, C3a); |
482 | 0 | let B312 = geomath::sin_cos_series(true, ssig2, csig2, C3a) |
483 | 0 | - geomath::sin_cos_series(true, ssig1, csig1, C3a); |
484 | 0 | let domg12 = -self.f * self._A3f(eps) * salp0 * (sig12 + B312); |
485 | 0 | let lam12 = eta + domg12; |
486 | | |
487 | | let mut dlam12: f64; |
488 | 0 | if diffp { |
489 | 0 | if calp2 == 0.0 { |
490 | 0 | dlam12 = -2.0 * self._f1 * dn1 / sbet1; |
491 | 0 | } else { |
492 | 0 | let res = self._Lengths( |
493 | 0 | eps, |
494 | 0 | sig12, |
495 | 0 | ssig1, |
496 | 0 | csig1, |
497 | 0 | dn1, |
498 | 0 | ssig2, |
499 | 0 | csig2, |
500 | 0 | dn2, |
501 | 0 | cbet1, |
502 | 0 | cbet2, |
503 | 0 | caps::REDUCEDLENGTH, |
504 | 0 | C1a, |
505 | 0 | C2a, |
506 | 0 | ); |
507 | 0 | dlam12 = res.1; |
508 | 0 | dlam12 *= self._f1 / (calp2 * cbet2); |
509 | 0 | } |
510 | 0 | } else { |
511 | 0 | dlam12 = f64::NAN; |
512 | 0 | } |
513 | 0 | ( |
514 | 0 | lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps, domg12, dlam12, |
515 | 0 | ) |
516 | 0 | } |
517 | | |
518 | | // returns (a12, s12, azi1, azi2, m12, M12, M21, S12) |
519 | 0 | pub fn _gen_inverse_azi( |
520 | 0 | &self, |
521 | 0 | lat1: f64, |
522 | 0 | lon1: f64, |
523 | 0 | lat2: f64, |
524 | 0 | lon2: f64, |
525 | 0 | outmask: u64, |
526 | 0 | ) -> (f64, f64, f64, f64, f64, f64, f64, f64) { |
527 | 0 | let mut azi1 = f64::NAN; |
528 | 0 | let mut azi2 = f64::NAN; |
529 | 0 | let outmask = outmask & caps::OUT_MASK; |
530 | | |
531 | 0 | let (a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12) = |
532 | 0 | self._gen_inverse(lat1, lon1, lat2, lon2, outmask); |
533 | 0 | if outmask & caps::AZIMUTH != 0 { |
534 | 0 | azi1 = geomath::atan2d(salp1, calp1); |
535 | 0 | azi2 = geomath::atan2d(salp2, calp2); |
536 | 0 | } |
537 | 0 | (a12, s12, azi1, azi2, m12, M12, M21, S12) |
538 | 0 | } |
539 | | |
540 | | // returns (a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12) |
541 | 0 | pub fn _gen_inverse( |
542 | 0 | &self, |
543 | 0 | lat1: f64, |
544 | 0 | lon1: f64, |
545 | 0 | lat2: f64, |
546 | 0 | lon2: f64, |
547 | 0 | outmask: u64, |
548 | 0 | ) -> (f64, f64, f64, f64, f64, f64, f64, f64, f64, f64) { |
549 | 0 | let mut lat1 = lat1; |
550 | 0 | let mut lat2 = lat2; |
551 | 0 | let mut a12 = f64::NAN; |
552 | 0 | let mut s12 = f64::NAN; |
553 | 0 | let mut m12 = f64::NAN; |
554 | 0 | let mut M12 = f64::NAN; |
555 | 0 | let mut M21 = f64::NAN; |
556 | 0 | let mut S12 = f64::NAN; |
557 | 0 | let outmask = outmask & caps::OUT_MASK; |
558 | | |
559 | 0 | let (mut lon12, mut lon12s) = geomath::ang_diff(lon1, lon2); |
560 | 0 | let mut lonsign = if lon12 >= 0.0 { 1.0 } else { -1.0 }; |
561 | | |
562 | 0 | lon12 = lonsign * geomath::ang_round(lon12); |
563 | 0 | lon12s = geomath::ang_round((180.0 - lon12) - lonsign * lon12s); |
564 | 0 | let lam12 = lon12.to_radians(); |
565 | | let slam12: f64; |
566 | | let mut clam12: f64; |
567 | 0 | if lon12 > 90.0 { |
568 | 0 | let res = geomath::sincosd(lon12s); |
569 | 0 | slam12 = res.0; |
570 | 0 | clam12 = res.1; |
571 | 0 | clam12 = -clam12; |
572 | 0 | } else { |
573 | 0 | let res = geomath::sincosd(lon12); |
574 | 0 | slam12 = res.0; |
575 | 0 | clam12 = res.1; |
576 | 0 | }; |
577 | 0 | lat1 = geomath::ang_round(geomath::lat_fix(lat1)); |
578 | 0 | lat2 = geomath::ang_round(geomath::lat_fix(lat2)); |
579 | | |
580 | 0 | let swapp = if lat1.abs() < lat2.abs() { -1.0 } else { 1.0 }; |
581 | 0 | if swapp < 0.0 { |
582 | 0 | lonsign *= -1.0; |
583 | 0 | std::mem::swap(&mut lat2, &mut lat1); |
584 | 0 | } |
585 | 0 | let latsign = if lat1 < 0.0 { 1.0 } else { -1.0 }; |
586 | 0 | lat1 *= latsign; |
587 | 0 | lat2 *= latsign; |
588 | | |
589 | 0 | let (mut sbet1, mut cbet1) = geomath::sincosd(lat1); |
590 | 0 | sbet1 *= self._f1; |
591 | | |
592 | 0 | geomath::norm(&mut sbet1, &mut cbet1); |
593 | 0 | cbet1 = cbet1.max(self.tiny_); |
594 | | |
595 | 0 | let (mut sbet2, mut cbet2) = geomath::sincosd(lat2); |
596 | 0 | sbet2 *= self._f1; |
597 | | |
598 | 0 | geomath::norm(&mut sbet2, &mut cbet2); |
599 | 0 | cbet2 = cbet2.max(self.tiny_); |
600 | | |
601 | 0 | if cbet1 < -sbet1 { |
602 | 0 | if cbet2 == cbet1 { |
603 | 0 | sbet2 = if sbet2 < 0.0 { sbet1 } else { -sbet1 }; |
604 | 0 | } |
605 | 0 | } else if sbet2.abs() == -sbet1 { |
606 | 0 | cbet2 = cbet1; |
607 | 0 | } |
608 | | |
609 | 0 | let dn1 = (1.0 + self._ep2 * geomath::sq(sbet1)).sqrt(); |
610 | 0 | let dn2 = (1.0 + self._ep2 * geomath::sq(sbet2)).sqrt(); |
611 | | |
612 | 0 | let mut C1a: [f64; CARR_SIZE] = [0.0; CARR_SIZE]; |
613 | 0 | let mut C2a: [f64; CARR_SIZE] = [0.0; CARR_SIZE]; |
614 | 0 | let mut C3a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER]; |
615 | | |
616 | 0 | let mut meridian = lat1 == -90.0 || slam12 == 0.0; |
617 | 0 | let mut calp1 = 0.0; |
618 | 0 | let mut salp1 = 0.0; |
619 | 0 | let mut calp2 = 0.0; |
620 | 0 | let mut salp2 = 0.0; |
621 | 0 | let mut ssig1 = 0.0; |
622 | 0 | let mut csig1 = 0.0; |
623 | 0 | let mut ssig2 = 0.0; |
624 | 0 | let mut csig2 = 0.0; |
625 | | let mut sig12: f64; |
626 | 0 | let mut s12x = 0.0; |
627 | 0 | let mut m12x = 0.0; |
628 | | |
629 | 0 | if meridian { |
630 | 0 | calp1 = clam12; |
631 | 0 | salp1 = slam12; |
632 | 0 | calp2 = 1.0; |
633 | 0 | salp2 = 0.0; |
634 | | |
635 | 0 | ssig1 = sbet1; |
636 | 0 | csig1 = calp1 * cbet1; |
637 | 0 | ssig2 = sbet2; |
638 | 0 | csig2 = calp2 * cbet2; |
639 | | |
640 | 0 | sig12 = ((csig1 * ssig2 - ssig1 * csig2).max(0.0)).atan2(csig1 * csig2 + ssig1 * ssig2); |
641 | 0 | let res = self._Lengths( |
642 | 0 | self._n, |
643 | 0 | sig12, |
644 | 0 | ssig1, |
645 | 0 | csig1, |
646 | 0 | dn1, |
647 | 0 | ssig2, |
648 | 0 | csig2, |
649 | 0 | dn2, |
650 | 0 | cbet1, |
651 | 0 | cbet2, |
652 | 0 | outmask | caps::DISTANCE | caps::REDUCEDLENGTH, |
653 | 0 | &mut C1a, |
654 | 0 | &mut C2a, |
655 | | ); |
656 | 0 | s12x = res.0; |
657 | 0 | m12x = res.1; |
658 | 0 | M12 = res.3; |
659 | 0 | M21 = res.4; |
660 | | |
661 | 0 | if sig12 < 1.0 || m12x >= 0.0 { |
662 | 0 | if sig12 < 3.0 * self.tiny_ { |
663 | 0 | sig12 = 0.0; |
664 | 0 | m12x = 0.0; |
665 | 0 | s12x = 0.0; |
666 | 0 | } |
667 | 0 | m12x *= self._b; |
668 | 0 | s12x *= self._b; |
669 | 0 | a12 = sig12.to_degrees(); |
670 | 0 | } else { |
671 | 0 | meridian = false; |
672 | 0 | } |
673 | 0 | } |
674 | | |
675 | 0 | let mut somg12 = 2.0; |
676 | 0 | let mut comg12 = 0.0; |
677 | 0 | let mut omg12 = 0.0; |
678 | | let dnm: f64; |
679 | 0 | let mut eps = 0.0; |
680 | 0 | if !meridian && sbet1 == 0.0 && (self.f <= 0.0 || lon12s >= self.f * 180.0) { |
681 | 0 | calp1 = 0.0; |
682 | 0 | calp2 = 0.0; |
683 | 0 | salp1 = 1.0; |
684 | 0 | salp2 = 1.0; |
685 | | |
686 | 0 | s12x = self.a * lam12; |
687 | 0 | sig12 = lam12 / self._f1; |
688 | 0 | omg12 = lam12 / self._f1; |
689 | 0 | m12x = self._b * sig12.sin(); |
690 | 0 | if outmask & caps::GEODESICSCALE != 0 { |
691 | 0 | M12 = sig12.cos(); |
692 | 0 | M21 = sig12.cos(); |
693 | 0 | } |
694 | 0 | a12 = lon12 / self._f1; |
695 | 0 | } else if !meridian { |
696 | 0 | let res = self._InverseStart( |
697 | 0 | sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12, slam12, clam12, &mut C1a, &mut C2a, |
698 | | ); |
699 | 0 | sig12 = res.0; |
700 | 0 | salp1 = res.1; |
701 | 0 | calp1 = res.2; |
702 | 0 | salp2 = res.3; |
703 | 0 | calp2 = res.4; |
704 | 0 | dnm = res.5; |
705 | | |
706 | 0 | if sig12 >= 0.0 { |
707 | 0 | s12x = sig12 * self._b * dnm; |
708 | 0 | m12x = geomath::sq(dnm) * self._b * (sig12 / dnm).sin(); |
709 | 0 | if outmask & caps::GEODESICSCALE != 0 { |
710 | 0 | M12 = (sig12 / dnm).cos(); |
711 | 0 | M21 = (sig12 / dnm).cos(); |
712 | 0 | } |
713 | 0 | a12 = sig12.to_degrees(); |
714 | 0 | omg12 = lam12 / (self._f1 * dnm); |
715 | | } else { |
716 | 0 | let mut tripn = false; |
717 | 0 | let mut tripb = false; |
718 | 0 | let mut salp1a = self.tiny_; |
719 | 0 | let mut calp1a = 1.0; |
720 | 0 | let mut salp1b = self.tiny_; |
721 | 0 | let mut calp1b = -1.0; |
722 | 0 | let mut domg12 = 0.0; |
723 | 0 | for numit in 0..self.maxit2_ { |
724 | 0 | let res = self._Lambda12( |
725 | 0 | sbet1, |
726 | 0 | cbet1, |
727 | 0 | dn1, |
728 | 0 | sbet2, |
729 | 0 | cbet2, |
730 | 0 | dn2, |
731 | 0 | salp1, |
732 | 0 | calp1, |
733 | 0 | slam12, |
734 | 0 | clam12, |
735 | 0 | numit < self.maxit1_, |
736 | 0 | &mut C1a, |
737 | 0 | &mut C2a, |
738 | 0 | &mut C3a, |
739 | | ); |
740 | 0 | let v = res.0; |
741 | 0 | salp2 = res.1; |
742 | 0 | calp2 = res.2; |
743 | 0 | sig12 = res.3; |
744 | 0 | ssig1 = res.4; |
745 | 0 | csig1 = res.5; |
746 | 0 | ssig2 = res.6; |
747 | 0 | csig2 = res.7; |
748 | 0 | eps = res.8; |
749 | 0 | domg12 = res.9; |
750 | 0 | let dv = res.10; |
751 | | |
752 | 0 | if tripb |
753 | 0 | || v.abs() < if tripn { 8.0 } else { 1.0 } * self.tol0_ |
754 | 0 | || v.abs().is_nan() |
755 | | { |
756 | 0 | break; |
757 | 0 | }; |
758 | 0 | if v > 0.0 && (numit > self.maxit1_ || calp1 / salp1 > calp1b / salp1b) { |
759 | 0 | salp1b = salp1; |
760 | 0 | calp1b = calp1; |
761 | 0 | } else if v < 0.0 && (numit > self.maxit1_ || calp1 / salp1 < calp1a / salp1a) { |
762 | 0 | salp1a = salp1; |
763 | 0 | calp1a = calp1; |
764 | 0 | } |
765 | 0 | if numit < self.maxit1_ && dv > 0.0 { |
766 | 0 | let dalp1 = -v / dv; |
767 | 0 | let sdalp1 = dalp1.sin(); |
768 | 0 | let cdalp1 = dalp1.cos(); |
769 | 0 | let nsalp1 = salp1 * cdalp1 + calp1 * sdalp1; |
770 | 0 | if nsalp1 > 0.0 && dalp1.abs() < PI { |
771 | 0 | calp1 = calp1 * cdalp1 - salp1 * sdalp1; |
772 | 0 | salp1 = nsalp1; |
773 | 0 | geomath::norm(&mut salp1, &mut calp1); |
774 | 0 | tripn = v.abs() <= 16.0 * self.tol0_; |
775 | 0 | continue; |
776 | 0 | } |
777 | 0 | } |
778 | | |
779 | 0 | salp1 = (salp1a + salp1b) / 2.0; |
780 | 0 | calp1 = (calp1a + calp1b) / 2.0; |
781 | 0 | geomath::norm(&mut salp1, &mut calp1); |
782 | 0 | tripn = false; |
783 | 0 | tripb = (salp1a - salp1).abs() + (calp1a - calp1) < self.tolb_ |
784 | 0 | || (salp1 - salp1b).abs() + (calp1 - calp1b) < self.tolb_; |
785 | | } |
786 | 0 | let lengthmask = outmask |
787 | 0 | | if outmask & (caps::REDUCEDLENGTH | caps::GEODESICSCALE) != 0 { |
788 | 0 | caps::DISTANCE |
789 | | } else { |
790 | 0 | caps::EMPTY |
791 | | }; |
792 | 0 | let res = self._Lengths( |
793 | 0 | eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, lengthmask, |
794 | 0 | &mut C1a, &mut C2a, |
795 | | ); |
796 | 0 | s12x = res.0; |
797 | 0 | m12x = res.1; |
798 | 0 | M12 = res.3; |
799 | 0 | M21 = res.4; |
800 | | |
801 | 0 | m12x *= self._b; |
802 | 0 | s12x *= self._b; |
803 | 0 | a12 = sig12.to_degrees(); |
804 | 0 | if outmask & caps::AREA != 0 { |
805 | 0 | let sdomg12 = domg12.sin(); |
806 | 0 | let cdomg12 = domg12.cos(); |
807 | 0 | somg12 = slam12 * cdomg12 - clam12 * sdomg12; |
808 | 0 | comg12 = clam12 * cdomg12 + slam12 * sdomg12; |
809 | 0 | } |
810 | | } |
811 | 0 | } |
812 | 0 | if outmask & caps::DISTANCE != 0 { |
813 | 0 | s12 = 0.0 + s12x; |
814 | 0 | } |
815 | 0 | if outmask & caps::REDUCEDLENGTH != 0 { |
816 | 0 | m12 = 0.0 + m12x; |
817 | 0 | } |
818 | 0 | if outmask & caps::AREA != 0 { |
819 | 0 | let salp0 = salp1 * cbet1; |
820 | 0 | let calp0 = calp1.hypot(salp1 * sbet1); |
821 | 0 | if calp0 != 0.0 && salp0 != 0.0 { |
822 | 0 | ssig1 = sbet1; |
823 | 0 | csig1 = calp1 * cbet1; |
824 | 0 | ssig2 = sbet2; |
825 | 0 | csig2 = calp2 * cbet2; |
826 | 0 | let k2 = geomath::sq(calp0) * self._ep2; |
827 | 0 | eps = k2 / (2.0 * (1.0 + (1.0 + k2).sqrt()) + k2); |
828 | 0 | let A4 = geomath::sq(self.a) * calp0 * salp0 * self._e2; |
829 | 0 | geomath::norm(&mut ssig1, &mut csig1); |
830 | 0 | geomath::norm(&mut ssig2, &mut csig2); |
831 | 0 | let mut C4a: [f64; GEODESIC_ORDER] = [0.0; GEODESIC_ORDER]; |
832 | 0 | self._C4f(eps, &mut C4a); |
833 | 0 | let B41 = geomath::sin_cos_series(false, ssig1, csig1, &C4a); |
834 | 0 | let B42 = geomath::sin_cos_series(false, ssig2, csig2, &C4a); |
835 | 0 | S12 = A4 * (B42 - B41); |
836 | 0 | } else { |
837 | 0 | S12 = 0.0; |
838 | 0 | } |
839 | | |
840 | 0 | if !meridian && somg12 > 1.0 { |
841 | 0 | somg12 = omg12.sin(); |
842 | 0 | comg12 = omg12.cos(); |
843 | 0 | } |
844 | | |
845 | | // We're diverging from Karney's implementation here |
846 | | // which uses the hardcoded constant: -0.7071 for FRAC_1_SQRT_2 |
847 | | let alp12: f64; |
848 | 0 | if !meridian && comg12 > -FRAC_1_SQRT_2 && sbet2 - sbet1 < 1.75 { |
849 | 0 | let domg12 = 1.0 + comg12; |
850 | 0 | let dbet1 = 1.0 + cbet1; |
851 | 0 | let dbet2 = 1.0 + cbet2; |
852 | 0 | alp12 = 2.0 |
853 | 0 | * (somg12 * (sbet1 * dbet2 + sbet2 * dbet1)) |
854 | 0 | .atan2(domg12 * (sbet1 * sbet2 + dbet1 * dbet2)); |
855 | 0 | } else { |
856 | 0 | let mut salp12 = salp2 * calp1 - calp2 * salp1; |
857 | 0 | let mut calp12 = calp2 * calp1 + salp2 * salp1; |
858 | | |
859 | 0 | if salp12 == 0.0 && calp12 < 0.0 { |
860 | 0 | salp12 = self.tiny_ * calp1; |
861 | 0 | calp12 = -1.0; |
862 | 0 | } |
863 | 0 | alp12 = salp12.atan2(calp12); |
864 | | } |
865 | 0 | S12 += self._c2 * alp12; |
866 | 0 | S12 *= swapp * lonsign * latsign; |
867 | 0 | S12 += 0.0; |
868 | 0 | } |
869 | | |
870 | 0 | if swapp < 0.0 { |
871 | 0 | std::mem::swap(&mut salp2, &mut salp1); |
872 | | |
873 | 0 | std::mem::swap(&mut calp2, &mut calp1); |
874 | | |
875 | 0 | if outmask & caps::GEODESICSCALE != 0 { |
876 | 0 | std::mem::swap(&mut M21, &mut M12); |
877 | 0 | } |
878 | 0 | } |
879 | 0 | salp1 *= swapp * lonsign; |
880 | 0 | calp1 *= swapp * latsign; |
881 | 0 | salp2 *= swapp * lonsign; |
882 | 0 | calp2 *= swapp * latsign; |
883 | 0 | (a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12) |
884 | 0 | } |
885 | | |
886 | | /// returns (a12, lat2, lon2, azi2, s12, m12, M12, M21, S12) |
887 | 0 | pub fn _gen_direct( |
888 | 0 | &self, |
889 | 0 | lat1: f64, |
890 | 0 | lon1: f64, |
891 | 0 | azi1: f64, |
892 | 0 | arcmode: bool, |
893 | 0 | s12_a12: f64, |
894 | 0 | mut outmask: u64, |
895 | 0 | ) -> (f64, f64, f64, f64, f64, f64, f64, f64, f64) { |
896 | 0 | if !arcmode { |
897 | 0 | outmask |= caps::DISTANCE_IN |
898 | 0 | }; |
899 | | |
900 | 0 | let line = |
901 | 0 | geodesic_line::GeodesicLine::new(self, lat1, lon1, azi1, Some(outmask), None, None); |
902 | 0 | line._gen_position(arcmode, s12_a12, outmask) |
903 | 0 | } |
904 | | |
905 | | /// Get the area of the geodesic in square meters |
906 | 0 | pub fn area(&self) -> f64 { |
907 | 0 | self._c2 * 4.0 * std::f64::consts::PI |
908 | 0 | } |
909 | | } |
910 | | |
911 | | /// Place a second point, given the first point, an azimuth, and a distance. |
912 | | /// |
913 | | /// # Arguments |
914 | | /// - lat1 - Latitude of 1st point (degrees) [-90.,90.] |
915 | | /// - lon1 - Longitude of 1st point (degrees) [-180., 180.] |
916 | | /// - azi1 - Azimuth at 1st point (degrees) [-180., 180.] |
917 | | /// - s12 - Distance from 1st to 2nd point (meters) Value may be negative |
918 | | /// |
919 | | /// # Returns |
920 | | /// |
921 | | /// There are a variety of outputs associated with this calculation. We save computation by |
922 | | /// only calculating the outputs you need. See the following impls which return different subsets of |
923 | | /// the following outputs: |
924 | | /// |
925 | | /// - lat2 latitude of point 2 (degrees). |
926 | | /// - lon2 longitude of point 2 (degrees). |
927 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
928 | | /// - m12 reduced length of geodesic (meters). |
929 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
930 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
931 | | /// - S12 area under the geodesic (meters<sup>2</sup>). |
932 | | /// - a12 arc length between point 1 and point 2 (degrees). |
933 | | /// |
934 | | /// If either point is at a pole, the azimuth is defined by keeping the |
935 | | /// longitude fixed, writing lat = ±(90° − ε), and taking the limit ε → 0+. |
936 | | /// An arc length greater that 180° signifies a geodesic which is not a |
937 | | /// shortest path. (For a prolate ellipsoid, an additional condition is |
938 | | /// necessary for a shortest path: the longitudinal extent must not |
939 | | /// exceed of 180°.) |
940 | | /// ```rust |
941 | | /// // Example, determine the point 10000 km NE of JFK: |
942 | | /// use geographiclib_rs::{Geodesic, DirectGeodesic}; |
943 | | /// |
944 | | /// let g = Geodesic::wgs84(); |
945 | | /// let (lat, lon, az) = g.direct(40.64, -73.78, 45.0, 10e6); |
946 | | /// |
947 | | /// use approx::assert_relative_eq; |
948 | | /// assert_relative_eq!(lat, 32.621100463725796); |
949 | | /// assert_relative_eq!(lon, 49.052487092959836); |
950 | | /// assert_relative_eq!(az, 140.4059858768007); |
951 | | /// ``` |
952 | | pub trait DirectGeodesic<T> { |
953 | | fn direct(&self, lat1: f64, lon1: f64, azi1: f64, s12: f64) -> T; |
954 | | } |
955 | | |
956 | | impl DirectGeodesic<(f64, f64)> for Geodesic { |
957 | | /// See the documentation for the DirectGeodesic trait. |
958 | | /// |
959 | | /// # Returns |
960 | | /// - lat2 latitude of point 2 (degrees). |
961 | | /// - lon2 longitude of point 2 (degrees). |
962 | 0 | fn direct(&self, lat1: f64, lon1: f64, azi1: f64, s12: f64) -> (f64, f64) { |
963 | 0 | let capabilities = caps::LATITUDE | caps::LONGITUDE; |
964 | 0 | let (_a12, lat2, lon2, _azi2, _s12, _m12, _M12, _M21, _S12) = |
965 | 0 | self._gen_direct(lat1, lon1, azi1, false, s12, capabilities); |
966 | | |
967 | 0 | (lat2, lon2) |
968 | 0 | } |
969 | | } |
970 | | |
971 | | impl DirectGeodesic<(f64, f64, f64)> for Geodesic { |
972 | | /// See the documentation for the DirectGeodesic trait. |
973 | | /// |
974 | | /// # Returns |
975 | | /// - lat2 latitude of point 2 (degrees). |
976 | | /// - lon2 longitude of point 2 (degrees). |
977 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
978 | 0 | fn direct(&self, lat1: f64, lon1: f64, azi1: f64, s12: f64) -> (f64, f64, f64) { |
979 | 0 | let capabilities = caps::LATITUDE | caps::LONGITUDE | caps::AZIMUTH; |
980 | 0 | let (_a12, lat2, lon2, azi2, _s12, _m12, _M12, _M21, _S12) = |
981 | 0 | self._gen_direct(lat1, lon1, azi1, false, s12, capabilities); |
982 | | |
983 | 0 | (lat2, lon2, azi2) |
984 | 0 | } |
985 | | } |
986 | | |
987 | | impl DirectGeodesic<(f64, f64, f64, f64)> for Geodesic { |
988 | | /// See the documentation for the DirectGeodesic trait. |
989 | | /// |
990 | | /// # Returns |
991 | | /// - lat2 latitude of point 2 (degrees). |
992 | | /// - lon2 longitude of point 2 (degrees). |
993 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
994 | | /// - m12 reduced length of geodesic (meters). |
995 | 0 | fn direct(&self, lat1: f64, lon1: f64, azi1: f64, s12: f64) -> (f64, f64, f64, f64) { |
996 | 0 | let capabilities = caps::LATITUDE | caps::LONGITUDE | caps::AZIMUTH | caps::REDUCEDLENGTH; |
997 | 0 | let (_a12, lat2, lon2, azi2, _s12, m12, _M12, _M21, _S12) = |
998 | 0 | self._gen_direct(lat1, lon1, azi1, false, s12, capabilities); |
999 | | |
1000 | 0 | (lat2, lon2, azi2, m12) |
1001 | 0 | } |
1002 | | } |
1003 | | |
1004 | | impl DirectGeodesic<(f64, f64, f64, f64, f64)> for Geodesic { |
1005 | | /// See the documentation for the DirectGeodesic trait. |
1006 | | /// |
1007 | | /// # Returns |
1008 | | /// - lat2 latitude of point 2 (degrees). |
1009 | | /// - lon2 longitude of point 2 (degrees). |
1010 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1011 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
1012 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
1013 | 0 | fn direct(&self, lat1: f64, lon1: f64, azi1: f64, s12: f64) -> (f64, f64, f64, f64, f64) { |
1014 | 0 | let capabilities = caps::LATITUDE | caps::LONGITUDE | caps::AZIMUTH | caps::GEODESICSCALE; |
1015 | 0 | let (_a12, lat2, lon2, azi2, _s12, _m12, M12, M21, _S12) = |
1016 | 0 | self._gen_direct(lat1, lon1, azi1, false, s12, capabilities); |
1017 | | |
1018 | 0 | (lat2, lon2, azi2, M12, M21) |
1019 | 0 | } |
1020 | | } |
1021 | | |
1022 | | impl DirectGeodesic<(f64, f64, f64, f64, f64, f64)> for Geodesic { |
1023 | | /// See the documentation for the DirectGeodesic trait. |
1024 | | /// |
1025 | | /// # Returns |
1026 | | /// - lat2 latitude of point 2 (degrees). |
1027 | | /// - lon2 longitude of point 2 (degrees). |
1028 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1029 | | /// - m12 reduced length of geodesic (meters). |
1030 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
1031 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
1032 | 0 | fn direct(&self, lat1: f64, lon1: f64, azi1: f64, s12: f64) -> (f64, f64, f64, f64, f64, f64) { |
1033 | 0 | let capabilities = caps::LATITUDE |
1034 | 0 | | caps::LONGITUDE |
1035 | 0 | | caps::AZIMUTH |
1036 | 0 | | caps::REDUCEDLENGTH |
1037 | 0 | | caps::GEODESICSCALE; |
1038 | 0 | let (_a12, lat2, lon2, azi2, _s12, m12, M12, M21, _S12) = |
1039 | 0 | self._gen_direct(lat1, lon1, azi1, false, s12, capabilities); |
1040 | | |
1041 | 0 | (lat2, lon2, azi2, m12, M12, M21) |
1042 | 0 | } |
1043 | | } |
1044 | | |
1045 | | impl DirectGeodesic<(f64, f64, f64, f64, f64, f64, f64, f64)> for Geodesic { |
1046 | | /// See the documentation for the DirectGeodesic trait. |
1047 | | /// |
1048 | | /// # Returns |
1049 | | /// - lat2 latitude of point 2 (degrees). |
1050 | | /// - lon2 longitude of point 2 (degrees). |
1051 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1052 | | /// - m12 reduced length of geodesic (meters). |
1053 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
1054 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
1055 | | /// - S12 area under the geodesic (meters<sup>2</sup>). |
1056 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1057 | 0 | fn direct( |
1058 | 0 | &self, |
1059 | 0 | lat1: f64, |
1060 | 0 | lon1: f64, |
1061 | 0 | azi1: f64, |
1062 | 0 | s12: f64, |
1063 | 0 | ) -> (f64, f64, f64, f64, f64, f64, f64, f64) { |
1064 | 0 | let capabilities = caps::LATITUDE |
1065 | 0 | | caps::LONGITUDE |
1066 | 0 | | caps::AZIMUTH |
1067 | 0 | | caps::REDUCEDLENGTH |
1068 | 0 | | caps::GEODESICSCALE |
1069 | 0 | | caps::AREA; |
1070 | 0 | let (a12, lat2, lon2, azi2, _s12, m12, M12, M21, S12) = |
1071 | 0 | self._gen_direct(lat1, lon1, azi1, false, s12, capabilities); |
1072 | | |
1073 | 0 | (lat2, lon2, azi2, m12, M12, M21, S12, a12) |
1074 | 0 | } |
1075 | | } |
1076 | | |
1077 | | /// Measure the distance (and other values) between two points. |
1078 | | /// |
1079 | | /// # Arguments |
1080 | | /// - lat1 latitude of point 1 (degrees). |
1081 | | /// - lon1 longitude of point 1 (degrees). |
1082 | | /// - lat2 latitude of point 2 (degrees). |
1083 | | /// - lon2 longitude of point 2 (degrees). |
1084 | | /// |
1085 | | /// # Returns |
1086 | | /// |
1087 | | /// There are a variety of outputs associated with this calculation. We save computation by |
1088 | | /// only calculating the outputs you need. See the following impls which return different subsets of |
1089 | | /// the following outputs: |
1090 | | /// |
1091 | | /// - s12 distance between point 1 and point 2 (meters). |
1092 | | /// - azi1 azimuth at point 1 (degrees). |
1093 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1094 | | /// - m12 reduced length of geodesic (meters). |
1095 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
1096 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
1097 | | /// - S12 area under the geodesic (meters<sup>2</sup>). |
1098 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1099 | | /// |
1100 | | /// `lat1` and `lat2` should be in the range [−90°, 90°]. |
1101 | | /// The values of `azi1` and `azi2` returned are in the range |
1102 | | /// [−180°, 180°]. |
1103 | | /// |
1104 | | /// If either point is at a pole, the azimuth is defined by keeping the |
1105 | | /// longitude fixed, writing `lat` = ±(90° − ε), |
1106 | | /// and taking the limit ε → 0+. |
1107 | | /// |
1108 | | /// The solution to the inverse problem is found using Newton's method. If |
1109 | | /// this fails to converge (this is very unlikely in geodetic applications |
1110 | | /// but does occur for very eccentric ellipsoids), then the bisection method |
1111 | | /// is used to refine the solution. |
1112 | | /// |
1113 | | /// ```rust |
1114 | | /// // Example, determine the distance between two points |
1115 | | /// use geographiclib_rs::{Geodesic, InverseGeodesic}; |
1116 | | /// |
1117 | | /// let g = Geodesic::wgs84(); |
1118 | | /// let p1 = (34.095925, -118.2884237); |
1119 | | /// let p2 = (59.4323439, 24.7341649); |
1120 | | /// let s12: f64 = g.inverse(p1.0, p1.1, p2.0, p2.1); |
1121 | | /// |
1122 | | /// use approx::assert_relative_eq; |
1123 | | /// assert_relative_eq!(s12, 9094718.72751138); |
1124 | | /// ``` |
1125 | | pub trait InverseGeodesic<T> { |
1126 | | fn inverse(&self, lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> T; |
1127 | | } |
1128 | | |
1129 | | impl InverseGeodesic<f64> for Geodesic { |
1130 | | /// See the documentation for the InverseGeodesic trait. |
1131 | | /// |
1132 | | /// # Returns |
1133 | | /// - s12 distance between point 1 and point 2 (meters). |
1134 | 0 | fn inverse(&self, lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 { |
1135 | 0 | let capabilities = caps::DISTANCE; |
1136 | 0 | let (_a12, s12, _azi1, _azi2, _m12, _M12, _M21, _S12) = |
1137 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1138 | | |
1139 | 0 | s12 |
1140 | 0 | } |
1141 | | } |
1142 | | |
1143 | | impl InverseGeodesic<(f64, f64)> for Geodesic { |
1144 | | /// See the documentation for the InverseGeodesic trait. |
1145 | | /// |
1146 | | /// # Returns |
1147 | | /// - s12 distance between point 1 and point 2 (meters). |
1148 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1149 | 0 | fn inverse(&self, lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> (f64, f64) { |
1150 | 0 | let capabilities = caps::DISTANCE; |
1151 | 0 | let (a12, s12, _azi1, _azi2, _m12, _M12, _M21, _S12) = |
1152 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1153 | | |
1154 | 0 | (s12, a12) |
1155 | 0 | } |
1156 | | } |
1157 | | |
1158 | | impl InverseGeodesic<(f64, f64, f64)> for Geodesic { |
1159 | | /// See the documentation for the InverseGeodesic trait. |
1160 | | /// |
1161 | | /// # Returns |
1162 | | /// - azi1 azimuth at point 1 (degrees). |
1163 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1164 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1165 | 0 | fn inverse(&self, lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> (f64, f64, f64) { |
1166 | 0 | let capabilities = caps::AZIMUTH; |
1167 | 0 | let (a12, _s12, azi1, azi2, _m12, _M12, _M21, _S12) = |
1168 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1169 | | |
1170 | 0 | (azi1, azi2, a12) |
1171 | 0 | } |
1172 | | } |
1173 | | |
1174 | | impl InverseGeodesic<(f64, f64, f64, f64)> for Geodesic { |
1175 | | /// See the documentation for the InverseGeodesic trait. |
1176 | | /// |
1177 | | /// # Returns |
1178 | | /// - s12 distance between point 1 and point 2 (meters). |
1179 | | /// - azi1 azimuth at point 1 (degrees). |
1180 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1181 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1182 | 0 | fn inverse(&self, lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> (f64, f64, f64, f64) { |
1183 | 0 | let capabilities = caps::DISTANCE | caps::AZIMUTH; |
1184 | 0 | let (a12, s12, azi1, azi2, _m12, _M12, _M21, _S12) = |
1185 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1186 | | |
1187 | 0 | (s12, azi1, azi2, a12) |
1188 | 0 | } |
1189 | | } |
1190 | | |
1191 | | impl InverseGeodesic<(f64, f64, f64, f64, f64)> for Geodesic { |
1192 | | /// See the documentation for the InverseGeodesic trait. |
1193 | | /// |
1194 | | /// # Returns |
1195 | | /// - s12 distance between point 1 and point 2 (meters). |
1196 | | /// - azi1 azimuth at point 1 (degrees). |
1197 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1198 | | /// - m12 reduced length of geodesic (meters). |
1199 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1200 | 0 | fn inverse(&self, lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> (f64, f64, f64, f64, f64) { |
1201 | 0 | let capabilities = caps::DISTANCE | caps::AZIMUTH | caps::REDUCEDLENGTH; |
1202 | 0 | let (a12, s12, azi1, azi2, m12, _M12, _M21, _S12) = |
1203 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1204 | | |
1205 | 0 | (s12, azi1, azi2, m12, a12) |
1206 | 0 | } |
1207 | | } |
1208 | | |
1209 | | impl InverseGeodesic<(f64, f64, f64, f64, f64, f64)> for Geodesic { |
1210 | | /// See the documentation for the InverseGeodesic trait. |
1211 | | /// |
1212 | | /// # Returns |
1213 | | /// - s12 distance between point 1 and point 2 (meters). |
1214 | | /// - azi1 azimuth at point 1 (degrees). |
1215 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1216 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
1217 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
1218 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1219 | 0 | fn inverse( |
1220 | 0 | &self, |
1221 | 0 | lat1: f64, |
1222 | 0 | lon1: f64, |
1223 | 0 | lat2: f64, |
1224 | 0 | lon2: f64, |
1225 | 0 | ) -> (f64, f64, f64, f64, f64, f64) { |
1226 | 0 | let capabilities = caps::DISTANCE | caps::AZIMUTH | caps::GEODESICSCALE; |
1227 | 0 | let (a12, s12, azi1, azi2, _m12, M12, M21, _S12) = |
1228 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1229 | | |
1230 | 0 | (s12, azi1, azi2, M12, M21, a12) |
1231 | 0 | } |
1232 | | } |
1233 | | |
1234 | | impl InverseGeodesic<(f64, f64, f64, f64, f64, f64, f64)> for Geodesic { |
1235 | | /// See the documentation for the InverseGeodesic trait. |
1236 | | /// |
1237 | | /// # Returns |
1238 | | /// - s12 distance between point 1 and point 2 (meters). |
1239 | | /// - azi1 azimuth at point 1 (degrees). |
1240 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1241 | | /// - m12 reduced length of geodesic (meters). |
1242 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
1243 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
1244 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1245 | 0 | fn inverse( |
1246 | 0 | &self, |
1247 | 0 | lat1: f64, |
1248 | 0 | lon1: f64, |
1249 | 0 | lat2: f64, |
1250 | 0 | lon2: f64, |
1251 | 0 | ) -> (f64, f64, f64, f64, f64, f64, f64) { |
1252 | 0 | let capabilities = |
1253 | 0 | caps::DISTANCE | caps::AZIMUTH | caps::REDUCEDLENGTH | caps::GEODESICSCALE; |
1254 | 0 | let (a12, s12, azi1, azi2, m12, M12, M21, _S12) = |
1255 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1256 | | |
1257 | 0 | (s12, azi1, azi2, m12, M12, M21, a12) |
1258 | 0 | } |
1259 | | } |
1260 | | |
1261 | | impl InverseGeodesic<(f64, f64, f64, f64, f64, f64, f64, f64)> for Geodesic { |
1262 | | /// See the documentation for the InverseGeodesic trait. |
1263 | | /// |
1264 | | /// # Returns |
1265 | | /// - s12 distance between point 1 and point 2 (meters). |
1266 | | /// - azi1 azimuth at point 1 (degrees). |
1267 | | /// - azi2 (forward) azimuth at point 2 (degrees). |
1268 | | /// - m12 reduced length of geodesic (meters). |
1269 | | /// - M12 geodesic scale of point 2 relative to point 1 (dimensionless). |
1270 | | /// - M21 geodesic scale of point 1 relative to point 2 (dimensionless). |
1271 | | /// - S12 area under the geodesic (meters<sup>2</sup>). |
1272 | | /// - a12 arc length between point 1 and point 2 (degrees). |
1273 | 0 | fn inverse( |
1274 | 0 | &self, |
1275 | 0 | lat1: f64, |
1276 | 0 | lon1: f64, |
1277 | 0 | lat2: f64, |
1278 | 0 | lon2: f64, |
1279 | 0 | ) -> (f64, f64, f64, f64, f64, f64, f64, f64) { |
1280 | 0 | let capabilities = |
1281 | 0 | caps::DISTANCE | caps::AZIMUTH | caps::REDUCEDLENGTH | caps::GEODESICSCALE | caps::AREA; |
1282 | 0 | let (a12, s12, azi1, azi2, m12, M12, M21, S12) = |
1283 | 0 | self._gen_inverse_azi(lat1, lon1, lat2, lon2, capabilities); |
1284 | | |
1285 | 0 | (s12, azi1, azi2, m12, M12, M21, S12, a12) |
1286 | 0 | } |
1287 | | } |
1288 | | |
1289 | | #[cfg(test)] |
1290 | | mod tests { |
1291 | | use super::*; |
1292 | | use crate::geodesic_line::GeodesicLine; |
1293 | | use approx::assert_relative_eq; |
1294 | | use std::io::BufRead; |
1295 | | |
1296 | | #[allow(clippy::type_complexity)] |
1297 | | const TESTCASES: &[(f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64)] = &[ |
1298 | | ( |
1299 | | 35.60777, |
1300 | | -139.44815, |
1301 | | 111.098748429560326, |
1302 | | -11.17491, |
1303 | | -69.95921, |
1304 | | 129.289270889708762, |
1305 | | 8935244.5604818305, |
1306 | | 80.50729714281974, |
1307 | | 6273170.2055303837, |
1308 | | 0.16606318447386067, |
1309 | | 0.16479116945612937, |
1310 | | 12841384694976.432, |
1311 | | ), |
1312 | | ( |
1313 | | 55.52454, |
1314 | | 106.05087, |
1315 | | 22.020059880982801, |
1316 | | 77.03196, |
1317 | | 197.18234, |
1318 | | 109.112041110671519, |
1319 | | 4105086.1713924406, |
1320 | | 36.892740690445894, |
1321 | | 3828869.3344387607, |
1322 | | 0.80076349608092607, |
1323 | | 0.80101006984201008, |
1324 | | 61674961290615.615, |
1325 | | ), |
1326 | | ( |
1327 | | -21.97856, |
1328 | | 142.59065, |
1329 | | -32.44456876433189, |
1330 | | 41.84138, |
1331 | | 98.56635, |
1332 | | -41.84359951440466, |
1333 | | 8394328.894657671, |
1334 | | 75.62930491011522, |
1335 | | 6161154.5773110616, |
1336 | | 0.24816339233950381, |
1337 | | 0.24930251203627892, |
1338 | | -6637997720646.717, |
1339 | | ), |
1340 | | ( |
1341 | | -66.99028, |
1342 | | 112.2363, |
1343 | | 173.73491240878403, |
1344 | | -12.70631, |
1345 | | 285.90344, |
1346 | | 2.512956620913668, |
1347 | | 11150344.2312080241, |
1348 | | 100.278634181155759, |
1349 | | 6289939.5670446687, |
1350 | | -0.17199490274700385, |
1351 | | -0.17722569526345708, |
1352 | | -121287239862139.744, |
1353 | | ), |
1354 | | ( |
1355 | | -17.42761, |
1356 | | 173.34268, |
1357 | | -159.033557661192928, |
1358 | | -15.84784, |
1359 | | 5.93557, |
1360 | | -20.787484651536988, |
1361 | | 16076603.1631180673, |
1362 | | 144.640108810286253, |
1363 | | 3732902.1583877189, |
1364 | | -0.81273638700070476, |
1365 | | -0.81299800519154474, |
1366 | | 97825992354058.708, |
1367 | | ), |
1368 | | ( |
1369 | | 32.84994, |
1370 | | 48.28919, |
1371 | | 150.492927788121982, |
1372 | | -56.28556, |
1373 | | 202.29132, |
1374 | | 48.113449399816759, |
1375 | | 16727068.9438164461, |
1376 | | 150.565799985466607, |
1377 | | 3147838.1910180939, |
1378 | | -0.87334918086923126, |
1379 | | -0.86505036767110637, |
1380 | | -72445258525585.010, |
1381 | | ), |
1382 | | ( |
1383 | | 6.96833, |
1384 | | 52.74123, |
1385 | | 92.581585386317712, |
1386 | | -7.39675, |
1387 | | 206.17291, |
1388 | | 90.721692165923907, |
1389 | | 17102477.2496958388, |
1390 | | 154.147366239113561, |
1391 | | 2772035.6169917581, |
1392 | | -0.89991282520302447, |
1393 | | -0.89986892177110739, |
1394 | | -1311796973197.995, |
1395 | | ), |
1396 | | ( |
1397 | | -50.56724, |
1398 | | -16.30485, |
1399 | | -105.439679907590164, |
1400 | | -33.56571, |
1401 | | -94.97412, |
1402 | | -47.348547835650331, |
1403 | | 6455670.5118668696, |
1404 | | 58.083719495371259, |
1405 | | 5409150.7979815838, |
1406 | | 0.53053508035997263, |
1407 | | 0.52988722644436602, |
1408 | | 41071447902810.047, |
1409 | | ), |
1410 | | ( |
1411 | | -58.93002, |
1412 | | -8.90775, |
1413 | | 140.965397902500679, |
1414 | | -8.91104, |
1415 | | 133.13503, |
1416 | | 19.255429433416599, |
1417 | | 11756066.0219864627, |
1418 | | 105.755691241406877, |
1419 | | 6151101.2270708536, |
1420 | | -0.26548622269867183, |
1421 | | -0.27068483874510741, |
1422 | | -86143460552774.735, |
1423 | | ), |
1424 | | ( |
1425 | | -68.82867, |
1426 | | -74.28391, |
1427 | | 93.774347763114881, |
1428 | | -50.63005, |
1429 | | -8.36685, |
1430 | | 34.65564085411343, |
1431 | | 3956936.926063544, |
1432 | | 35.572254987389284, |
1433 | | 3708890.9544062657, |
1434 | | 0.81443963736383502, |
1435 | | 0.81420859815358342, |
1436 | | -41845309450093.787, |
1437 | | ), |
1438 | | ( |
1439 | | -10.62672, |
1440 | | -32.0898, |
1441 | | -86.426713286747751, |
1442 | | 5.883, |
1443 | | -134.31681, |
1444 | | -80.473780971034875, |
1445 | | 11470869.3864563009, |
1446 | | 103.387395634504061, |
1447 | | 6184411.6622659713, |
1448 | | -0.23138683500430237, |
1449 | | -0.23155097622286792, |
1450 | | 4198803992123.548, |
1451 | | ), |
1452 | | ( |
1453 | | -21.76221, |
1454 | | 166.90563, |
1455 | | 29.319421206936428, |
1456 | | 48.72884, |
1457 | | 213.97627, |
1458 | | 43.508671946410168, |
1459 | | 9098627.3986554915, |
1460 | | 81.963476716121964, |
1461 | | 6299240.9166992283, |
1462 | | 0.13965943368590333, |
1463 | | 0.14152969707656796, |
1464 | | 10024709850277.476, |
1465 | | ), |
1466 | | ( |
1467 | | -19.79938, |
1468 | | -174.47484, |
1469 | | 71.167275780171533, |
1470 | | -11.99349, |
1471 | | -154.35109, |
1472 | | 65.589099775199228, |
1473 | | 2319004.8601169389, |
1474 | | 20.896611684802389, |
1475 | | 2267960.8703918325, |
1476 | | 0.93427001867125849, |
1477 | | 0.93424887135032789, |
1478 | | -3935477535005.785, |
1479 | | ), |
1480 | | ( |
1481 | | -11.95887, |
1482 | | -116.94513, |
1483 | | 92.712619830452549, |
1484 | | 4.57352, |
1485 | | 7.16501, |
1486 | | 78.64960934409585, |
1487 | | 13834722.5801401374, |
1488 | | 124.688684161089762, |
1489 | | 5228093.177931598, |
1490 | | -0.56879356755666463, |
1491 | | -0.56918731952397221, |
1492 | | -9919582785894.853, |
1493 | | ), |
1494 | | ( |
1495 | | -87.85331, |
1496 | | 85.66836, |
1497 | | -65.120313040242748, |
1498 | | 66.48646, |
1499 | | 16.09921, |
1500 | | -4.888658719272296, |
1501 | | 17286615.3147144645, |
1502 | | 155.58592449699137, |
1503 | | 2635887.4729110181, |
1504 | | -0.90697975771398578, |
1505 | | -0.91095608883042767, |
1506 | | 42667211366919.534, |
1507 | | ), |
1508 | | ( |
1509 | | 1.74708, |
1510 | | 128.32011, |
1511 | | -101.584843631173858, |
1512 | | -11.16617, |
1513 | | 11.87109, |
1514 | | -86.325793296437476, |
1515 | | 12942901.1241347408, |
1516 | | 116.650512484301857, |
1517 | | 5682744.8413270572, |
1518 | | -0.44857868222697644, |
1519 | | -0.44824490340007729, |
1520 | | 10763055294345.653, |
1521 | | ), |
1522 | | ( |
1523 | | -25.72959, |
1524 | | -144.90758, |
1525 | | -153.647468693117198, |
1526 | | -57.70581, |
1527 | | -269.17879, |
1528 | | -48.343983158876487, |
1529 | | 9413446.7452453107, |
1530 | | 84.664533838404295, |
1531 | | 6356176.6898881281, |
1532 | | 0.09492245755254703, |
1533 | | 0.09737058264766572, |
1534 | | 74515122850712.444, |
1535 | | ), |
1536 | | ( |
1537 | | -41.22777, |
1538 | | 122.32875, |
1539 | | 14.285113402275739, |
1540 | | -7.57291, |
1541 | | 130.37946, |
1542 | | 10.805303085187369, |
1543 | | 3812686.035106021, |
1544 | | 34.34330804743883, |
1545 | | 3588703.8812128856, |
1546 | | 0.82605222593217889, |
1547 | | 0.82572158200920196, |
1548 | | -2456961531057.857, |
1549 | | ), |
1550 | | ( |
1551 | | 11.01307, |
1552 | | 138.25278, |
1553 | | 79.43682622782374, |
1554 | | 6.62726, |
1555 | | 247.05981, |
1556 | | 103.708090215522657, |
1557 | | 11911190.819018408, |
1558 | | 107.341669954114577, |
1559 | | 6070904.722786735, |
1560 | | -0.29767608923657404, |
1561 | | -0.29785143390252321, |
1562 | | 17121631423099.696, |
1563 | | ), |
1564 | | ( |
1565 | | -29.47124, |
1566 | | 95.14681, |
1567 | | -163.779130441688382, |
1568 | | -27.46601, |
1569 | | -69.15955, |
1570 | | -15.909335945554969, |
1571 | | 13487015.8381145492, |
1572 | | 121.294026715742277, |
1573 | | 5481428.9945736388, |
1574 | | -0.51527225545373252, |
1575 | | -0.51556587964721788, |
1576 | | 104679964020340.318, |
1577 | | ), |
1578 | | ]; |
1579 | | |
1580 | | #[test] |
1581 | | fn test_inverse_and_direct() -> Result<(), String> { |
1582 | | // See python/test_geodesic.py |
1583 | | let geod = Geodesic::wgs84(); |
1584 | | let (_a12, s12, _azi1, _azi2, _m12, _M12, _M21, _S12) = |
1585 | | geod._gen_inverse_azi(0.0, 0.0, 1.0, 1.0, caps::STANDARD); |
1586 | | assert_eq!(s12, 156899.56829134026); |
1587 | | |
1588 | | // Test inverse |
1589 | | for (lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12) in TESTCASES.iter() { |
1590 | | let ( |
1591 | | computed_a12, |
1592 | | computed_s12, |
1593 | | computed_azi1, |
1594 | | computed_azi2, |
1595 | | computed_m12, |
1596 | | computed_M12, |
1597 | | computed_M21, |
1598 | | computed_S12, |
1599 | | ) = geod._gen_inverse_azi(*lat1, *lon1, *lat2, *lon2, caps::ALL | caps::LONG_UNROLL); |
1600 | | assert_relative_eq!(computed_azi1, azi1, epsilon = 1e-13f64); |
1601 | | assert_relative_eq!(computed_azi2, azi2, epsilon = 1e-13f64); |
1602 | | assert_relative_eq!(computed_s12, s12, epsilon = 1e-8f64); |
1603 | | assert_relative_eq!(computed_a12, a12, epsilon = 1e-13f64); |
1604 | | assert_relative_eq!(computed_m12, m12, epsilon = 1e-8f64); |
1605 | | assert_relative_eq!(computed_M12, M12, epsilon = 1e-15f64); |
1606 | | assert_relative_eq!(computed_M21, M21, epsilon = 1e-15f64); |
1607 | | assert_relative_eq!(computed_S12, S12, epsilon = 0.1f64); |
1608 | | } |
1609 | | |
1610 | | // Test direct |
1611 | | for (lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12) in TESTCASES.iter() { |
1612 | | let ( |
1613 | | computed_a12, |
1614 | | computed_lat2, |
1615 | | computed_lon2, |
1616 | | computed_azi2, |
1617 | | _computed_s12, |
1618 | | computed_m12, |
1619 | | computed_M12, |
1620 | | computed_M21, |
1621 | | computed_S12, |
1622 | | ) = geod._gen_direct( |
1623 | | *lat1, |
1624 | | *lon1, |
1625 | | *azi1, |
1626 | | false, |
1627 | | *s12, |
1628 | | caps::ALL | caps::LONG_UNROLL, |
1629 | | ); |
1630 | | assert_relative_eq!(computed_lat2, lat2, epsilon = 1e-13f64); |
1631 | | assert_relative_eq!(computed_lon2, lon2, epsilon = 1e-13f64); |
1632 | | assert_relative_eq!(computed_azi2, azi2, epsilon = 1e-13f64); |
1633 | | assert_relative_eq!(computed_a12, a12, epsilon = 1e-13f64); |
1634 | | assert_relative_eq!(computed_m12, m12, epsilon = 1e-8f64); |
1635 | | assert_relative_eq!(computed_M12, M12, epsilon = 1e-15f64); |
1636 | | assert_relative_eq!(computed_M21, M21, epsilon = 1e-15f64); |
1637 | | assert_relative_eq!(computed_S12, S12, epsilon = 0.1f64); |
1638 | | } |
1639 | | Ok(()) |
1640 | | } |
1641 | | |
1642 | | #[test] |
1643 | | fn test_arcdirect() { |
1644 | | // Corresponds with ArcDirectCheck from Java, or test_arcdirect from Python |
1645 | | let geod = Geodesic::wgs84(); |
1646 | | for (lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12) in TESTCASES.iter() { |
1647 | | let ( |
1648 | | _computed_a12, |
1649 | | computed_lat2, |
1650 | | computed_lon2, |
1651 | | computed_azi2, |
1652 | | computed_s12, |
1653 | | computed_m12, |
1654 | | computed_M12, |
1655 | | computed_M21, |
1656 | | computed_S12, |
1657 | | ) = geod._gen_direct( |
1658 | | *lat1, |
1659 | | *lon1, |
1660 | | *azi1, |
1661 | | true, |
1662 | | *a12, |
1663 | | caps::ALL | caps::LONG_UNROLL, |
1664 | | ); |
1665 | | assert_relative_eq!(computed_lat2, lat2, epsilon = 1e-13); |
1666 | | assert_relative_eq!(computed_lon2, lon2, epsilon = 1e-13); |
1667 | | assert_relative_eq!(computed_azi2, azi2, epsilon = 1e-13); |
1668 | | assert_relative_eq!(computed_s12, s12, epsilon = 1e-8); |
1669 | | assert_relative_eq!(computed_m12, m12, epsilon = 1e-8); |
1670 | | assert_relative_eq!(computed_M12, M12, epsilon = 1e-15); |
1671 | | assert_relative_eq!(computed_M21, M21, epsilon = 1e-15); |
1672 | | assert_relative_eq!(computed_S12, S12, epsilon = 0.1); |
1673 | | } |
1674 | | } |
1675 | | |
1676 | | #[test] |
1677 | | fn test_geninverse() { |
1678 | | let geod = Geodesic::wgs84(); |
1679 | | let res = geod._gen_inverse(0.0, 0.0, 1.0, 1.0, caps::STANDARD); |
1680 | | assert_eq!(res.0, 1.4141938478710363); |
1681 | | assert_eq!(res.1, 156899.56829134026); |
1682 | | assert_eq!(res.2, 0.7094236375834774); |
1683 | | assert_eq!(res.3, 0.7047823085448635); |
1684 | | assert_eq!(res.4, 0.7095309793242709); |
1685 | | assert_eq!(res.5, 0.7046742434480923); |
1686 | | assert!(res.6.is_nan()); |
1687 | | assert!(res.7.is_nan()); |
1688 | | assert!(res.8.is_nan()); |
1689 | | assert!(res.9.is_nan()); |
1690 | | } |
1691 | | |
1692 | | #[test] |
1693 | | fn test_inverse_start() { |
1694 | | let geod = Geodesic::wgs84(); |
1695 | | let res = geod._InverseStart( |
1696 | | -0.017393909556108908, |
1697 | | 0.9998487145115275, |
1698 | | 1.0000010195104125, |
1699 | | -0.0, |
1700 | | 1.0, |
1701 | | 1.0, |
1702 | | 0.017453292519943295, |
1703 | | 0.01745240643728351, |
1704 | | 0.9998476951563913, |
1705 | | &mut [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], |
1706 | | &mut [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], |
1707 | | ); |
1708 | | assert_eq!(res.0, -1.0); |
1709 | | assert_relative_eq!(res.1, 0.7095310092765433, epsilon = 1e-13); |
1710 | | assert_relative_eq!(res.2, 0.7046742132893822, epsilon = 1e-13); |
1711 | | assert!(res.3.is_nan()); |
1712 | | assert!(res.4.is_nan()); |
1713 | | assert_eq!(res.5, 1.0000002548969817); |
1714 | | |
1715 | | let res = geod._InverseStart( |
1716 | | -0.017393909556108908, |
1717 | | 0.9998487145115275, |
1718 | | 1.0000010195104125, |
1719 | | -0.0, |
1720 | | 1.0, |
1721 | | 1.0, |
1722 | | 0.017453292519943295, |
1723 | | 0.01745240643728351, |
1724 | | 0.9998476951563913, |
1725 | | &mut [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], |
1726 | | &mut [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], |
1727 | | ); |
1728 | | assert_eq!(res.0, -1.0); |
1729 | | assert_relative_eq!(res.1, 0.7095310092765433, epsilon = 1e-13); |
1730 | | assert_relative_eq!(res.2, 0.7046742132893822, epsilon = 1e-13); |
1731 | | assert!(res.3.is_nan()); |
1732 | | assert!(res.4.is_nan()); |
1733 | | assert_eq!(res.5, 1.0000002548969817); |
1734 | | } |
1735 | | |
1736 | | #[test] |
1737 | | fn test_lambda12() { |
1738 | | let geod = Geodesic::wgs84(); |
1739 | | let res1 = geod._Lambda12( |
1740 | | -0.017393909556108908, |
1741 | | 0.9998487145115275, |
1742 | | 1.0000010195104125, |
1743 | | -0.0, |
1744 | | 1.0, |
1745 | | 1.0, |
1746 | | 0.7095310092765433, |
1747 | | 0.7046742132893822, |
1748 | | 0.01745240643728351, |
1749 | | 0.9998476951563913, |
1750 | | true, |
1751 | | &mut [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], |
1752 | | &mut [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], |
1753 | | &mut [0.0, 1.0, 2.0, 3.0, 4.0, 5.0], |
1754 | | ); |
1755 | | assert_eq!(res1.0, 1.4834408705897495e-09); |
1756 | | assert_eq!(res1.1, 0.7094236675312185); |
1757 | | assert_eq!(res1.2, 0.7047822783999007); |
1758 | | assert_eq!(res1.3, 0.024682339962725352); |
1759 | | assert_eq!(res1.4, -0.024679833885152578); |
1760 | | assert_eq!(res1.5, 0.9996954065111039); |
1761 | | assert_eq!(res1.6, -0.0); |
1762 | | assert_eq!(res1.7, 1.0); |
1763 | | assert_relative_eq!(res1.8, 0.0008355095326524276, epsilon = 1e-13); |
1764 | | assert_eq!(res1.9, -5.8708496511415445e-05); |
1765 | | assert_eq!(res1.10, 0.034900275148485); |
1766 | | |
1767 | | let res2 = geod._Lambda12( |
1768 | | -0.017393909556108908, |
1769 | | 0.9998487145115275, |
1770 | | 1.0000010195104125, |
1771 | | -0.0, |
1772 | | 1.0, |
1773 | | 1.0, |
1774 | | 0.7095309793242709, |
1775 | | 0.7046742434480923, |
1776 | | 0.01745240643728351, |
1777 | | 0.9998476951563913, |
1778 | | true, |
1779 | | &mut [ |
1780 | | 0.0, |
1781 | | -0.00041775465696698233, |
1782 | | -4.362974596862037e-08, |
1783 | | -1.2151022357848552e-11, |
1784 | | -4.7588881620421004e-15, |
1785 | | -2.226614930167366e-18, |
1786 | | -1.1627237498131586e-21, |
1787 | | ], |
1788 | | &mut [ |
1789 | | 0.0, |
1790 | | -0.0008355098973052918, |
1791 | | -1.7444619952659748e-07, |
1792 | | -7.286557795511902e-11, |
1793 | | -3.80472772706481e-14, |
1794 | | -2.2251271876594078e-17, |
1795 | | 1.2789961247944744e-20, |
1796 | | ], |
1797 | | &mut [ |
1798 | | 0.0, |
1799 | | 0.00020861391868413911, |
1800 | | 4.3547247296823945e-08, |
1801 | | 1.515432276542012e-11, |
1802 | | 6.645637323698485e-15, |
1803 | | 3.3399223952510497e-18, |
1804 | | ], |
1805 | | ); |
1806 | | assert_eq!(res2.0, 6.046459990680098e-17); |
1807 | | assert_eq!(res2.1, 0.7094236375834774); |
1808 | | assert_eq!(res2.2, 0.7047823085448635); |
1809 | | assert_eq!(res2.3, 0.024682338906797385); |
1810 | | assert_eq!(res2.4, -0.02467983282954624); |
1811 | | assert_eq!(res2.5, 0.9996954065371639); |
1812 | | assert_eq!(res2.6, -0.0); |
1813 | | assert_eq!(res2.7, 1.0); |
1814 | | assert_relative_eq!(res2.8, 0.0008355096040059597, epsilon = 1e-18); |
1815 | | assert_eq!(res2.9, -5.870849152149326e-05); |
1816 | | assert_eq!(res2.10, 0.03490027216297455); |
1817 | | } |
1818 | | |
1819 | | #[test] |
1820 | | fn test_lengths() { |
1821 | | // Results taken from the python implementation |
1822 | | let geod = Geodesic::wgs84(); |
1823 | | let mut c1a = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]; |
1824 | | let mut c2a = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]; |
1825 | | let res1 = geod._Lengths( |
1826 | | 0.0008355095326524276, |
1827 | | 0.024682339962725352, |
1828 | | -0.024679833885152578, |
1829 | | 0.9996954065111039, |
1830 | | 1.0000010195104125, |
1831 | | -0.0, |
1832 | | 1.0, |
1833 | | 1.0, |
1834 | | 0.9998487145115275, |
1835 | | 1.0, |
1836 | | 4101, |
1837 | | &mut c1a, |
1838 | | &mut c2a, |
1839 | | ); |
1840 | | assert!(res1.0.is_nan()); |
1841 | | assert_eq!(res1.1, 0.024679842274314294); |
1842 | | assert_eq!(res1.2, 0.0016717180169067588); |
1843 | | assert!(res1.3.is_nan()); |
1844 | | assert!(res1.4.is_nan()); |
1845 | | |
1846 | | let res2 = geod._Lengths( |
1847 | | 0.0008355096040059597, |
1848 | | 0.024682338906797385, |
1849 | | -0.02467983282954624, |
1850 | | 0.9996954065371639, |
1851 | | 1.0000010195104125, |
1852 | | -0.0, |
1853 | | 1.0, |
1854 | | 1.0, |
1855 | | 0.9998487145115275, |
1856 | | 1.0, |
1857 | | 4101, |
1858 | | &mut [ |
1859 | | 0.0, |
1860 | | -0.00041775465696698233, |
1861 | | -4.362974596862037e-08, |
1862 | | -1.2151022357848552e-11, |
1863 | | -4.7588881620421004e-15, |
1864 | | -2.226614930167366e-18, |
1865 | | -1.1627237498131586e-21, |
1866 | | ], |
1867 | | &mut [ |
1868 | | 0.0, |
1869 | | -0.0008355098973052918, |
1870 | | -1.7444619952659748e-07, |
1871 | | -7.286557795511902e-11, |
1872 | | -3.80472772706481e-14, |
1873 | | -2.2251271876594078e-17, |
1874 | | 1.2789961247944744e-20, |
1875 | | ], |
1876 | | ); |
1877 | | assert!(res2.0.is_nan()); |
1878 | | assert_eq!(res2.1, 0.02467984121870759); |
1879 | | assert_eq!(res2.2, 0.0016717181597332804); |
1880 | | assert!(res2.3.is_nan()); |
1881 | | assert!(res2.4.is_nan()); |
1882 | | |
1883 | | let res3 = geod._Lengths( |
1884 | | 0.0008355096040059597, |
1885 | | 0.024682338906797385, |
1886 | | -0.02467983282954624, |
1887 | | 0.9996954065371639, |
1888 | | 1.0000010195104125, |
1889 | | -0.0, |
1890 | | 1.0, |
1891 | | 1.0, |
1892 | | 0.9998487145115275, |
1893 | | 1.0, |
1894 | | 1920, |
1895 | | &mut [ |
1896 | | 0.0, |
1897 | | -0.00041775469264372037, |
1898 | | -4.362975342068502e-08, |
1899 | | -1.215102547098435e-11, |
1900 | | -4.758889787701359e-15, |
1901 | | -2.2266158809456692e-18, |
1902 | | -1.1627243456014359e-21, |
1903 | | ], |
1904 | | &mut [ |
1905 | | 0.0, |
1906 | | -0.0008355099686589174, |
1907 | | -1.744462293162189e-07, |
1908 | | -7.286559662008413e-11, |
1909 | | -3.804729026574989e-14, |
1910 | | -2.2251281376754273e-17, |
1911 | | 1.2789967801615795e-20, |
1912 | | ], |
1913 | | ); |
1914 | | assert_eq!(res3.0, 0.024682347295447677); |
1915 | | assert!(res3.1.is_nan()); |
1916 | | assert!(res3.2.is_nan()); |
1917 | | assert!(res3.3.is_nan()); |
1918 | | assert!(res3.4.is_nan()); |
1919 | | |
1920 | | let res = geod._Lengths( |
1921 | | 0.0007122620325664751, |
1922 | | 1.405117407023628, |
1923 | | -0.8928657853278468, |
1924 | | 0.45032287238256896, |
1925 | | 1.0011366173804046, |
1926 | | 0.2969032234925426, |
1927 | | 0.9549075745221299, |
1928 | | 1.0001257451360057, |
1929 | | 0.8139459053827204, |
1930 | | 0.9811634781422108, |
1931 | | 1920, |
1932 | | &mut [ |
1933 | | 0.0, |
1934 | | -0.0003561309485314716, |
1935 | | -3.170731714689771e-08, |
1936 | | -7.527972480734327e-12, |
1937 | | -2.5133854116682488e-15, |
1938 | | -1.0025061462383107e-18, |
1939 | | -4.462794158625518e-22, |
1940 | | ], |
1941 | | &mut [ |
1942 | | 0.0, |
1943 | | -0.0007122622584701569, |
1944 | | -1.2678416507678478e-07, |
1945 | | -4.514641118748122e-11, |
1946 | | -2.0096353119518367e-14, |
1947 | | -1.0019350865558619e-17, |
1948 | | 4.90907357448807e-21, |
1949 | | ], |
1950 | | ); |
1951 | | assert_eq!(res.0, 1.4056304412645388); |
1952 | | assert!(res.1.is_nan()); |
1953 | | assert!(res.2.is_nan()); |
1954 | | assert!(res.3.is_nan()); |
1955 | | assert!(res.4.is_nan()); |
1956 | | } |
1957 | | |
1958 | | #[test] |
1959 | | fn test_goed__C4f() { |
1960 | | let geod = Geodesic::wgs84(); |
1961 | | let mut c = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]; |
1962 | | geod._C4f(0.12, &mut c); |
1963 | | assert_eq!( |
1964 | | c, |
1965 | | [ |
1966 | | 0.6420952961066771, |
1967 | | 0.0023680700061156517, |
1968 | | 9.96704067834604e-05, |
1969 | | 5.778187189466089e-06, |
1970 | | 3.9979026199316593e-07, |
1971 | | 3.2140078103714466e-08, |
1972 | | ] |
1973 | | ); |
1974 | | } |
1975 | | |
1976 | | #[test] |
1977 | | fn test_goed__C3f() { |
1978 | | let geod = Geodesic::wgs84(); |
1979 | | let mut c = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]; |
1980 | | geod._C3f(0.12, &mut c); |
1981 | | assert_eq!( |
1982 | | c, |
1983 | | [ |
1984 | | 1.0, |
1985 | | 0.031839442894193756, |
1986 | | 0.0009839921354137713, |
1987 | | 5.0055242248766214e-05, |
1988 | | 3.1656788204092044e-06, |
1989 | | 2.0412e-07, |
1990 | | ] |
1991 | | ); |
1992 | | } |
1993 | | |
1994 | | #[test] |
1995 | | fn test_goed__A3f() { |
1996 | | let geod = Geodesic::wgs84(); |
1997 | | assert_eq!(geod._A3f(0.12), 0.9363788874000158); |
1998 | | } |
1999 | | |
2000 | | #[test] |
2001 | | fn test_geod_init() { |
2002 | | // Check that after the init the variables are correctly set. |
2003 | | // Actual values are taken from the python implementation |
2004 | | let geod = Geodesic::wgs84(); |
2005 | | assert_eq!(geod.a, 6378137.0, "geod.a wrong"); |
2006 | | assert_eq!(geod.f, 0.0033528106647474805, "geod.f wrong"); |
2007 | | assert_eq!(geod._f1, 0.9966471893352525, "geod._f1 wrong"); |
2008 | | assert_eq!(geod._e2, 0.0066943799901413165, "geod._e2 wrong"); |
2009 | | assert_eq!(geod._ep2, 0.006739496742276434, "geod._ep2 wrong"); |
2010 | | assert_eq!(geod._n, 0.0016792203863837047, "geod._n wrong"); |
2011 | | assert_eq!(geod._b, 6356752.314245179, "geod._b wrong"); |
2012 | | assert_eq!(geod._c2, 40589732499314.76, "geod._c2 wrong"); |
2013 | | assert_eq!(geod._etol2, 3.6424611488788524e-08, "geod._etol2 wrong"); |
2014 | | assert_eq!( |
2015 | | geod._A3x, |
2016 | | [ |
2017 | | -0.0234375, |
2018 | | -0.046927475637074494, |
2019 | | -0.06281503005876607, |
2020 | | -0.2502088451303832, |
2021 | | -0.49916038980680816, |
2022 | | 1.0 |
2023 | | ], |
2024 | | "geod._A3x wrong" |
2025 | | ); |
2026 | | |
2027 | | assert_eq!( |
2028 | | geod._C3x, |
2029 | | [ |
2030 | | 0.0234375, |
2031 | | 0.03908873781853724, |
2032 | | 0.04695366939653196, |
2033 | | 0.12499964752736174, |
2034 | | 0.24958019490340408, |
2035 | | 0.01953125, |
2036 | | 0.02345061890926862, |
2037 | | 0.046822392185686165, |
2038 | | 0.062342661206936094, |
2039 | | 0.013671875, |
2040 | | 0.023393770302437927, |
2041 | | 0.025963026642854565, |
2042 | | 0.013671875, |
2043 | | 0.01362595881755982, |
2044 | | 0.008203125 |
2045 | | ], |
2046 | | "geod._C3x wrong" |
2047 | | ); |
2048 | | assert_eq!( |
2049 | | geod._C4x, |
2050 | | [ |
2051 | | 0.00646020646020646, |
2052 | | 0.0035037627212872787, |
2053 | | 0.034742279454780166, |
2054 | | -0.01921732223244865, |
2055 | | -0.19923321555984239, |
2056 | | 0.6662190894642603, |
2057 | | 0.000111000111000111, |
2058 | | 0.003426620602971002, |
2059 | | -0.009510765372597735, |
2060 | | -0.01893413691235592, |
2061 | | 0.0221370239510936, |
2062 | | 0.0007459207459207459, |
2063 | | -0.004142006291321442, |
2064 | | -0.00504225176309005, |
2065 | | 0.007584982177746079, |
2066 | | -0.0021565735851450138, |
2067 | | -0.001962613370670692, |
2068 | | 0.0036104265913438913, |
2069 | | -0.0009472009472009472, |
2070 | | 0.0020416649913317735, |
2071 | | 0.0012916376552740189 |
2072 | | ], |
2073 | | "geod._C4x wrong" |
2074 | | ); |
2075 | | } |
2076 | | |
2077 | | // The test_std_geodesic_* tests below are based on Karney's GeodSolve unit |
2078 | | // tests, found in many geographiclib variants. |
2079 | | // The versions below are mostly adapted from their Java counterparts, |
2080 | | // which use a testing structure more similar to Rust than do the C++ versions. |
2081 | | // Note that the Java tests often incorporate more than one of the C++ tests, |
2082 | | // and take their name from the lowest-numbered test in the set. |
2083 | | // These tests use that convention as well. |
2084 | | |
2085 | | #[test] |
2086 | | fn test_std_geodesic_geodsolve0() { |
2087 | | let geod = Geodesic::wgs84(); |
2088 | | let (s12, azi1, azi2, _a12) = geod.inverse(40.6, -73.8, 49.01666667, 2.55); |
2089 | | assert_relative_eq!(azi1, 53.47022, epsilon = 0.5e-5); |
2090 | | assert_relative_eq!(azi2, 111.59367, epsilon = 0.5e-5); |
2091 | | assert_relative_eq!(s12, 5853226.0, epsilon = 0.5); |
2092 | | } |
2093 | | |
2094 | | #[test] |
2095 | | fn test_std_geodesic_geodsolve1() { |
2096 | | let geod = Geodesic::wgs84(); |
2097 | | let (lat2, lon2, azi2) = geod.direct(40.63972222, -73.77888889, 53.5, 5850e3); |
2098 | | assert_relative_eq!(lat2, 49.01467, epsilon = 0.5e-5); |
2099 | | assert_relative_eq!(lon2, 2.56106, epsilon = 0.5e-5); |
2100 | | assert_relative_eq!(azi2, 111.62947, epsilon = 0.5e-5); |
2101 | | } |
2102 | | |
2103 | | #[test] |
2104 | | fn test_std_geodesic_geodsolve2() { |
2105 | | // Check fix for antipodal prolate bug found 2010-09-04 |
2106 | | let geod = Geodesic::new(6.4e6, -1f64 / 150.0); |
2107 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.07476, 0.0, -0.07476, 180.0); |
2108 | | assert_relative_eq!(azi1, 90.00078, epsilon = 0.5e-5); |
2109 | | assert_relative_eq!(azi2, 90.00078, epsilon = 0.5e-5); |
2110 | | assert_relative_eq!(s12, 20106193.0, epsilon = 0.5); |
2111 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.1, 0.0, -0.1, 180.0); |
2112 | | assert_relative_eq!(azi1, 90.00105, epsilon = 0.5e-5); |
2113 | | assert_relative_eq!(azi2, 90.00105, epsilon = 0.5e-5); |
2114 | | assert_relative_eq!(s12, 20106193.0, epsilon = 0.5); |
2115 | | } |
2116 | | |
2117 | | #[test] |
2118 | | fn test_std_geodesic_geodsolve4() { |
2119 | | // Check fix for short line bug found 2010-05-21 |
2120 | | let geod = Geodesic::wgs84(); |
2121 | | let s12: f64 = geod.inverse(36.493349428792, 0.0, 36.49334942879201, 0.0000008); |
2122 | | assert_relative_eq!(s12, 0.072, epsilon = 0.5e-3); |
2123 | | } |
2124 | | |
2125 | | #[test] |
2126 | | fn test_std_geodesic_geodsolve5() { |
2127 | | // Check fix for point2=pole bug found 2010-05-03 |
2128 | | let geod = Geodesic::wgs84(); |
2129 | | let (lat2, lon2, azi2) = geod.direct(0.01777745589997, 30.0, 0.0, 10e6); |
2130 | | assert_relative_eq!(lat2, 90.0, epsilon = 0.5e-5); |
2131 | | if lon2 < 0.0 { |
2132 | | assert_relative_eq!(lon2, -150.0, epsilon = 0.5e-5); |
2133 | | assert_relative_eq!(azi2.abs(), 180.0, epsilon = 0.5e-5); |
2134 | | } else { |
2135 | | assert_relative_eq!(lon2, 30.0, epsilon = 0.5e-5); |
2136 | | assert_relative_eq!(azi2, 0.0, epsilon = 0.5e-5); |
2137 | | } |
2138 | | } |
2139 | | |
2140 | | #[test] |
2141 | | fn test_std_geodesic_geodsolve6() { |
2142 | | // Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4.4 |
2143 | | // x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6.1). |
2144 | | let geod = Geodesic::wgs84(); |
2145 | | let s12: f64 = geod.inverse( |
2146 | | 88.202499451857, |
2147 | | 0.0, |
2148 | | -88.202499451857, |
2149 | | 179.981022032992859592, |
2150 | | ); |
2151 | | assert_relative_eq!(s12, 20003898.214, epsilon = 0.5e-3); |
2152 | | let s12: f64 = geod.inverse( |
2153 | | 89.333123580033, |
2154 | | 0.0, |
2155 | | -89.333123580032997687, |
2156 | | 179.99295812360148422, |
2157 | | ); |
2158 | | assert_relative_eq!(s12, 20003926.881, epsilon = 0.5e-3); |
2159 | | } |
2160 | | |
2161 | | #[test] |
2162 | | fn test_std_geodesic_geodsolve9() { |
2163 | | // Check fix for volatile x bug found 2011-06-25 (gcc 4.4.4 x86 -O3) |
2164 | | let geod = Geodesic::wgs84(); |
2165 | | let s12: f64 = geod.inverse( |
2166 | | 56.320923501171, |
2167 | | 0.0, |
2168 | | -56.320923501171, |
2169 | | 179.664747671772880215, |
2170 | | ); |
2171 | | assert_relative_eq!(s12, 19993558.287, epsilon = 0.5e-3); |
2172 | | } |
2173 | | |
2174 | | #[test] |
2175 | | fn test_std_geodesic_geodsolve10() { |
2176 | | // Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio |
2177 | | // 10 rel + debug) |
2178 | | let geod = Geodesic::wgs84(); |
2179 | | let s12: f64 = geod.inverse( |
2180 | | 52.784459512564, |
2181 | | 0.0, |
2182 | | -52.784459512563990912, |
2183 | | 179.634407464943777557, |
2184 | | ); |
2185 | | assert_relative_eq!(s12, 19991596.095, epsilon = 0.5e-3); |
2186 | | } |
2187 | | |
2188 | | #[test] |
2189 | | fn test_std_geodesic_geodsolve11() { |
2190 | | // Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio |
2191 | | // 10 rel + debug) |
2192 | | let geod = Geodesic::wgs84(); |
2193 | | let s12: f64 = geod.inverse( |
2194 | | 48.522876735459, |
2195 | | 0.0, |
2196 | | -48.52287673545898293, |
2197 | | 179.599720456223079643, |
2198 | | ); |
2199 | | assert_relative_eq!(s12, 19989144.774, epsilon = 0.5e-3); |
2200 | | } |
2201 | | |
2202 | | #[test] |
2203 | | fn test_std_geodesic_geodsolve12() { |
2204 | | // Check fix for inverse geodesics on extreme prolate/oblate |
2205 | | // ellipsoids Reported 2012-08-29 Stefan Guenther |
2206 | | // <stefan.gunther@embl.de>; fixed 2012-10-07 |
2207 | | let geod = Geodesic::new(89.8, -1.83); |
2208 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, -10.0, 160.0); |
2209 | | assert_relative_eq!(azi1, 120.27, epsilon = 1e-2); |
2210 | | assert_relative_eq!(azi2, 105.15, epsilon = 1e-2); |
2211 | | assert_relative_eq!(s12, 266.7, epsilon = 1e-1); |
2212 | | } |
2213 | | |
2214 | | #[test] |
2215 | | fn test_std_geodesic_geodsolve14() { |
2216 | | // Check fix for inverse ignoring lon12 = nan |
2217 | | let geod = Geodesic::wgs84(); |
2218 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 1.0, f64::NAN); |
2219 | | assert!(azi1.is_nan()); |
2220 | | assert!(azi2.is_nan()); |
2221 | | assert!(s12.is_nan()); |
2222 | | } |
2223 | | |
2224 | | #[test] |
2225 | | fn test_std_geodesic_geodsolve15() { |
2226 | | // Initial implementation of Math::eatanhe was wrong for e^2 < 0. This |
2227 | | // checks that this is fixed. |
2228 | | let geod = Geodesic::new(6.4e6, -1f64 / 150.0); |
2229 | | let (_lat2, _lon2, _azi2, _m12, _M12, _M21, S12, _a12) = geod.direct(1.0, 2.0, 3.0, 4.0); |
2230 | | assert_relative_eq!(S12, 23700.0, epsilon = 0.5); |
2231 | | } |
2232 | | |
2233 | | #[test] |
2234 | | fn test_std_geodesic_geodsolve17() { |
2235 | | // Check fix for LONG_UNROLL bug found on 2015-05-07 |
2236 | | let geod = Geodesic::new(6.4e6, -1f64 / 150.0); |
2237 | | let (_a12, lat2, lon2, azi2, _s12, _m12, _M12, _M21, _S12) = geod._gen_direct( |
2238 | | 40.0, |
2239 | | -75.0, |
2240 | | -10.0, |
2241 | | false, |
2242 | | 2e7, |
2243 | | caps::STANDARD | caps::LONG_UNROLL, |
2244 | | ); |
2245 | | assert_relative_eq!(lat2, -39.0, epsilon = 1.0); |
2246 | | assert_relative_eq!(lon2, -254.0, epsilon = 1.0); |
2247 | | assert_relative_eq!(azi2, -170.0, epsilon = 1.0); |
2248 | | |
2249 | | let line = GeodesicLine::new(&geod, 40.0, -75.0, -10.0, None, None, None); |
2250 | | let (_a12, lat2, lon2, azi2, _s12, _m12, _M12, _M21, _S12) = |
2251 | | line._gen_position(false, 2e7, caps::STANDARD | caps::LONG_UNROLL); |
2252 | | assert_relative_eq!(lat2, -39.0, epsilon = 1.0); |
2253 | | assert_relative_eq!(lon2, -254.0, epsilon = 1.0); |
2254 | | assert_relative_eq!(azi2, -170.0, epsilon = 1.0); |
2255 | | |
2256 | | let (lat2, lon2, azi2) = geod.direct(40.0, -75.0, -10.0, 2e7); |
2257 | | assert_relative_eq!(lat2, -39.0, epsilon = 1.0); |
2258 | | assert_relative_eq!(lon2, 105.0, epsilon = 1.0); |
2259 | | assert_relative_eq!(azi2, -170.0, epsilon = 1.0); |
2260 | | |
2261 | | let (_a12, lat2, lon2, azi2, _s12, _m12, _M12, _M21, _S12) = |
2262 | | line._gen_position(false, 2e7, caps::STANDARD); |
2263 | | assert_relative_eq!(lat2, -39.0, epsilon = 1.0); |
2264 | | assert_relative_eq!(lon2, 105.0, epsilon = 1.0); |
2265 | | assert_relative_eq!(azi2, -170.0, epsilon = 1.0); |
2266 | | } |
2267 | | |
2268 | | #[test] |
2269 | | fn test_std_geodesic_geodsolve26() { |
2270 | | // Check 0/0 problem with area calculation on sphere 2015-09-08 |
2271 | | let geod = Geodesic::new(6.4e6, 0.0); |
2272 | | let (_a12, _s12, _salp1, _calp1, _salp2, _calp2, _m12, _M12, _M21, S12) = |
2273 | | geod._gen_inverse(1.0, 2.0, 3.0, 4.0, caps::AREA); |
2274 | | assert_relative_eq!(S12, 49911046115.0, epsilon = 0.5); |
2275 | | } |
2276 | | |
2277 | | #[test] |
2278 | | fn test_std_geodesic_geodsolve28() { |
2279 | | // Check for bad placement of assignment of r.a12 with |f| > 0.01 (bug in |
2280 | | // Java implementation fixed on 2015-05-19). |
2281 | | let geod = Geodesic::new(6.4e6, 0.1); |
2282 | | let (a12, _lat2, _lon2, _azi2, _s12, _m12, _M12, _M21, _S12) = |
2283 | | geod._gen_direct(1.0, 2.0, 10.0, false, 5e6, caps::STANDARD); |
2284 | | assert_relative_eq!(a12, 48.55570690, epsilon = 0.5e-8); |
2285 | | } |
2286 | | |
2287 | | #[test] |
2288 | | fn test_std_geodesic_geodsolve29() { |
2289 | | // Check longitude unrolling with inverse calculation 2015-09-16 |
2290 | | let geod = Geodesic::wgs84(); |
2291 | | let (_a12, s12, _salp1, _calp1, _salp2, _calp2, _m12, _M12, _M21, _S12) = |
2292 | | geod._gen_inverse(0.0, 539.0, 0.0, 181.0, caps::STANDARD); |
2293 | | // Note: This is also supposed to check adjusted longitudes, but geographiclib-rs |
2294 | | // doesn't seem to support that as of 2021/01/18. |
2295 | | // assert_relative_eq!(lon1, 179, epsilon = 1e-10); |
2296 | | // assert_relative_eq!(lon2, -179, epsilon = 1e-10); |
2297 | | assert_relative_eq!(s12, 222639.0, epsilon = 0.5); |
2298 | | let (_a12, s12, _salp1, _calp1, _salp2, _calp2, _m12, _M12, _M21, _S12) = |
2299 | | geod._gen_inverse(0.0, 539.0, 0.0, 181.0, caps::STANDARD | caps::LONG_UNROLL); |
2300 | | // assert_relative_eq!(lon1, 539, epsilon = 1e-10); |
2301 | | // assert_relative_eq!(lon2, 541, epsilon = 1e-10); |
2302 | | assert_relative_eq!(s12, 222639.0, epsilon = 0.5); |
2303 | | } |
2304 | | |
2305 | | #[test] |
2306 | | fn test_std_geodesic_geodsolve33() { |
2307 | | // Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave -- |
2308 | | // sind(-0.0) = +0.0 -- and in some version of Visual Studio -- |
2309 | | // fmod(-0.0, 360.0) = +0.0. |
2310 | | let geod = Geodesic::wgs84(); |
2311 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.0, 179.0); |
2312 | | assert_relative_eq!(azi1, 90.0, epsilon = 0.5e-5); |
2313 | | assert_relative_eq!(azi2, 90.0, epsilon = 0.5e-5); |
2314 | | assert_relative_eq!(s12, 19926189.0, epsilon = 0.5); |
2315 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.0, 179.5); |
2316 | | assert_relative_eq!(azi1, 55.96650, epsilon = 0.5e-5); |
2317 | | assert_relative_eq!(azi2, 124.03350, epsilon = 0.5e-5); |
2318 | | assert_relative_eq!(s12, 19980862.0, epsilon = 0.5); |
2319 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.0, 180.0); |
2320 | | assert_relative_eq!(azi1, 0.0, epsilon = 0.5e-5); |
2321 | | assert_relative_eq!(azi2.abs(), 180.0, epsilon = 0.5e-5); |
2322 | | assert_relative_eq!(s12, 20003931.0, epsilon = 0.5); |
2323 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 1.0, 180.0); |
2324 | | assert_relative_eq!(azi1, 0.0, epsilon = 0.5e-5); |
2325 | | assert_relative_eq!(azi2.abs(), 180.0, epsilon = 0.5e-5); |
2326 | | assert_relative_eq!(s12, 19893357.0, epsilon = 0.5); |
2327 | | |
2328 | | let geod = Geodesic::new(6.4e6, 0.0); |
2329 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.0, 179.0); |
2330 | | assert_relative_eq!(azi1, 90.0, epsilon = 0.5e-5); |
2331 | | assert_relative_eq!(azi2, 90.0, epsilon = 0.5e-5); |
2332 | | assert_relative_eq!(s12, 19994492.0, epsilon = 0.5); |
2333 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.0, 180.0); |
2334 | | assert_relative_eq!(azi1, 0.0, epsilon = 0.5e-5); |
2335 | | assert_relative_eq!(azi2.abs(), 180.0, epsilon = 0.5e-5); |
2336 | | assert_relative_eq!(s12, 20106193.0, epsilon = 0.5); |
2337 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 1.0, 180.0); |
2338 | | assert_relative_eq!(azi1, 0.0, epsilon = 0.5e-5); |
2339 | | assert_relative_eq!(azi2.abs(), 180.0, epsilon = 0.5e-5); |
2340 | | assert_relative_eq!(s12, 19994492.0, epsilon = 0.5); |
2341 | | |
2342 | | let geod = Geodesic::new(6.4e6, -1.0 / 300.0); |
2343 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.0, 179.0); |
2344 | | assert_relative_eq!(azi1, 90.0, epsilon = 0.5e-5); |
2345 | | assert_relative_eq!(azi2, 90.0, epsilon = 0.5e-5); |
2346 | | assert_relative_eq!(s12, 19994492.0, epsilon = 0.5); |
2347 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.0, 180.0); |
2348 | | assert_relative_eq!(azi1, 90.0, epsilon = 0.5e-5); |
2349 | | assert_relative_eq!(azi2, 90.0, epsilon = 0.5e-5); |
2350 | | assert_relative_eq!(s12, 20106193.0, epsilon = 0.5); |
2351 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 0.5, 180.0); |
2352 | | assert_relative_eq!(azi1, 33.02493, epsilon = 0.5e-5); |
2353 | | assert_relative_eq!(azi2, 146.97364, epsilon = 0.5e-5); |
2354 | | assert_relative_eq!(s12, 20082617.0, epsilon = 0.5); |
2355 | | let (s12, azi1, azi2, _a12) = geod.inverse(0.0, 0.0, 1.0, 180.0); |
2356 | | assert_relative_eq!(azi1, 0.0, epsilon = 0.5e-5); |
2357 | | assert_relative_eq!(azi2.abs(), 180.0, epsilon = 0.5e-5); |
2358 | | assert_relative_eq!(s12, 20027270.0, epsilon = 0.5); |
2359 | | } |
2360 | | |
2361 | | #[test] |
2362 | | fn test_std_geodesic_geodsolve55() { |
2363 | | // Check fix for nan + point on equator or pole not returning all nans in |
2364 | | // Geodesic::Inverse, found 2015-09-23. |
2365 | | let geod = Geodesic::wgs84(); |
2366 | | let (s12, azi1, azi2, _a12) = geod.inverse(f64::NAN, 0.0, 0.0, 90.0); |
2367 | | assert!(azi1.is_nan()); |
2368 | | assert!(azi2.is_nan()); |
2369 | | assert!(s12.is_nan()); |
2370 | | let (s12, azi1, azi2, _a12) = geod.inverse(f64::NAN, 0.0, 90.0, 3.0); |
2371 | | assert!(azi1.is_nan()); |
2372 | | assert!(azi2.is_nan()); |
2373 | | assert!(s12.is_nan()); |
2374 | | } |
2375 | | |
2376 | | #[test] |
2377 | | fn test_std_geodesic_geodsolve59() { |
2378 | | // Check for points close with longitudes close to 180 deg apart. |
2379 | | let geod = Geodesic::wgs84(); |
2380 | | let (s12, azi1, azi2, _a12) = geod.inverse(5.0, 0.00000000000001, 10.0, 180.0); |
2381 | | assert_relative_eq!(azi1, 0.000000000000035, epsilon = 1.5e-14); |
2382 | | assert_relative_eq!(azi2, 179.99999999999996, epsilon = 1.5e-14); |
2383 | | assert_relative_eq!(s12, 18345191.174332713, epsilon = 5e-9); |
2384 | | } |
2385 | | |
2386 | | #[test] |
2387 | | fn test_std_geodesic_geodsolve61() { |
2388 | | // Make sure small negative azimuths are west-going |
2389 | | let geod = Geodesic::wgs84(); |
2390 | | let (_a12, lat2, lon2, azi2, _s12, _m12, _M12, _M21, _S12) = geod._gen_direct( |
2391 | | 45.0, |
2392 | | 0.0, |
2393 | | -0.000000000000000003, |
2394 | | false, |
2395 | | 1e7, |
2396 | | caps::STANDARD | caps::LONG_UNROLL, |
2397 | | ); |
2398 | | assert_relative_eq!(lat2, 45.30632, epsilon = 0.5e-5); |
2399 | | assert_relative_eq!(lon2, -180.0, epsilon = 0.5e-5); |
2400 | | assert_relative_eq!(azi2.abs(), 180.0, epsilon = 0.5e-5); |
2401 | | // geographiclib-rs does not appear to support Geodesic.inverse_line or |
2402 | | // or GeodesicLine.position as of 2021/01/18. |
2403 | | // let line = geod.inverse_line(45, 0, 80, -0.000000000000000003); |
2404 | | // let res = line.position(1e7, caps::STANDARD | caps::LONG_UNROLL); |
2405 | | // assert_relative_eq!(lat2, 45.30632, epsilon = 0.5e-5); |
2406 | | // assert_relative_eq!(lon2, -180, epsilon = 0.5e-5); |
2407 | | // assert_relative_eq!(azi2.abs(), 180, epsilon = 0.5e-5); |
2408 | | } |
2409 | | |
2410 | | // #[test] |
2411 | | // fn test_std_geodesic_geodsolve65() { |
2412 | | // // Check for bug in east-going check in GeodesicLine (needed to check for |
2413 | | // // sign of 0) and sign error in area calculation due to a bogus override |
2414 | | // // of the code for alp12. Found/fixed on 2015-12-19. |
2415 | | // // These tests rely on Geodesic.inverse_line, which is not supported by |
2416 | | // // geographiclib-rs as of 2021/01/18. |
2417 | | // } |
2418 | | |
2419 | | // #[test] |
2420 | | // fn test_std_geodesic_geodsolve69() { |
2421 | | // // Check for InverseLine if line is slightly west of S and that s13 is |
2422 | | // // correctly set. |
2423 | | // // These tests rely on Geodesic.inverse_line, which is not supported by |
2424 | | // // geographiclib-rs as of 2021/01/18. |
2425 | | // } |
2426 | | |
2427 | | // #[test] |
2428 | | // fn test_std_geodesic_geodsolve71() { |
2429 | | // // Check that DirectLine sets s13. |
2430 | | // // These tests rely on Geodesic.direct_line, which is not supported by |
2431 | | // // geographiclib-rs as of 2021/01/18. |
2432 | | // } |
2433 | | |
2434 | | #[test] |
2435 | | fn test_std_geodesic_geodsolve73() { |
2436 | | // Check for backwards from the pole bug reported by Anon on 2016-02-13. |
2437 | | // This only affected the Java implementation. It was introduced in Java |
2438 | | // version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17. |
2439 | | // Also the + sign on azi2 is a check on the normalizing of azimuths |
2440 | | // (converting -0.0 to +0.0). |
2441 | | let geod = Geodesic::wgs84(); |
2442 | | let (lat2, lon2, azi2) = geod.direct(90.0, 10.0, 180.0, -1e6); |
2443 | | assert_relative_eq!(lat2, 81.04623, epsilon = 0.5e-5); |
2444 | | assert_relative_eq!(lon2, -170.0, epsilon = 0.5e-5); |
2445 | | assert_relative_eq!(azi2, 0.0, epsilon = 0.5e-5); |
2446 | | assert!(azi2.is_sign_positive()); |
2447 | | } |
2448 | | |
2449 | | #[test] |
2450 | | fn test_std_geodesic_geodsolve74() { |
2451 | | // Check fix for inaccurate areas, bug introduced in v1.46, fixed |
2452 | | // 2015-10-16. |
2453 | | let geod = Geodesic::wgs84(); |
2454 | | let (a12, s12, azi1, azi2, m12, M12, M21, S12) = |
2455 | | geod._gen_inverse_azi(54.1589, 15.3872, 54.1591, 15.3877, caps::ALL); |
2456 | | assert_relative_eq!(azi1, 55.723110355, epsilon = 5e-9); |
2457 | | assert_relative_eq!(azi2, 55.723515675, epsilon = 5e-9); |
2458 | | assert_relative_eq!(s12, 39.527686385, epsilon = 5e-9); |
2459 | | assert_relative_eq!(a12, 0.000355495, epsilon = 5e-9); |
2460 | | assert_relative_eq!(m12, 39.527686385, epsilon = 5e-9); |
2461 | | assert_relative_eq!(M12, 0.999999995, epsilon = 5e-9); |
2462 | | assert_relative_eq!(M21, 0.999999995, epsilon = 5e-9); |
2463 | | assert_relative_eq!(S12, 286698586.30197, epsilon = 5e-4); |
2464 | | } |
2465 | | |
2466 | | #[test] |
2467 | | fn test_std_geodesic_geodsolve76() { |
2468 | | // The distance from Wellington and Salamanca (a classic failure of |
2469 | | // Vincenty) |
2470 | | let geod = Geodesic::wgs84(); |
2471 | | let (s12, azi1, azi2, _a12) = geod.inverse( |
2472 | | -(41.0 + 19.0 / 60.0), |
2473 | | 174.0 + 49.0 / 60.0, |
2474 | | 40.0 + 58.0 / 60.0, |
2475 | | -(5.0 + 30.0 / 60.0), |
2476 | | ); |
2477 | | assert_relative_eq!(azi1, 160.39137649664, epsilon = 0.5e-11); |
2478 | | assert_relative_eq!(azi2, 19.50042925176, epsilon = 0.5e-11); |
2479 | | assert_relative_eq!(s12, 19960543.857179, epsilon = 0.5e-6); |
2480 | | } |
2481 | | |
2482 | | #[test] |
2483 | | fn test_std_geodesic_geodsolve78() { |
2484 | | // An example where the NGS calculator fails to converge |
2485 | | let geod = Geodesic::wgs84(); |
2486 | | let (s12, azi1, azi2, _a12) = geod.inverse(27.2, 0.0, -27.1, 179.5); |
2487 | | assert_relative_eq!(azi1, 45.82468716758, epsilon = 0.5e-11); |
2488 | | assert_relative_eq!(azi2, 134.22776532670, epsilon = 0.5e-11); |
2489 | | assert_relative_eq!(s12, 19974354.765767, epsilon = 0.5e-6); |
2490 | | } |
2491 | | |
2492 | | #[test] |
2493 | | fn test_std_geodesic_geodsolve80() { |
2494 | | // Some tests to add code coverage: computing scale in special cases + zero |
2495 | | // length geodesic (includes GeodSolve80 - GeodSolve83). |
2496 | | let geod = Geodesic::wgs84(); |
2497 | | let (_a12, _s12, _salp1, _calp1, _salp2, _calp2, _m12, M12, M21, _S12) = |
2498 | | geod._gen_inverse(0.0, 0.0, 0.0, 90.0, caps::GEODESICSCALE); |
2499 | | assert_relative_eq!(M12, -0.00528427534, epsilon = 0.5e-10); |
2500 | | assert_relative_eq!(M21, -0.00528427534, epsilon = 0.5e-10); |
2501 | | |
2502 | | let (_a12, _s12, _salp1, _calp1, _salp2, _calp2, _m12, M12, M21, _S12) = |
2503 | | geod._gen_inverse(0.0, 0.0, 1e-6, 1e-6, caps::GEODESICSCALE); |
2504 | | assert_relative_eq!(M12, 1.0, epsilon = 0.5e-10); |
2505 | | assert_relative_eq!(M21, 1.0, epsilon = 0.5e-10); |
2506 | | |
2507 | | let (a12, s12, azi1, azi2, m12, M12, M21, S12) = |
2508 | | geod._gen_inverse_azi(20.001, 0.0, 20.001, 0.0, caps::ALL); |
2509 | | assert_relative_eq!(a12, 0.0, epsilon = 1e-13); |
2510 | | assert_relative_eq!(s12, 0.0, epsilon = 1e-8); |
2511 | | assert_relative_eq!(azi1, 180.0, epsilon = 1e-13); |
2512 | | assert_relative_eq!(azi2, 180.0, epsilon = 1e-13); |
2513 | | assert_relative_eq!(m12, 0.0, epsilon = 1e-8); |
2514 | | assert_relative_eq!(M12, 1.0, epsilon = 1e-15); |
2515 | | assert_relative_eq!(M21, 1.0, epsilon = 1e-15); |
2516 | | assert_relative_eq!(S12, 0.0, epsilon = 1e-10); |
2517 | | |
2518 | | let (a12, s12, azi1, azi2, m12, M12, M21, S12) = |
2519 | | geod._gen_inverse_azi(90.0, 0.0, 90.0, 180.0, caps::ALL); |
2520 | | assert_relative_eq!(a12, 0.0, epsilon = 1e-13); |
2521 | | assert_relative_eq!(s12, 0.0, epsilon = 1e-8); |
2522 | | assert_relative_eq!(azi1, 0.0, epsilon = 1e-13); |
2523 | | assert_relative_eq!(azi2, 180.0, epsilon = 1e-13); |
2524 | | assert_relative_eq!(m12, 0.0, epsilon = 1e-8); |
2525 | | assert_relative_eq!(M12, 1.0, epsilon = 1e-15); |
2526 | | assert_relative_eq!(M21, 1.0, epsilon = 1e-15); |
2527 | | assert_relative_eq!(S12, 127516405431022.0, epsilon = 0.5); |
2528 | | |
2529 | | // An incapable line which can't take distance as input |
2530 | | let line = GeodesicLine::new(&geod, 1.0, 2.0, 90.0, Some(caps::LATITUDE), None, None); |
2531 | | let (a12, _lat2, _lon2, _azi2, _s12, _m12, _M12, _M21, _S12) = |
2532 | | line._gen_position(false, 1000.0, caps::CAP_NONE); |
2533 | | assert!(a12.is_nan()); |
2534 | | } |
2535 | | |
2536 | | #[test] |
2537 | | fn test_std_geodesic_geodsolve84() { |
2538 | | // Tests for python implementation to check fix for range errors with |
2539 | | // {fmod,sin,cos}(inf) (includes GeodSolve84 - GeodSolve91). |
2540 | | let geod = Geodesic::wgs84(); |
2541 | | let (lat2, lon2, azi2) = geod.direct(0.0, 0.0, 90.0, f64::INFINITY); |
2542 | | assert!(lat2.is_nan()); |
2543 | | assert!(lon2.is_nan()); |
2544 | | assert!(azi2.is_nan()); |
2545 | | let (lat2, lon2, azi2) = geod.direct(0.0, 0.0, 90.0, f64::NAN); |
2546 | | assert!(lat2.is_nan()); |
2547 | | assert!(lon2.is_nan()); |
2548 | | assert!(azi2.is_nan()); |
2549 | | let (lat2, lon2, azi2) = geod.direct(0.0, 0.0, f64::INFINITY, 1000.0); |
2550 | | assert!(lat2.is_nan()); |
2551 | | assert!(lon2.is_nan()); |
2552 | | assert!(azi2.is_nan()); |
2553 | | let (lat2, lon2, azi2) = geod.direct(0.0, 0.0, f64::NAN, 1000.0); |
2554 | | assert!(lat2.is_nan()); |
2555 | | assert!(lon2.is_nan()); |
2556 | | assert!(azi2.is_nan()); |
2557 | | let (lat2, lon2, azi2) = geod.direct(0.0, f64::INFINITY, 90.0, 1000.0); |
2558 | | assert_eq!(lat2, 0.0); |
2559 | | assert!(lon2.is_nan()); |
2560 | | assert_eq!(azi2, 90.0); |
2561 | | let (lat2, lon2, azi2) = geod.direct(0.0, f64::NAN, 90.0, 1000.0); |
2562 | | assert_eq!(lat2, 0.0); |
2563 | | assert!(lon2.is_nan()); |
2564 | | assert_eq!(azi2, 90.0); |
2565 | | let (lat2, lon2, azi2) = geod.direct(f64::INFINITY, 0.0, 90.0, 1000.0); |
2566 | | assert!(lat2.is_nan()); |
2567 | | assert!(lon2.is_nan()); |
2568 | | assert!(azi2.is_nan()); |
2569 | | let (lat2, lon2, azi2) = geod.direct(f64::NAN, 0.0, 90.0, 1000.0); |
2570 | | assert!(lat2.is_nan()); |
2571 | | assert!(lon2.is_nan()); |
2572 | | assert!(azi2.is_nan()); |
2573 | | } |
2574 | | |
2575 | | // *_geodtest_* tests are based on Karney's GeodTest*.dat test datasets. |
2576 | | // A description of these files' content can be found at: |
2577 | | // https://geographiclib.sourceforge.io/html/geodesic.html#testgeod |
2578 | | // Here are some key excerpts... |
2579 | | // This consists of a set of geodesics for the WGS84 ellipsoid. |
2580 | | // Each line of the test set gives 10 space delimited numbers |
2581 | | // latitude at point 1, lat1 (degrees, exact) |
2582 | | // longitude at point 1, lon1 (degrees, always 0) |
2583 | | // azimuth at point 1, azi1 (clockwise from north in degrees, exact) |
2584 | | // latitude at point 2, lat2 (degrees, accurate to 10−18 deg) |
2585 | | // longitude at point 2, lon2 (degrees, accurate to 10−18 deg) |
2586 | | // azimuth at point 2, azi2 (degrees, accurate to 10−18 deg) |
2587 | | // geodesic distance from point 1 to point 2, s12 (meters, exact) |
2588 | | // arc distance on the auxiliary sphere, a12 (degrees, accurate to 10−18 deg) |
2589 | | // reduced length of the geodesic, m12 (meters, accurate to 0.1 pm) |
2590 | | // the area under the geodesic, S12 (m2, accurate to 1 mm2) |
2591 | | |
2592 | | static FULL_TEST_PATH: &str = "test_fixtures/test_data_unzipped/GeodTest.dat"; |
2593 | | static SHORT_TEST_PATH: &str = "test_fixtures/test_data_unzipped/GeodTest-short.dat"; |
2594 | | static BUILTIN_TEST_PATH: &str = "test_fixtures/GeodTest-100.dat"; |
2595 | | fn test_input_path() -> &'static str { |
2596 | | if cfg!(feature = "test_full") { |
2597 | | FULL_TEST_PATH |
2598 | | } else if cfg!(feature = "test_short") { |
2599 | | SHORT_TEST_PATH |
2600 | | } else { |
2601 | | BUILTIN_TEST_PATH |
2602 | | } |
2603 | | } |
2604 | | |
2605 | | fn geodtest_basic<T>(path: &str, f: T) |
2606 | | where |
2607 | | T: Fn(usize, &(f64, f64, f64, f64, f64, f64, f64, f64, f64, f64)), |
2608 | | { |
2609 | | let dir_base = std::env::current_dir().expect("Failed to determine current directory"); |
2610 | | let path_base = dir_base.as_path(); |
2611 | | let pathbuf = std::path::Path::new(path_base).join(path); |
2612 | | let path = pathbuf.as_path(); |
2613 | | let file = match std::fs::File::open(path) { |
2614 | | Ok(val) => val, |
2615 | | Err(_error) => { |
2616 | | let path_str = path |
2617 | | .to_str() |
2618 | | .expect("Failed to convert GeodTest path to string during error reporting"); |
2619 | | panic!("Failed to open test input file. Run `script/download-test-data.sh` to download test input to: {}\nFor details see https://geographiclib.sourceforge.io/html/geodesic.html#testgeod", path_str) |
2620 | | } |
2621 | | }; |
2622 | | let reader = std::io::BufReader::new(file); |
2623 | | reader.lines().enumerate().for_each(|(i, line)| { |
2624 | | let line_safe = line.expect("Failed to read line"); |
2625 | | let items: Vec<f64> = line_safe |
2626 | | .split(' ') |
2627 | | .enumerate() |
2628 | | .map(|(j, item)| match item.parse::<f64>() { |
2629 | | Ok(parsed) => parsed, |
2630 | | Err(_error) => { |
2631 | | panic!("Error parsing item {} on line {}: {}", j + 1, i + 1, item) |
2632 | | } |
2633 | | }) |
2634 | | .collect(); |
2635 | | assert_eq!(items.len(), 10); |
2636 | | let tuple = ( |
2637 | | items[0], items[1], items[2], items[3], items[4], items[5], items[6], items[7], |
2638 | | items[8], items[9], |
2639 | | ); |
2640 | | f(i + 1, &tuple); // report 1-based line number rather than 0-based |
2641 | | }); |
2642 | | } |
2643 | | |
2644 | | #[test] |
2645 | | fn test_geodtest_geodesic_direct12() { |
2646 | | let g = std::sync::Arc::new(std::sync::Mutex::new(Geodesic::wgs84())); |
2647 | | |
2648 | | geodtest_basic( |
2649 | | test_input_path(), |
2650 | | |_line_num, &(lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, S12)| { |
2651 | | let g = g.lock().unwrap(); |
2652 | | let (lat2_out, lon2_out, azi2_out, m12_out, _M12_out, _M21_out, S12_out, a12_out) = |
2653 | | g.direct(lat1, lon1, azi1, s12); |
2654 | | assert_relative_eq!(lat2, lat2_out, epsilon = 1e-13); |
2655 | | assert_relative_eq!(lon2, lon2_out, epsilon = 2e-8); |
2656 | | assert_relative_eq!(azi2, azi2_out, epsilon = 2e-8); |
2657 | | assert_relative_eq!(m12, m12_out, epsilon = 9e-9); |
2658 | | assert_relative_eq!(S12, S12_out, epsilon = 2e4); // Note: unreasonable tolerance |
2659 | | assert_relative_eq!(a12, a12_out, epsilon = 9e-14); |
2660 | | }, |
2661 | | ); |
2662 | | } |
2663 | | |
2664 | | #[test] |
2665 | | fn test_geodtest_geodesic_direct21() { |
2666 | | let g = std::sync::Arc::new(std::sync::Mutex::new(Geodesic::wgs84())); |
2667 | | |
2668 | | geodtest_basic( |
2669 | | test_input_path(), |
2670 | | |_line_num, &(lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, S12)| { |
2671 | | let g = g.lock().unwrap(); |
2672 | | // Reverse some values for 2->1 instead of 1->2 |
2673 | | let (lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, S12) = |
2674 | | (lat2, lon2, azi2, lat1, lon1, azi1, -s12, -a12, -m12, -S12); |
2675 | | let (lat2_out, lon2_out, azi2_out, m12_out, _M12_out, _M21_out, S12_out, a12_out) = |
2676 | | g.direct(lat1, lon1, azi1, s12); |
2677 | | assert_relative_eq!(lat2, lat2_out, epsilon = 8e-14); |
2678 | | assert_relative_eq!(lon2, lon2_out, epsilon = 4e-6); |
2679 | | assert_relative_eq!(azi2, azi2_out, epsilon = 4e-6); |
2680 | | assert_relative_eq!(m12, m12_out, epsilon = 1e-8); |
2681 | | assert_relative_eq!(S12, S12_out, epsilon = 3e6); // Note: unreasonable tolerance |
2682 | | assert_relative_eq!(a12, a12_out, epsilon = 9e-14); |
2683 | | }, |
2684 | | ); |
2685 | | } |
2686 | | |
2687 | | #[test] |
2688 | | fn test_geodtest_geodesic_inverse12() { |
2689 | | let g = std::sync::Arc::new(std::sync::Mutex::new(Geodesic::wgs84())); |
2690 | | |
2691 | | geodtest_basic( |
2692 | | test_input_path(), |
2693 | | |line_num, &(lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, S12)| { |
2694 | | let g = g.lock().unwrap(); |
2695 | | let (s12_out, azi1_out, azi2_out, m12_out, _M12_out, _M21_out, S12_out, a12_out) = |
2696 | | g.inverse(lat1, lon1, lat2, lon2); |
2697 | | assert_relative_eq!(s12, s12_out, epsilon = 8e-9); |
2698 | | assert_relative_eq!(azi1, azi1_out, epsilon = 2e-2); |
2699 | | assert_relative_eq!(azi2, azi2_out, epsilon = 2e-2); |
2700 | | assert_relative_eq!(m12, m12_out, epsilon = 5e-5); |
2701 | | // Our area calculation differs significantly (~1e7) from the value in GeodTest.dat for |
2702 | | // line 400001, BUT our value also perfectly matches the value returned by GeographicLib |
2703 | | // (C++) 1.51. Here's the problem line, for reference: |
2704 | | // 4.199535552987 0 90 -4.199535552987 179.398106343454992238 90 19970505.608097404994 180 0 |
2705 | | if line_num != 400001 { |
2706 | | assert_relative_eq!(S12, S12_out, epsilon = 3e10); // Note: unreasonable tolerance |
2707 | | } |
2708 | | assert_relative_eq!(a12, a12_out, epsilon = 2e-10); |
2709 | | }, |
2710 | | ); |
2711 | | } |
2712 | | |
2713 | | #[test] |
2714 | | fn test_turnaround() { |
2715 | | let g = Geodesic::wgs84(); |
2716 | | |
2717 | | let start = (0.0, 0.0); |
2718 | | let destination = (0.0, 1.0); |
2719 | | |
2720 | | let (distance, azi1, _, _) = g.inverse(start.0, start.1, destination.0, destination.1); |
2721 | | |
2722 | | // Confirm that we've gone due-east |
2723 | | assert_eq!(azi1, 90.0); |
2724 | | |
2725 | | // Turn around by adding 180 degrees to the azimuth |
2726 | | let turn_around = azi1 + 180.0; |
2727 | | |
2728 | | // Confirm that turn around is due west |
2729 | | assert_eq!(turn_around, 270.0); |
2730 | | |
2731 | | // Test that we can turn around and get back to the starting point. |
2732 | | let (lat, lon) = g.direct(destination.0, destination.1, turn_around, distance); |
2733 | | assert_relative_eq!(lat, start.0, epsilon = 1.0e-3); |
2734 | | assert_relative_eq!(lon, start.1, epsilon = 1.0e-3); |
2735 | | } |
2736 | | } |