GeographicLib: GeodesicExact.cpp Source File (original) (raw)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
31#include
32
33#if defined(_MSC_VER)
34
35# pragma warning (disable: 4701)
36#endif
37
39
40 using namespace std;
41
42 GeodesicExact::GeodesicExact(real a, real f)
43 : maxit2_(maxit1_ + Math::digits() + 10)
44
45
46
47 , tiny_(sqrt(numeric_limits::min()))
48 , tol0_(numeric_limits::epsilon())
49
50
51
52 , tol1_(200 * tol0_)
53 , tol2_(sqrt(tol0_))
54 , tolb_(tol0_)
55 , xthresh_(1000 * tol2_)
56 , _a(a)
57 , _f(f)
58 , _f1(1 - _f)
59 , _e2(_f * (2 - _f))
60 , _ep2(_e2 / Math::sq(_f1))
61 , _n(_f / ( 2 - _f))
62 , _b(_a * _f1)
63
64
65
66 , _c2((Math::sq(_a) + Math::sq(_b) *
67 (_f == 0 ? 1 :
68 (_f > 0 ? asinh(sqrt(_ep2)) : atan(sqrt(-_e2))) /
69 sqrt(fabs(_e2))))/2)
70
71
72
73
74
75
76
77
78
79
80 , _etol2(real(0.1) * tol2_ /
81 sqrt( fmax(real(0.001), fabs(_f)) * fmin(real(1), 1 - _f/2) / 2 ))
82 {
83 if (!(isfinite(_a) && _a > 0))
84 throw GeographicErr("Equatorial radius is not positive");
85 if (!(isfinite(_b) && _b > 0))
86 throw GeographicErr("Polar semi-axis is not positive");
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299 static const int ndiv = 100;
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319#if GEOGRAPHICLIB_PRECISION == 1
320 static const unsigned char narr[2*ndiv+1] = {
321 19,18,16,15,14,13,13,13,12,12,11,11,11,11,10,10,10,10,9,9,9,9,9,9,9,9,8,
322 8,8,8,8,8,8,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,
323 5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,
324 2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,
325 4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,
326 6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8
327 };
328#elif GEOGRAPHICLIB_PRECISION == 2
329 static const unsigned char narr[2*ndiv+1] = {
330 22,21,19,18,17,17,16,15,15,15,14,14,14,13,13,13,13,13,13,12,12,12,12,12,
331 12,11,11,11,11,11,11,11,11,11,11,10,10,10,10,10,10,10,10,9,9,9,9,9,9,9,9,
332 9,9,9,9,9,9,8,8,8,8,8,8,8,8,8,8,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,
333 6,6,5,5,5,5,5,5,5,5,4,4,3,2,3,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,7,7,
334 7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,9,9,
335 9,9,9,9,9,10,10,10,10,10,10,10,10,10,11,11,11,11,11,11,11,11,11,11,12,12,
336 12,12,12,13,13,13,13,13,14,14,15,15,16,17,18,19
337 };
338#elif GEOGRAPHICLIB_PRECISION == 3
339 static const unsigned char narr[2*ndiv+1] = {
340 23,22,20,19,18,17,17,16,16,15,15,15,15,14,14,14,14,13,13,13,13,13,13,13,
341 12,12,12,12,12,12,11,11,11,11,11,11,11,11,11,11,11,10,10,10,10,10,10,10,
342 10,10,10,9,9,9,9,9,9,9,9,9,9,9,9,9,9,8,8,8,8,8,8,8,8,8,8,8,7,7,7,7,7,7,7,
343 7,7,7,7,7,6,6,6,6,6,6,5,5,5,5,5,4,2,4,5,5,5,5,5,6,6,6,6,6,6,6,7,7,7,7,7,
344 7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,
345 10,10,10,10,10,10,10,10,10,10,11,11,11,11,11,11,11,11,11,11,11,12,12,12,
346 12,12,12,13,13,13,13,13,13,13,14,14,14,15,15,15,16,16,17,18,19,20
347 };
348#elif GEOGRAPHICLIB_PRECISION == 4
349 static const unsigned char narr[2*ndiv+1] = {
350 25,24,22,21,20,19,19,18,18,17,17,17,17,16,16,16,15,15,15,15,15,15,15,14,
351 14,14,14,14,14,13,13,13,13,13,13,13,13,13,13,13,13,12,12,12,12,12,12,12,
352 12,12,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,10,10,10,10,10,10,10,
353 10,10,10,9,9,9,9,9,9,9,9,9,9,9,9,9,8,8,8,8,8,8,7,7,7,7,7,6,2,6,7,7,7,7,7,
354 8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,
355 11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,12,
356 12,13,13,13,13,13,13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,
357 15,16,16,16,17,17,17,17,18,18,19,20,21,23,24
358 };
359#elif GEOGRAPHICLIB_PRECISION == 5
360 static const unsigned char narr[2*ndiv+1] = {
361 27,26,24,23,22,22,21,21,20,20,20,19,19,19,19,18,18,18,18,18,17,17,17,17,
362 17,17,17,17,16,16,16,16,16,16,16,15,15,15,15,15,15,15,15,15,15,15,15,14,
363 14,14,14,14,14,14,14,14,14,14,13,13,13,13,13,13,13,13,13,13,13,13,13,13,
364 12,12,12,12,12,12,12,12,12,12,11,11,11,11,11,11,11,11,11,11,11,10,10,10,
365 10,9,9,9,2,9,9,9,10,10,10,10,10,11,11,11,11,11,11,11,11,11,11,12,12,12,
366 12,12,12,12,12,12,12,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,14,14,
367 14,14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15,15,15,15,16,16,16,
368 16,16,16,16,17,17,17,17,17,17,17,17,18,18,18,18,18,19,19,19,19,20,20,21,
369 21,21,22,23,24,26,27
370 };
371#else
372#error "Bad value for GEOGRAPHICLIB_PRECISION"
373#endif
374 real n = ndiv * _n;
375 int j = ndiv + int(n < 0 ? floor(n) : ceil(n));
376 int N = int(narr[j]);
377
378 N = (N % 2 == 0 ? 2 : 3) * (1 << (N/2));
379#if GEOGRAPHICLIB_PRECISION == 5
381
382
384 while (N < M) N = N % 3 == 0 ? 4*N/3 : 3*N/2;
385 }
386#endif
388 _nC4 = N;
389 }
390
396
398 unsigned caps) const {
400 }
401
403 bool arcmode, real s12_a12,
404 unsigned outmask,
405 real& lat2, real& lon2, real& azi2,
406 real& s12, real& m12,
407 real& M12, real& M21,
408 real& S12) const {
409
412 .
413 GenPosition(arcmode, s12_a12, outmask,
414 lat2, lon2, azi2, s12, m12, M12, M21, S12);
415 }
416
418 real azi1,
419 bool arcmode, real s12_a12,
420 unsigned caps) const {
422 real salp1, calp1;
423
425
428 caps, arcmode, s12_a12);
429 }
430
432 real azi1, real s12,
433 unsigned caps) const {
434 return GenDirectLine(lat1, lon1, azi1, false, s12, caps);
435 }
436
438 real azi1, real a12,
439 unsigned caps) const {
440 return GenDirectLine(lat1, lon1, azi1, true, a12, caps);
441 }
442
445 unsigned outmask, real& s12,
449 real& S12) const {
450
451
452
453 using std::isnan;
455
456 int lonsign = signbit(lon12) ? -1 : 1;
457 lon12 *= lonsign; lon12s *= lonsign;
460 slam12, clam12;
461
463
464 lon12s = (Math::hd - lon12) - lon12s;
465
466
469
470
471 int swapp = fabs(lat1) < fabs(lat2) || isnan(lat2) ? -1 : 1;
472 if (swapp < 0) {
473 lonsign *= -1;
474 swap(lat1, lat2);
475 }
476
477 int latsign = signbit(lat1) ? 1 : -1;
478 lat1 *= latsign;
479 lat2 *= latsign;
480
481
482
483
484
485
486
487
488
489
490
491
492 real sbet1, cbet1, sbet2, cbet2, s12x, m12x = Math::NaN();
493
494
495 EllipticFunction E(-_ep2);
496
498
499
500 Math::norm(sbet1, cbet1); cbet1 = fmax(tiny_, cbet1);
501
503
504 Math::norm(sbet2, cbet2); cbet2 = fmax(tiny_, cbet2);
505
506
507
508
509
510
511
512
513
514 if (cbet1 < -sbet1) {
515 if (cbet2 == cbet1)
516 sbet2 = copysign(sbet1, sbet2);
517 } else {
518 if (fabs(sbet2) == -sbet1)
519 cbet2 = cbet1;
520 }
521
523 dn1 = (_f >= 0 ? sqrt(1 + _ep2 * Math::sq(sbet1)) :
524 sqrt(1 - _e2 * Math::sq(cbet1)) / _f1),
525 dn2 = (_f >= 0 ? sqrt(1 + _ep2 * Math::sq(sbet2)) :
526 sqrt(1 - _e2 * Math::sq(cbet2)) / _f1);
527
528 real a12, sig12;
529
530 bool meridian = lat1 == -Math::qd || slam12 == 0;
531
532 if (meridian) {
533
534
535
536
537 calp1 = clam12; salp1 = slam12;
538 calp2 = 1; salp2 = 0;
539
541
542 ssig1 = sbet1, csig1 = calp1 * cbet1,
543 ssig2 = sbet2, csig2 = calp2 * cbet2;
544
545
546 sig12 = atan2(fmax(real(0), csig1 * ssig2 - ssig1 * csig2),
547 csig1 * csig2 + ssig1 * ssig2);
548 {
550 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
552 s12x, m12x, dummy, M12, M21);
553 }
554
555
556
557
558
559
560
561 if (sig12 < 1 || m12x >= 0) {
562
563 if (sig12 < 3 * tiny_ ||
564
565 (sig12 < tol0_ && (s12x < 0 || m12x < 0)))
566 sig12 = m12x = s12x = 0;
567 m12x *= _b;
568 s12x *= _b;
570 } else
571
572 meridian = false;
573 }
574
575
576 real omg12 = 0, somg12 = 2, comg12 = 0;
577 if (!meridian &&
578 sbet1 == 0 &&
579 (_f <= 0 || lon12s >= _f * Math::hd)) {
580
581
582 calp1 = calp2 = 0; salp1 = salp2 = 1;
583 s12x = _a * lam12;
584 sig12 = omg12 = lam12 / _f1;
585 m12x = _b * sin(sig12);
587 M12 = M21 = cos(sig12);
588 a12 = lon12 / _f1;
589
590 } else if (!meridian) {
591
592
593
594
595
597 sig12 = InverseStart(E, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
598 lam12, slam12, clam12,
599 salp1, calp1, salp2, calp2, dnm);
600
601 if (sig12 >= 0) {
602
603 s12x = sig12 * _b * dnm;
604 m12x = Math::sq(dnm) * _b * sin(sig12 / dnm);
606 M12 = M21 = cos(sig12 / dnm);
608 omg12 = lam12 / (_f1 * dnm);
609 } else {
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, domg12 = 0;
625 unsigned numit = 0;
626
627 real salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
628 for (bool tripn = false, tripb = false;; ++numit) {
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
651 real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
652 slam12, clam12,
653 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
654 E, domg12, numit < maxit1_, dv);
655 if (tripb ||
656
657 !(fabs(v) >= (tripn ? 8 : 1) * tol0_) ||
658
659 numit == maxit2_)
660 break;
661
662 if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
663 { salp1b = salp1; calp1b = calp1; }
664 else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
665 { salp1a = salp1; calp1a = calp1; }
666 if (numit < maxit1_ && dv > 0) {
668 dalp1 = -v/dv;
669
670
671
672 if (fabs(dalp1) < Math::pi()) {
674 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
675 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
676 if (nsalp1 > 0) {
677 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
678 salp1 = nsalp1;
680
681
682
683 tripn = fabs(v) <= 16 * tol0_;
684 continue;
685 }
686 }
687 }
688
689
690
691
692
693
694
695
696 salp1 = (salp1a + salp1b)/2;
697 calp1 = (calp1a + calp1b)/2;
699 tripn = false;
700 tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
701 fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
702 }
703 {
705 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
706 cbet1, cbet2, outmask, s12x, m12x, dummy, M12, M21);
707 }
708 m12x *= _b;
709 s12x *= _b;
711 if (outmask & AREA) {
712
713 real sdomg12 = sin(domg12), cdomg12 = cos(domg12);
714 somg12 = slam12 * cdomg12 - clam12 * sdomg12;
715 comg12 = clam12 * cdomg12 + slam12 * sdomg12;
716 }
717 }
718 }
719
721 s12 = real(0) + s12x;
722
724 m12 = real(0) + m12x;
725
726 if (outmask & AREA) {
728
729 salp0 = salp1 * cbet1,
730 calp0 = hypot(calp1, salp1 * sbet1);
732
733 A4 = Math::sq(_a) * calp0 * salp0 * _e2;
734 if (A4 != 0) {
737
738 ssig1 = sbet1, csig1 = calp1 * cbet1,
739 ssig2 = sbet2, csig2 = calp2 * cbet2;
742 I4Integrand i4(_ep2, k2);
743 vector C4a(_nC4);
745 S12 = A4 * DST::integral(ssig1, csig1, ssig2, csig2, C4a.data(), _nC4);
746 } else
747
748 S12 = 0;
749
750 if (!meridian && somg12 == 2) {
751 somg12 = sin(omg12); comg12 = cos(omg12);
752 }
753
754 if (!meridian &&
755
756 comg12 > -real(0.7071) &&
757 sbet2 - sbet1 < real(1.75)) {
758
759
760
761 real domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
762 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
763 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
764 } else {
765
767 salp12 = salp2 * calp1 - calp2 * salp1,
768 calp12 = calp2 * calp1 + salp2 * salp1;
769
770
771
772
773 if (salp12 == 0 && calp12 < 0) {
774 salp12 = tiny_ * calp1;
775 calp12 = -1;
776 }
777 alp12 = atan2(salp12, calp12);
778 }
779 S12 += _c2 * alp12;
780 S12 *= swapp * lonsign * latsign;
781
782 S12 += 0;
783 }
784
785
786 if (swapp < 0) {
787 swap(salp1, salp2);
788 swap(calp1, calp2);
790 swap(M12, M21);
791 }
792
793 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
794 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
795
796
797 return a12;
798 }
799
800 Math::real GeodesicExact::GenInverse(real lat1, real lon1,
801 real lat2, real lon2,
802 unsigned outmask,
803 real& s12, real& azi1, real& azi2,
804 real& m12, real& M12, real& M21,
805 real& S12) const {
806 outmask &= OUT_MASK;
807 real salp1, calp1, salp2, calp2,
808 a12 = GenInverse(lat1, lon1, lat2, lon2,
809 outmask, s12, salp1, calp1, salp2, calp2,
810 m12, M12, M21, S12);
814 }
815 return a12;
816 }
817
819 real lat2, real lon2,
820 unsigned caps) const {
821 real t, salp1, calp1, salp2, calp2,
822 a12 = GenInverse(lat1, lon1, lat2, lon2,
823
824 0u, t, salp1, calp1, salp2, calp2,
825 t, t, t, t),
827
829 return GeodesicLineExact(*this, lat1, lon1, azi1, salp1, calp1, caps,
830 true, a12);
831 }
832
837 real cbet1, real cbet2, unsigned outmask,
839 real& M12, real& M21) const {
840
841
842
843 outmask &= OUT_ALL;
844
845
846
847
848
849
850
852
853 s12b = E.E() / (Math::pi() / 2) *
854 (sig12 + (E.deltaE(ssig2, csig2, dn2) - E.deltaE(ssig1, csig1, dn1)));
857 m0x = - E.k2() * E.D() / (Math::pi() / 2),
858 J12 = m0x *
859 (sig12 + (E.deltaD(ssig2, csig2, dn2) - E.deltaD(ssig1, csig1, dn1)));
861 m0 = m0x;
862
863
864
865 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
866 csig1 * csig2 * J12;
867 }
869 real csig12 = csig1 * csig2 + ssig1 * ssig2;
870 real t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
871 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
872 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
873 }
874 }
875 }
876
878
879
884 r = (p + q - 1) / 6;
885 if ( !(q == 0 && r <= 0) ) {
887
888
889 S = p * q / 4,
891 r3 = r * r2,
892
893
894 disc = S * (S + 2 * r3);
896 if (disc >= 0) {
897 real T3 = S + r3;
898
899
900
901 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
902
903 real T = cbrt(T3);
904
905 u += T + (T != 0 ? r2 / T : 0);
906 } else {
907
908 real ang = atan2(sqrt(-disc), -(S + r3));
909
910
911 u += 2 * r * cos(ang / 3);
912 }
914 v = sqrt(Math::sq(u) + q),
915
916 uv = u < 0 ? q / (v - u) : u + v,
917 w = (uv - q) / (2 * v);
918
919
920 k = uv / (sqrt(uv + Math::sq(w)) + w);
921 } else {
922
923
924 k = 0;
925 }
926 return k;
927 }
928
929 Math::real GeodesicExact::InverseStart(EllipticFunction& E,
934
936
937 real& dnm) const {
938
939
940
942 sig12 = -1,
943
944 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
945 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
946 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
947 bool shortline = cbet12 >= 0 && sbet12 < real(0.5) &&
948 cbet2 * lam12 < real(0.5);
949 real somg12, comg12;
950 if (shortline) {
952
953
954 sbetm2 /= sbetm2 + Math::sq(cbet1 + cbet2);
955 dnm = sqrt(1 + _ep2 * sbetm2);
956 real omg12 = lam12 / (_f1 * dnm);
957 somg12 = sin(omg12); comg12 = cos(omg12);
958 } else {
959 somg12 = slam12; comg12 = clam12;
960 }
961
962 salp1 = cbet2 * somg12;
963 calp1 = comg12 >= 0 ?
964 sbet12 + cbet2 * sbet1 * Math::sq(somg12) / (1 + comg12) :
965 sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
966
968 ssig12 = hypot(salp1, calp1),
969 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
970
971 if (shortline && ssig12 < _etol2) {
972
973 salp2 = cbet1 * somg12;
974 calp2 = sbet12 - cbet1 * sbet2 *
975 (comg12 >= 0 ? Math::sq(somg12) / (1 + comg12) : 1 - comg12);
977
978 sig12 = atan2(ssig12, csig12);
979 } else if (fabs(_n) > real(0.1) ||
980 csig12 >= 0 ||
982
983 } else {
984
985
986 real x, y, lamscale, betscale;
987 real lam12x = atan2(-slam12, -clam12);
988 if (_f >= 0) {
989
990 {
992 E.Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
993 lamscale = _e2/_f1 * cbet1 * 2 * E.H();
994 }
995 betscale = lamscale * cbet1;
996
997 x = lam12x / lamscale;
998 y = sbet12a / betscale;
999 } else {
1000
1002 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1003 bet12a = atan2(sbet12a, cbet12a);
1004 real m12b, m0, dummy;
1005
1006
1007 Lengths(E, Math::pi() + bet12a,
1008 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1009 cbet1, cbet2, REDUCEDLENGTH, dummy, m12b, m0, dummy, dummy);
1010 x = -1 + m12b / (cbet1 * cbet2 * m0 * Math::pi());
1011 betscale = x < -real(0.01) ? sbet12a / x :
1013 lamscale = betscale / cbet1;
1014 y = lam12x / lamscale;
1015 }
1016
1017 if (y > -tol1_ && x > -1 - xthresh_) {
1018
1019
1020 if (_f >= 0) {
1021 salp1 = fmin(real(1), -x); calp1 = - sqrt(1 - Math::sq(salp1));
1022 } else {
1023 calp1 = fmax(real(x > -tol1_ ? 0 : -1), x);
1024 salp1 = sqrt(1 - Math::sq(calp1));
1025 }
1026 } else {
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061 real k = Astroid(x, y);
1063 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
1064 somg12 = sin(omg12a); comg12 = -cos(omg12a);
1065
1066 salp1 = cbet2 * somg12;
1067 calp1 = sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
1068 }
1069 }
1070
1071 if (!(salp1 <= 0))
1073 else {
1074 salp1 = 1; calp1 = 0;
1075 }
1076 return sig12;
1077 }
1078
1084 real& sig12,
1087 EllipticFunction& E,
1088 real& domg12,
1089 bool diffp, real& dlam12) const
1090 {
1091
1092 if (sbet1 == 0 && calp1 == 0)
1093
1094
1095 calp1 = -tiny_;
1096
1098
1099 salp0 = salp1 * cbet1,
1100 calp0 = hypot(calp1, salp1 * sbet1);
1101
1102 real somg1, comg1, somg2, comg2, somg12, comg12, cchi1, cchi2, lam12;
1103
1104
1105 ssig1 = sbet1; somg1 = salp0 * sbet1;
1106 csig1 = comg1 = calp1 * cbet1;
1107
1108 cchi1 = _f1 * dn1 * comg1;
1110
1111
1112
1113
1114
1115
1116
1117 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
1118
1119
1120
1121
1122 calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ?
1123 sqrt(Math::sq(calp1 * cbet1) +
1124 (cbet1 < -sbet1 ?
1125 (cbet2 - cbet1) * (cbet1 + cbet2) :
1126 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1127 fabs(calp1);
1128
1129
1130 ssig2 = sbet2; somg2 = salp0 * sbet2;
1131 csig2 = comg2 = calp2 * cbet2;
1132
1133 cchi2 = _f1 * dn2 * comg2;
1135
1136
1137
1138
1139 sig12 = atan2(fmax(real(0), csig1 * ssig2 - ssig1 * csig2),
1140 csig1 * csig2 + ssig1 * ssig2);
1141
1142
1143 somg12 = fmax(real(0), comg1 * somg2 - somg1 * comg2);
1144 comg12 = comg1 * comg2 + somg1 * somg2;
1146 E.Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
1147
1149 schi12 = fmax(real(0), cchi1 * somg2 - somg1 * cchi2),
1150 cchi12 = cchi1 * cchi2 + somg1 * somg2;
1151
1152 real eta = atan2(schi12 * clam120 - cchi12 * slam120,
1153 cchi12 * clam120 + schi12 * slam120);
1154 real deta12 = -_e2/_f1 * salp0 * E.H() / (Math::pi() / 2) *
1155 (sig12 + (E.deltaH(ssig2, csig2, dn2) - E.deltaH(ssig1, csig1, dn1)));
1156 lam12 = eta + deta12;
1157
1158 domg12 = deta12 + atan2(schi12 * comg12 - cchi12 * somg12,
1159 cchi12 * comg12 + schi12 * somg12);
1160 if (diffp) {
1161 if (calp2 == 0)
1162 dlam12 = - 2 * _f1 * dn1 / sbet1;
1163 else {
1165 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1167 dummy, dlam12, dummy, dummy, dummy);
1168 dlam12 *= _f1 / (calp2 * cbet2);
1169 }
1170 }
1171
1172 return lam12;
1173 }
1174
1175 Math::real GeodesicExact::I4Integrand::asinhsqrt(real x) {
1176
1177 return x == 0 ? 1 :
1178 (x > 0 ? asinh(sqrt(x))/sqrt(x) :
1179 asin(sqrt(-x))/sqrt(-x));
1180 }
1182
1183
1184
1185
1186
1187
1188 return x + (sqrt(1 + x) * asinhsqrt(x) - 1);
1189 }
1191
1192 return x == 0 ? 4/real(3) :
1193
1194 1 + (1 - asinhsqrt(x) / sqrt(1+x)) / (2*x);
1195 }
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211 Math::real GeodesicExact::I4Integrand::DtX(real y) const {
1212
1213
1214 if (X == y) return tdX;
1215 if (X * y <= 0) return ( tX - t(y) ) / (X - y);
1217 sy = sqrt(fabs(y)), sy1 = sqrt(1 + y),
1218 z = (X - y) / (sX * sy1 + sy * sX1),
1219 d1 = 2 * sX * sy,
1220 d2 = 2 * (X * sy * sy1 + y * sXX1);
1221 return X > 0 ?
1222 ( 1 + (asinh(z)/z) / d1 - (asinhsX + asinh(sy)) / d2 ) :
1223
1224 ( 1 - (asin (z)/z) / d1 - (asinhsX + asin (sy)) / d2 );
1225 }
1226 GeodesicExact::I4Integrand::I4Integrand(real ep2, real k2)
1227 : X( ep2 )
1228 , tX( t(X) )
1229 , tdX( td(X) )
1230 , _k2( k2 )
1231 {
1232 sX = sqrt(fabs(X));
1233 sX1 = sqrt(1 + X);
1234 sXX1 = sX * sX1;
1235 asinhsX = X > 0 ? asinh(sX) : asin(sX);
1236 }
1237 Math::real GeodesicExact::I4Integrand::operator()(real sig) const {
1238 real ssig = sin(sig);
1239 return - DtX(_k2 * Math::sq(ssig)) * ssig/2;
1240 }
1241
1242}
GeographicLib::Math::real real
Header for GeographicLib::GeodesicExact class.
Header for GeographicLib::GeodesicLineExact class.
void transform(std::function< real(real)> f, real F[]) const
static real integral(real sinx, real cosx, const real F[], int N)
Elliptic integrals and functions.
Math::real deltaE(real sn, real cn, real dn) const
Math::real deltaD(real sn, real cn, real dn) const
Exact geodesic calculations.
GeodesicLineExact InverseLine(real lat1, real lon1, real lat2, real lon2, unsigned caps=ALL) const
Definition GeodesicExact.cpp:818
GeodesicLineExact GenDirectLine(real lat1, real lon1, real azi1, bool arcmode, real s12_a12, unsigned caps=ALL) const
Definition GeodesicExact.cpp:417
GeodesicLineExact DirectLine(real lat1, real lon1, real azi1, real s12, unsigned caps=ALL) const
Definition GeodesicExact.cpp:431
Math::real GenDirect(real lat1, real lon1, real azi1, bool arcmode, real s12_a12, unsigned outmask, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
Definition GeodesicExact.cpp:402
friend class GeodesicLineExact
GeodesicLineExact Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
Definition GeodesicExact.cpp:397
static const GeodesicExact & WGS84()
Definition GeodesicExact.cpp:391
GeodesicLineExact ArcDirectLine(real lat1, real lon1, real azi1, real a12, unsigned caps=ALL) const
Definition GeodesicExact.cpp:437
Exception handling for GeographicLib.
Mathematical functions needed by GeographicLib.
static void sincosd(T x, T &sinx, T &cosx)
static T atan2d(T y, T x)
static void norm(T &x, T &y)
static constexpr int qd
degrees per quarter turn
static T AngNormalize(T x)
static void sincosde(T x, T t, T &sinx, T &cosx)
static T AngDiff(T x, T y, T &e)
static constexpr int hd
degrees per half turn
Namespace for GeographicLib.
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)