29 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
30 #define nC1 GEOGRAPHICLIB_GEODESIC_ORDER
31 #define nC1p GEOGRAPHICLIB_GEODESIC_ORDER
32 #define nC2 GEOGRAPHICLIB_GEODESIC_ORDER
33 #define nA3 GEOGRAPHICLIB_GEODESIC_ORDER
35 #define nC3 GEOGRAPHICLIB_GEODESIC_ORDER
36 #define nC3x ((nC3 * (nC3 - 1)) / 2)
37 #define nC4 GEOGRAPHICLIB_GEODESIC_ORDER
38 #define nC4x ((nC4 * (nC4 + 1)) / 2)
43 static unsigned init = 0;
44 static const int FALSE = 0;
45 static const int TRUE = 1;
46 static unsigned digits, maxit1, maxit2;
47 static real epsilon, realmin, pi, degree, NaN,
48 tiny, tol0, tol1, tol2, tolb, xthresh;
52 #if defined(__DBL_MANT_DIG__)
53 digits = __DBL_MANT_DIG__;
57 #if defined(__DBL_EPSILON__)
58 epsilon = __DBL_EPSILON__;
60 epsilon = pow(0.5, digits - 1);
62 #if defined(__DBL_MIN__)
63 realmin = __DBL_MIN__;
65 realmin = pow(0.5, 1022);
70 pi = atan2(0.0, -1.0);
73 maxit2 = maxit1 + digits + 10;
83 xthresh = 1000 * tol2;
101 static real sq(
real x) {
return x * x; }
110 return z == 0 ? x : x * log(y) / z;
115 y = log1px(2 * y/(1 - y))/2;
116 return x < 0 ? -y : y;
120 {
return sqrt(x * x + y * y); }
124 return x < 0 ? -y : y;
128 volatile real s = u + v;
129 volatile real up = s - v;
130 volatile real vpp = s - up;
141 {
return x < y ? x : y; }
144 {
return x > y ? x : y; }
147 {
real t = *x; *x = *y; *y = t; }
149 static void SinCosNorm(
real* sinx,
real* cosx) {
150 real r = hypotx(*sinx, *cosx);
156 {
return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); }
158 {
return AngNormalize(fmod(x, (
real)(360))); }
161 real t, d = sumx(-x, y, &t);
162 if ((d - (
real)(180)) + t > (
real)(0))
164 else if ((d + (
real)(180)) + t <= (
real)(0))
171 volatile real y = fabs(x);
173 y = y < z ? z - (z - y) : y;
174 return x < 0 ? -y : y;
180 static real SinCosSeries(boolx sinp,
182 const real c[],
int n);
189 boolx scalep,
real* pM12,
real* pM21,
213 boolx diffp,
real* pdlam12,
220 static void C1f(
real eps,
real c[]);
221 static void C1pf(
real eps,
real c[]);
223 static void C2f(
real eps,
real c[]);
224 static int transit(
real lon1,
real lon2);
225 static int transitdirect(
real lon1,
real lon2);
226 static void accini(
real s[]);
227 static void acccopy(
const real s[],
real t[]);
228 static void accadd(
real s[],
real y);
230 static void accneg(
real s[]);
235 g->
f = f <= 1 ? f : 1/
f;
237 g->e2 = g->
f * (2 - g->
f);
238 g->ep2 = g->e2 / sq(g->f1);
239 g->n = g->
f / ( 2 - g->
f);
241 g->c2 = (sq(g->
a) + sq(g->b) *
243 (g->e2 > 0 ? atanhx(sqrt(g->e2)) : atan(sqrt(-g->e2))) /
244 sqrt(fabs(g->e2))))/2;
254 g->etol2 = 0.1 * tol2 /
255 sqrt( maxx((
real)(0.001), fabs(g->
f)) * minx((
real)(1), 1 - g->
f/2) / 2 );
265 real alp1, cbet1, sbet1, phi, eps;
278 l->
azi1 = AngRound(AngNormalize(azi1));
280 alp1 = l->
azi1 * degree;
283 l->salp1 = l->
azi1 == -180 ? 0 : sin(alp1);
284 l->calp1 = fabs(l->
azi1) == 90 ? 0 : cos(alp1);
287 sbet1 = l->f1 * sin(phi);
288 cbet1 = fabs(lat1) == 90 ? tiny : cos(phi);
289 SinCosNorm(&sbet1, &cbet1);
290 l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));
293 l->salp0 = l->salp1 * cbet1;
296 l->calp0 = hypotx(l->calp1, l->salp1 * sbet1);
306 l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;
307 l->csig1 = l->comg1 = sbet1 != 0 || l->calp1 != 0 ? cbet1 * l->calp1 : 1;
308 SinCosNorm(&l->ssig1, &l->csig1);
311 l->k2 = sq(l->calp0) * g->ep2;
312 eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);
314 if (l->
caps & CAP_C1) {
316 l->A1m1 = A1m1f(eps);
318 l->B11 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C1a, nC1);
319 s = sin(l->B11); c = cos(l->B11);
321 l->stau1 = l->ssig1 * c + l->csig1 * s;
322 l->ctau1 = l->csig1 * c - l->ssig1 * s;
327 if (l->
caps & CAP_C1p)
330 if (l->
caps & CAP_C2) {
331 l->A2m1 = A2m1f(eps);
333 l->B21 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C2a, nC2);
336 if (l->
caps & CAP_C3) {
338 l->A3c = -l->
f * l->salp0 * A3f(g, eps);
339 l->B31 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C3a, nC3-1);
342 if (l->
caps & CAP_C4) {
345 l->A4 = sq(l->
a) * l->calp0 * l->salp0 * g->e2;
346 l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4);
351 unsigned flags,
real s12_a12,
356 real lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,
357 m12 = 0, M12 = 0, M21 = 0, S12 = 0;
359 real sig12, ssig12, csig12, B12 = 0, AB1 = 0;
360 real omg12, lam12, lon12;
361 real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
371 outmask &= l->
caps & OUT_ALL;
380 sig12 = s12_a12 * degree;
381 s12a = fabs(s12_a12);
382 s12a -= 180 * floor(s12a / 180);
383 ssig12 = s12a == 0 ? 0 : sin(sig12);
384 csig12 = s12a == 90 ? 0 : cos(sig12);
388 tau12 = s12_a12 / (l->b * (1 + l->A1m1)),
392 B12 = - SinCosSeries(TRUE,
393 l->stau1 * c + l->ctau1 * s,
394 l->ctau1 * c - l->stau1 * s,
396 sig12 = tau12 - (B12 - l->B11);
397 ssig12 = sin(sig12); csig12 = cos(sig12);
398 if (fabs(l->
f) > 0.01) {
421 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12,
422 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12,
424 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
425 serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b;
426 sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2));
427 ssig12 = sin(sig12); csig12 = cos(sig12);
433 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
434 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
435 dn2 = sqrt(1 + l->k2 * sq(ssig2));
437 if (flags & GEOD_ARCMODE || fabs(l->
f) > 0.01)
438 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
439 AB1 = (1 + l->A1m1) * (B12 - l->B11);
442 sbet2 = l->calp0 * ssig2;
444 cbet2 = hypotx(l->salp0, l->calp0 * csig2);
447 cbet2 = csig2 = tiny;
449 salp2 = l->salp0; calp2 = l->calp0 * csig2;
452 s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
456 somg2 = l->salp0 * ssig2; comg2 = csig2;
459 - (atan2(ssig2, csig2) - atan2(l->ssig1, l->csig1))
460 + (atan2(somg2, comg2) - atan2(l->somg1, l->comg1))
461 : atan2(somg2 * l->comg1 - comg2 * l->somg1,
462 comg2 * l->comg1 + somg2 * l->somg1);
463 lam12 = omg12 + l->A3c *
464 ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
466 lon12 = lam12 / degree;
470 AngNormalize(AngNormalize(l->
lon1) + AngNormalize2(lon12));
474 lat2 = atan2(sbet2, l->f1 * cbet2) / degree;
478 azi2 = 0 - atan2(-salp2, calp2) / degree;
482 B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),
483 AB2 = (1 + l->A2m1) * (B22 - l->B21),
484 J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);
488 m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
489 - l->csig1 * csig2 * J12);
491 real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / (l->dn1 + dn2);
492 M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1;
493 M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2;
499 B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
501 if (l->calp0 == 0 || l->salp0 == 0) {
503 salp12 = salp2 * l->calp1 - calp2 * l->salp1;
504 calp12 = calp2 * l->calp1 + salp2 * l->salp1;
509 if (salp12 == 0 && calp12 < 0) {
510 salp12 = tiny * l->calp1;
522 salp12 = l->calp0 * l->salp0 *
523 (csig12 <= 0 ? l->csig1 * (1 - csig12) + ssig12 * l->ssig1 :
524 ssig12 * (l->csig1 * ssig12 / (1 + csig12) + l->ssig1));
525 calp12 = sq(l->salp0) + sq(l->calp0) * l->csig1 * csig2;
527 S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
530 if (outmask & GEOD_LATITUDE)
532 if (outmask & GEOD_LONGITUDE)
534 if (outmask & GEOD_AZIMUTH)
536 if (outmask & GEOD_DISTANCE)
538 if (outmask & GEOD_REDUCEDLENGTH)
540 if (outmask & GEOD_GEODESICSCALE) {
541 if (pM12) *pM12 = M12;
542 if (pM21) *pM21 = M21;
544 if (outmask & GEOD_AREA)
547 return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree;
552 geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0);
557 unsigned flags,
real s12_a12,
563 (plat2 ? GEOD_LATITUDE : 0U) |
564 (plon2 ? GEOD_LONGITUDE : 0U) |
565 (pazi2 ? GEOD_AZIMUTH : 0U) |
566 (ps12 ? GEOD_DISTANCE : 0U) |
567 (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
568 (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
569 (pS12 ? GEOD_AREA : 0U);
576 plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
591 real s12 = 0, azi1 = 0, azi2 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
593 int latsign, lonsign, swapp;
594 real phi, sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0;
595 real dn1, dn2, lam12, slam12, clam12;
596 real a12 = 0, sig12, calp1 = 0, salp1 = 0, calp2 = 0, salp2 = 0;
598 real C1a[nC1 + 1], C2a[nC2 + 1], C3a[nC3];
603 (ps12 ? GEOD_DISTANCE : 0U) |
604 (pazi1 || pazi2 ? GEOD_AZIMUTH : 0U) |
605 (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
606 (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
607 (pS12 ? GEOD_AREA : 0U);
613 lon12 = AngDiff(AngNormalize(lon1), AngNormalize(lon2));
615 lon12 = AngRound(lon12);
617 lonsign = lon12 >= 0 ? 1 : -1;
620 lat1 = AngRound(lat1);
621 lat2 = AngRound(lat2);
623 swapp = fabs(lat1) >= fabs(lat2) ? 1 : -1;
629 latsign = lat1 < 0 ? 1 : -1;
646 sbet1 = g->f1 * sin(phi);
647 cbet1 = lat1 == -90 ? tiny : cos(phi);
648 SinCosNorm(&sbet1, &cbet1);
652 sbet2 = g->f1 * sin(phi);
653 cbet2 = fabs(lat2) == 90 ? tiny : cos(phi);
654 SinCosNorm(&sbet2, &cbet2);
664 if (cbet1 < -sbet1) {
666 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
668 if (fabs(sbet2) == -sbet1)
672 dn1 = sqrt(1 + g->ep2 * sq(sbet1));
673 dn2 = sqrt(1 + g->ep2 * sq(sbet2));
675 lam12 = lon12 * degree;
676 slam12 = lon12 == 180 ? 0 : sin(lam12);
679 meridian = lat1 == -90 || slam12 == 0;
686 real ssig1, csig1, ssig2, csig2;
687 calp1 = clam12; salp1 = slam12;
688 calp2 = 1; salp2 = 0;
691 ssig1 = sbet1; csig1 = calp1 * cbet1;
692 ssig2 = sbet2; csig2 = calp2 * cbet2;
695 sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (
real)(0)),
696 csig1 * csig2 + ssig1 * ssig2);
699 Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
700 cbet1, cbet2, &s12x, &m12x, &dummy,
701 (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
710 if (sig12 < 1 || m12x >= 0) {
713 a12 = sig12 / degree;
722 (g->
f <= 0 || lam12 <= pi - g->f * pi)) {
725 calp1 = calp2 = 0; salp1 = salp2 = 1;
727 sig12 = omg12 = lam12 / g->f1;
728 m12x = g->b * sin(sig12);
729 if (outmask & GEOD_GEODESICSCALE)
730 M12 = M21 = cos(sig12);
733 }
else if (!meridian) {
740 sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
742 &salp1, &calp1, &salp2, &calp2, &dnm,
747 s12x = sig12 * g->b * dnm;
748 m12x = sq(dnm) * g->b * sin(sig12 / dnm);
749 if (outmask & GEOD_GEODESICSCALE)
750 M12 = M21 = cos(sig12 / dnm);
751 a12 = sig12 / degree;
752 omg12 = lam12 / (g->f1 * dnm);
766 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0;
769 real salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
771 for (tripn = FALSE, tripb = FALSE; numit < maxit2; ++numit) {
775 v = (Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
776 &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
777 &eps, &omg12, numit < maxit1, &dv, C1a, C2a, C3a)
781 if (tripb || !(fabs(v) >= (tripn ? 8 : 2) * tol0))
break;
783 if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))
784 { salp1b = salp1; calp1b = calp1; }
785 else if (v < 0 && (numit > maxit1 || calp1/salp1 < calp1a/salp1a))
786 { salp1a = salp1; calp1a = calp1; }
787 if (numit < maxit1 && dv > 0) {
791 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
792 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
793 if (nsalp1 > 0 && fabs(dalp1) < pi) {
794 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
796 SinCosNorm(&salp1, &calp1);
800 tripn = fabs(v) <= 16 * tol0;
812 salp1 = (salp1a + salp1b)/2;
813 calp1 = (calp1a + calp1b)/2;
814 SinCosNorm(&salp1, &calp1);
816 tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb ||
817 fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb);
821 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
822 cbet1, cbet2, &s12x, &m12x, &dummy,
823 (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
827 a12 = sig12 / degree;
828 omg12 = lam12 - omg12;
832 if (outmask & GEOD_DISTANCE)
835 if (outmask & GEOD_REDUCEDLENGTH)
838 if (outmask & GEOD_AREA) {
841 salp0 = salp1 * cbet1,
842 calp0 = hypotx(calp1, salp1 * sbet1);
844 if (calp0 != 0 && salp0 != 0) {
847 ssig1 = sbet1, csig1 = calp1 * cbet1,
848 ssig2 = sbet2, csig2 = calp2 * cbet2,
849 k2 = sq(calp0) * g->ep2,
850 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
852 A4 = sq(g->
a) * calp0 * salp0 * g->e2;
855 SinCosNorm(&ssig1, &csig1);
856 SinCosNorm(&ssig2, &csig2);
858 B41 = SinCosSeries(FALSE, ssig1, csig1, C4a, nC4);
859 B42 = SinCosSeries(FALSE, ssig2, csig2, C4a, nC4);
860 S12 = A4 * (B42 - B41);
866 omg12 < (
real)(0.75) * pi &&
867 sbet2 - sbet1 < (
real)(1.75)) {
872 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
873 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
874 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
875 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
879 salp12 = salp2 * calp1 - calp2 * salp1,
880 calp12 = calp2 * calp1 + salp2 * salp1;
885 if (salp12 == 0 && calp12 < 0) {
886 salp12 = tiny * calp1;
889 alp12 = atan2(salp12, calp12);
891 S12 += g->c2 * alp12;
892 S12 *= swapp * lonsign * latsign;
899 swapx(&salp1, &salp2);
900 swapx(&calp1, &calp2);
901 if (outmask & GEOD_GEODESICSCALE)
905 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
906 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
908 if (outmask & GEOD_AZIMUTH) {
910 azi1 = 0 - atan2(-salp1, calp1) / degree;
911 azi2 = 0 - atan2(-salp2, calp2) / degree;
914 if (outmask & GEOD_DISTANCE)
916 if (outmask & GEOD_AZIMUTH) {
917 if (pazi1) *pazi1 =
azi1;
918 if (pazi2) *pazi2 = azi2;
920 if (outmask & GEOD_REDUCEDLENGTH)
922 if (outmask & GEOD_GEODESICSCALE) {
923 if (pM12) *pM12 = M12;
924 if (pM21) *pM21 = M21;
926 if (outmask & GEOD_AREA)
936 geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
939 real SinCosSeries(boolx sinp,
real sinx,
real cosx,
const real c[],
int n) {
947 ar = 2 * (cosx - sinx) * (cosx + sinx);
948 y0 = n & 1 ? *--c : 0; y1 = 0;
953 y1 = ar * y0 - y1 + *--c;
954 y0 = ar * y1 - y0 + *--c;
957 ? 2 * sinx * cosx * y0
967 boolx scalep,
real* pM12,
real* pM21,
970 real s12b = 0, m12b = 0, m0 = 0, M12 = 0, M21 = 0;
971 real A1m1, AB1, A2m1, AB2, J12;
978 AB1 = (1 + A1m1) * (SinCosSeries(TRUE, ssig2, csig2, C1a, nC1) -
979 SinCosSeries(TRUE, ssig1, csig1, C1a, nC1));
981 AB2 = (1 + A2m1) * (SinCosSeries(TRUE, ssig2, csig2, C2a, nC2) -
982 SinCosSeries(TRUE, ssig1, csig1, C2a, nC2));
984 J12 = m0 * sig12 + (AB1 - AB2);
988 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12;
990 s12b = (1 + A1m1) * sig12 + AB1;
992 real csig12 = csig1 * csig2 + ssig1 * ssig2;
993 real t = g->ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
994 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
995 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
1013 r = (p + q - 1) / 6;
1014 if ( !(q == 0 && r <= 0) ) {
1023 disc = S * (S + 2 * r3);
1027 real T3 = S + r3, T;
1031 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
1035 u += T + (T != 0 ? r2 / T : 0);
1038 real ang = atan2(sqrt(-disc), -(S + r3));
1041 u += 2 * r * cos(ang / 3);
1043 v = sqrt(sq(u) + q);
1045 uv = u < 0 ? q / (v - u) : u + v;
1046 w = (uv - q) / (2 * v);
1049 k = uv / (sqrt(uv + sq(w)) + w);
1069 real salp1 = 0, calp1 = 0, salp2 = 0, calp2 = 0, dnm = 0;
1077 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
1078 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
1079 #if defined(__GNUC__) && __GNUC__ == 4 && \
1080 (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
1089 volatile real xx1 = sbet2 * cbet1;
1090 volatile real xx2 = cbet2 * sbet1;
1091 sbet12a = xx1 + xx2;
1094 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
1096 boolx shortline = cbet12 >= 0 && sbet12 < (
real)(0.5) &&
1097 cbet2 * lam12 < (
real)(0.5);
1098 real omg12 = lam12, somg12, comg12, ssig12, csig12;
1100 real sbetm2 = sq(sbet1 + sbet2);
1103 sbetm2 /= sbetm2 + sq(cbet1 + cbet2);
1104 dnm = sqrt(1 + g->ep2 * sbetm2);
1105 omg12 /= g->f1 * dnm;
1107 somg12 = sin(omg12); comg12 = cos(omg12);
1109 salp1 = cbet2 * somg12;
1110 calp1 = comg12 >= 0 ?
1111 sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) :
1112 sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1114 ssig12 = hypotx(salp1, calp1);
1115 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
1117 if (shortline && ssig12 < g->etol2) {
1119 salp2 = cbet1 * somg12;
1120 calp2 = sbet12 - cbet1 * sbet2 *
1121 (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);
1122 SinCosNorm(&salp2, &calp2);
1124 sig12 = atan2(ssig12, csig12);
1125 }
else if (fabs(g->n) > (
real)(0.1) ||
1127 ssig12 >= 6 * fabs(g->n) * pi * sq(cbet1)) {
1132 real y, lamscale, betscale;
1141 k2 = sq(sbet1) * g->ep2,
1142 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1143 lamscale = g->
f * cbet1 * A3f(g, eps) * pi;
1145 betscale = lamscale * cbet1;
1147 x = (lam12 - pi) / lamscale;
1148 y = sbet12a / betscale;
1152 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1153 bet12a = atan2(sbet12a, cbet12a);
1154 real m12b, m0, dummy;
1157 Lengths(g, g->n, pi + bet12a,
1158 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1159 cbet1, cbet2, &dummy, &m12b, &m0, FALSE,
1160 &dummy, &dummy, C1a, C2a);
1161 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi);
1162 betscale = x < -(
real)(0.01) ? sbet12a / x :
1163 -g->
f * sq(cbet1) * pi;
1164 lamscale = betscale / cbet1;
1165 y = (lam12 - pi) / lamscale;
1168 if (y > -tol1 && x > -1 - xthresh) {
1171 salp1 = minx((
real)(1), -(
real)(x)); calp1 = - sqrt(1 - sq(salp1));
1173 calp1 = maxx((
real)(x > -tol1 ? 0 : -1), (
real)(x));
1174 salp1 = sqrt(1 - sq(calp1));
1211 real k = Astroid(x, y);
1213 omg12a = lamscale * ( g->
f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
1214 somg12 = sin(omg12a); comg12 = -cos(omg12a);
1216 salp1 = cbet2 * somg12;
1217 calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1222 SinCosNorm(&salp1, &calp1);
1224 salp1 = 1; calp1 = 0;
1247 boolx diffp,
real* pdlam12,
1250 real salp2 = 0, calp2 = 0, sig12 = 0,
1251 ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0, dlam12 = 0;
1253 real somg1, comg1, somg2, comg2, omg12, lam12;
1256 if (sbet1 == 0 && calp1 == 0)
1262 salp0 = salp1 * cbet1;
1263 calp0 = hypotx(calp1, salp1 * sbet1);
1267 ssig1 = sbet1; somg1 = salp0 * sbet1;
1268 csig1 = comg1 = calp1 * cbet1;
1269 SinCosNorm(&ssig1, &csig1);
1276 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
1281 calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ?
1282 sqrt(sq(calp1 * cbet1) +
1284 (cbet2 - cbet1) * (cbet1 + cbet2) :
1285 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1289 ssig2 = sbet2; somg2 = salp0 * sbet2;
1290 csig2 = comg2 = calp2 * cbet2;
1291 SinCosNorm(&ssig2, &csig2);
1295 sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (
real)(0)),
1296 csig1 * csig2 + ssig1 * ssig2);
1299 omg12 = atan2(maxx(comg1 * somg2 - somg1 * comg2, (
real)(0)),
1300 comg1 * comg2 + somg1 * somg2);
1301 k2 = sq(calp0) * g->ep2;
1302 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1304 B312 = (SinCosSeries(TRUE, ssig2, csig2, C3a, nC3-1) -
1305 SinCosSeries(TRUE, ssig1, csig1, C3a, nC3-1));
1306 h0 = -g->
f * A3f(g, eps);
1307 domg12 = salp0 * h0 * (sig12 + B312);
1308 lam12 = omg12 + domg12;
1312 dlam12 = - 2 * g->f1 * dn1 / sbet1;
1315 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1316 cbet1, cbet2, &dummy, &dlam12, &dummy,
1317 FALSE, &dummy, &dummy, C1a, C2a);
1318 dlam12 *= g->f1 / (calp2 * cbet2);
1342 v = eps * v + g->A3x[--i];
1351 for (j = nC3x, k = nC3 - 1; k; ) {
1353 for (i = nC3 - k; i; --i)
1354 t = eps * t + g->C3x[--j];
1358 for (k = 1; k < nC3; ) {
1369 for (j = nC4x, k = nC4; k; ) {
1371 for (i = nC4 - k + 1; i; --i)
1372 t = eps * t + g->C4x[--j];
1376 for (k = 1; k < nC4; ) {
1388 t = eps2*(eps2*(eps2+4)+64)/256;
1389 return (t + eps) / (1 - eps);
1397 c[1] = d*((6-eps2)*eps2-16)/32;
1399 c[2] = d*((64-9*eps2)*eps2-128)/2048;
1401 c[3] = d*(9*eps2-16)/768;
1403 c[4] = d*(3*eps2-5)/512;
1415 c[1] = d*(eps2*(205*eps2-432)+768)/1536;
1417 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
1419 c[3] = d*(116-225*eps2)/384;
1421 c[4] = d*(2695-7173*eps2)/7680;
1425 c[6] = 38081*d/61440;
1432 t = eps2*(eps2*(25*eps2+36)+64)/256;
1433 return t * (1 - eps) - eps;
1441 c[1] = d*(eps2*(eps2+2)+16)/32;
1443 c[2] = d*(eps2*(35*eps2+64)+384)/2048;
1445 c[3] = d*(15*eps2+80)/768;
1447 c[4] = d*(7*eps2+35)/512;
1457 g->A3x[1] = (g->n-1)/2;
1458 g->A3x[2] = (g->n*(3*g->n-1)-2)/8;
1459 g->A3x[3] = ((-g->n-3)*g->n-1)/16;
1460 g->A3x[4] = (-2*g->n-3)/64;
1461 g->A3x[5] = -3/(
real)(128);
1466 g->C3x[0] = (1-g->n)/4;
1467 g->C3x[1] = (1-g->n*g->n)/8;
1468 g->C3x[2] = ((3-g->n)*g->n+3)/64;
1469 g->C3x[3] = (2*g->n+5)/128;
1470 g->C3x[4] = 3/(
real)(128);
1471 g->C3x[5] = ((g->n-3)*g->n+2)/32;
1472 g->C3x[6] = ((-3*g->n-2)*g->n+3)/64;
1473 g->C3x[7] = (g->n+3)/128;
1474 g->C3x[8] = 5/(
real)(256);
1475 g->C3x[9] = (g->n*(5*g->n-9)+5)/192;
1476 g->C3x[10] = (9-10*g->n)/384;
1477 g->C3x[11] = 7/(
real)(512);
1478 g->C3x[12] = (7-14*g->n)/512;
1479 g->C3x[13] = 7/(
real)(512);
1480 g->C3x[14] = 21/(
real)(2560);
1487 g->C4x[0] = (g->n*(g->n*(g->n*(g->n*(100*g->n+208)+572)+3432)-12012)+30030)/
1489 g->C4x[1] = (g->n*(g->n*(g->n*(64*g->n+624)-4576)+6864)-3003)/15015;
1490 g->C4x[2] = (g->n*((14144-10656*g->n)*g->n-4576)-858)/45045;
1491 g->C4x[3] = ((-224*g->n-4784)*g->n+1573)/45045;
1492 g->C4x[4] = (1088*g->n+156)/45045;
1493 g->C4x[5] = 97/(
real)(15015);
1494 g->C4x[6] = (g->n*(g->n*((-64*g->n-624)*g->n+4576)-6864)+3003)/135135;
1495 g->C4x[7] = (g->n*(g->n*(5952*g->n-11648)+9152)-2574)/135135;
1496 g->C4x[8] = (g->n*(5792*g->n+1040)-1287)/135135;
1497 g->C4x[9] = (468-2944*g->n)/135135;
1498 g->C4x[10] = 1/(
real)(9009);
1499 g->C4x[11] = (g->n*((4160-1440*g->n)*g->n-4576)+1716)/225225;
1500 g->C4x[12] = ((4992-8448*g->n)*g->n-1144)/225225;
1501 g->C4x[13] = (1856*g->n-936)/225225;
1502 g->C4x[14] = 8/(
real)(10725);
1503 g->C4x[15] = (g->n*(3584*g->n-3328)+1144)/315315;
1504 g->C4x[16] = (1024*g->n-208)/105105;
1505 g->C4x[17] = -136/(
real)(63063);
1506 g->C4x[18] = (832-2560*g->n)/405405;
1507 g->C4x[19] = -128/(
real)(135135);
1508 g->C4x[20] = 128/(
real)(99099);
1511 int transit(
real lon1,
real lon2) {
1516 lon1 = AngNormalize(lon1);
1517 lon2 = AngNormalize(lon2);
1518 lon12 = AngDiff(lon1, lon2);
1519 return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
1520 (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
1523 int transitdirect(
real lon1,
real lon2) {
1524 lon1 = fmod(lon1, (
real)(720));
1525 lon2 = fmod(lon2, (
real)(720));
1526 return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
1527 ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
1530 void accini(
real s[]) {
1535 void acccopy(
const real s[],
real t[]) {
1537 t[0] = s[0]; t[1] = s[1];
1542 real u, z = sumx(y, s[1], &u);
1543 s[0] = sumx(z, s[0], &s[1]);
1558 void accneg(
real s[]) {
1560 s[0] = -s[0]; s[1] = -s[1];
1564 p->lat0 = p->lon0 = p->
lat = p->
lon = NaN;
1565 p->polyline = (polylinep != 0);
1568 p->
num = p->crossings = 0;
1574 lon = AngNormalize(lon);
1576 p->lat0 = p->
lat = lat;
1577 p->lon0 = p->
lon = lon;
1581 &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
1585 p->crossings += transit(p->
lon, lon);
1587 p->
lat = lat; p->
lon = lon;
1599 0, 0, 0, 0, p->polyline ? 0 : &S12);
1603 p->crossings += transitdirect(p->
lon, lon);
1605 p->
lat = lat; p->
lon = lon;
1612 boolx reverse, boolx sign,
1614 real s12, S12, t[2], area0;
1618 if (!p->polyline && pA) *pA = 0;
1622 if (pP) *pP = p->P[0];
1626 &s12, 0, 0, 0, 0, 0, &S12);
1627 if (pP) *pP = accsum(p->P, s12);
1630 crossings = p->crossings + transit(p->
lon, p->lon0);
1631 area0 = 4 * pi * g->c2;
1633 accadd(t, (t[0] < 0 ? 1 : -1) * area0/2);
1642 else if (t[0] <= -area0/2)
1650 if (pA) *pA = 0 + t[0];
1657 boolx reverse, boolx sign,
1659 real perimeter, tempsum, area0;
1661 unsigned num = p->
num + 1;
1664 if (!p->polyline && pA) *pA = 0;
1667 perimeter = p->P[0];
1668 tempsum = p->polyline ? 0 : p->A[0];
1669 crossings = p->crossings;
1670 for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
1673 i == 0 ? p->
lat : lat, i == 0 ? p->
lon : lon,
1674 i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
1675 &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
1679 crossings += transit(i == 0 ? p->
lon : lon,
1680 i != 0 ? p->lon0 : lon);
1684 if (pP) *pP = perimeter;
1688 area0 = 4 * pi * g->c2;
1690 tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
1697 if (tempsum > area0/2)
1699 else if (tempsum <= -area0/2)
1702 if (tempsum >= area0)
1704 else if (tempsum < 0)
1707 if (pA) *pA = 0 + tempsum;
1714 boolx reverse, boolx sign,
1716 real perimeter, tempsum, area0;
1718 unsigned num = p->
num + 1;
1721 if (!p->polyline && pA) *pA = NaN;
1724 perimeter = p->P[0] + s;
1726 if (pP) *pP = perimeter;
1731 crossings = p->crossings;
1733 real lat, lon, s12, S12;
1738 crossings += transitdirect(p->
lon, lon);
1740 &s12, 0, 0, 0, 0, 0, &S12);
1743 crossings += transit(lon, p->lon0);
1746 area0 = 4 * pi * g->c2;
1748 tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
1755 if (tempsum > area0/2)
1757 else if (tempsum <= -area0/2)
1760 if (tempsum >= area0)
1762 else if (tempsum < 0)
1765 if (pP) *pP = perimeter;
1766 if (pA) *pA = 0 + tempsum;
1776 for (i = 0; i < n; ++i)
unsigned geod_polygon_testedge(const struct geod_geodesic *g, const struct geod_polygon *p, double azi, double s, int reverse, int sign, double *pA, double *pP)
double geod_genposition(const struct geod_geodesicline *l, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
GeographicLib::Math::real real
void geod_polygon_addedge(const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)
void geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
void geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
double geod_geninverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2, double *pm12, double *pM12, double *pM21, double *pS12)
void geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
void geod_polygon_init(struct geod_polygon *p, int polylinep)
void geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
unsigned geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
void geod_polygonarea(const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)
double geod_gendirect(const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
unsigned geod_polygon_testpoint(const struct geod_geodesic *g, const struct geod_polygon *p, double lat, double lon, int reverse, int sign, double *pA, double *pP)
void geod_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
void geod_init(struct geod_geodesic *g, double a, double f)
Header for the geodesic routines in C.