29 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
30 #define nA1 GEOGRAPHICLIB_GEODESIC_ORDER
31 #define nC1 GEOGRAPHICLIB_GEODESIC_ORDER
32 #define nC1p GEOGRAPHICLIB_GEODESIC_ORDER
33 #define nA2 GEOGRAPHICLIB_GEODESIC_ORDER
34 #define nC2 GEOGRAPHICLIB_GEODESIC_ORDER
35 #define nA3 GEOGRAPHICLIB_GEODESIC_ORDER
37 #define nC3 GEOGRAPHICLIB_GEODESIC_ORDER
38 #define nC3x ((nC3 * (nC3 - 1)) / 2)
39 #define nC4 GEOGRAPHICLIB_GEODESIC_ORDER
40 #define nC4x ((nC4 * (nC4 + 1)) / 2)
41 #define nC (GEOGRAPHICLIB_GEODESIC_ORDER + 1)
46 static unsigned init = 0;
47 static const int FALSE = 0;
48 static const int TRUE = 1;
49 static unsigned digits, maxit1, maxit2;
50 static real epsilon, realmin, pi, degree, NaN,
51 tiny, tol0, tol1, tol2, tolb, xthresh;
55 #if defined(__DBL_MANT_DIG__)
56 digits = __DBL_MANT_DIG__;
60 #if defined(__DBL_EPSILON__)
61 epsilon = __DBL_EPSILON__;
63 epsilon = pow(0.5, digits - 1);
65 #if defined(__DBL_MIN__)
66 realmin = __DBL_MIN__;
68 realmin = pow(0.5, 1022);
73 pi = atan2(0.0, -1.0);
76 maxit2 = maxit1 + digits + 10;
86 xthresh = 1000 * tol2;
104 static real sq(
real x) {
return x * x; }
113 return z == 0 ? x : x * log(y) / z;
118 y = log1px(2 * y/(1 - y))/2;
119 return x < 0 ? -y : y;
123 {
return sqrt(x * x + y * y); }
127 return x < 0 ? -y : y;
131 volatile real s = u + v;
132 volatile real up = s - v;
133 volatile real vpp = s - up;
144 real y = N < 0 ? 0 : *p++;
145 while (--N >= 0) y = y * x + *p++;
150 {
return x < y ? x : y; }
153 {
return x > y ? x : y; }
156 {
real t = *x; *x = *y; *y = t; }
158 static void norm2(
real* sinx,
real* cosx) {
159 real r = hypotx(*sinx, *cosx);
165 x = fmod(x, (
real)(360));
166 return x < -180 ? x + 360 : (x < 180 ? x : x - 360);
170 {
return fabs(x) > 90 ? NaN : x; }
173 real t, d = - AngNormalize(sumx(AngNormalize(x), AngNormalize(-y), &t));
180 return (d == 180 && t < 0 ? -180 : d) - t;
185 volatile real y = fabs(x);
187 y = y < z ? z - (z - y) : y;
188 return x < 0 ? 0 - y : y;
195 r = fmod(x, (
real)(360));
196 q = (int)(floor(r / 90 + (
real)(0.5)));
201 s = sin(r); c = cos(r);
202 switch ((
unsigned)q & 3U) {
203 case 0U: *sinx = s; *cosx = c;
break;
204 case 1U: *sinx = c; *cosx = 0 - s;
break;
205 case 2U: *sinx = 0 - s; *cosx = 0 - c;
break;
206 default: *sinx = 0 - c; *cosx = s;
break;
216 if (fabs(y) > fabs(x)) { swapx(&x, &y); q = 2; }
217 if (x < 0) { x = -x; ++q; }
219 ang = atan2(y, x) / degree;
226 case 1: ang = (y > 0 ? 180 : -180) - ang;
break;
227 case 2: ang = 90 - ang;
break;
228 case 3: ang = -90 + ang;
break;
236 static real SinCosSeries(boolx sinp,
238 const real c[],
int n);
269 boolx diffp,
real* pdlam12,
276 static void C1f(
real eps,
real c[]);
277 static void C1pf(
real eps,
real c[]);
279 static void C2f(
real eps,
real c[]);
280 static int transit(
real lon1,
real lon2);
281 static int transitdirect(
real lon1,
real lon2);
282 static void accini(
real s[]);
283 static void acccopy(
const real s[],
real t[]);
284 static void accadd(
real s[],
real y);
286 static void accneg(
real s[]);
291 g->
f = f <= 1 ? f : 1/
f;
293 g->e2 = g->
f * (2 - g->
f);
294 g->ep2 = g->e2 / sq(g->f1);
295 g->n = g->
f / ( 2 - g->
f);
297 g->c2 = (sq(g->
a) + sq(g->b) *
299 (g->e2 > 0 ? atanhx(sqrt(g->e2)) : atan(sqrt(-g->e2))) /
300 sqrt(fabs(g->e2))))/2;
310 g->etol2 = 0.1 * tol2 /
311 sqrt( maxx((
real)(0.001), fabs(g->
f)) * minx((
real)(1), 1 - g->
f/2) / 2 );
321 real cbet1, sbet1, eps;
332 l->
lat1 = LatFix(lat1);
334 l->
azi1 = AngNormalize(azi1);
336 sincosdx(AngRound(l->
azi1), &l->salp1, &l->calp1);
338 sincosdx(AngRound(l->
lat1), &sbet1, &cbet1); sbet1 *= l->f1;
340 norm2(&sbet1, &cbet1); cbet1 = maxx(tiny, cbet1);
341 l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));
344 l->salp0 = l->salp1 * cbet1;
347 l->calp0 = hypotx(l->calp1, l->salp1 * sbet1);
357 l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;
358 l->csig1 = l->comg1 = sbet1 != 0 || l->calp1 != 0 ? cbet1 * l->calp1 : 1;
359 norm2(&l->ssig1, &l->csig1);
362 l->k2 = sq(l->calp0) * g->ep2;
363 eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);
365 if (l->
caps & CAP_C1) {
367 l->A1m1 = A1m1f(eps);
369 l->B11 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C1a, nC1);
370 s = sin(l->B11); c = cos(l->B11);
372 l->stau1 = l->ssig1 * c + l->csig1 * s;
373 l->ctau1 = l->csig1 * c - l->ssig1 * s;
378 if (l->
caps & CAP_C1p)
381 if (l->
caps & CAP_C2) {
382 l->A2m1 = A2m1f(eps);
384 l->B21 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C2a, nC2);
387 if (l->
caps & CAP_C3) {
389 l->A3c = -l->
f * l->salp0 * A3f(g, eps);
390 l->B31 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C3a, nC3-1);
393 if (l->
caps & CAP_C4) {
396 l->A4 = sq(l->
a) * l->calp0 * l->salp0 * g->e2;
397 l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4);
402 unsigned flags,
real s12_a12,
407 real lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,
408 m12 = 0, M12 = 0, M21 = 0, S12 = 0;
410 real sig12, ssig12, csig12, B12 = 0, AB1 = 0;
411 real omg12, lam12, lon12;
412 real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
422 outmask &= l->
caps & OUT_ALL;
430 sig12 = s12_a12 * degree;
431 sincosdx(s12_a12, &ssig12, &csig12);
435 tau12 = s12_a12 / (l->b * (1 + l->A1m1)),
439 B12 = - SinCosSeries(TRUE,
440 l->stau1 * c + l->ctau1 * s,
441 l->ctau1 * c - l->stau1 * s,
443 sig12 = tau12 - (B12 - l->B11);
444 ssig12 = sin(sig12); csig12 = cos(sig12);
445 if (fabs(l->
f) > 0.01) {
468 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12,
469 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12,
471 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
472 serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b;
473 sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2));
474 ssig12 = sin(sig12); csig12 = cos(sig12);
480 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
481 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
482 dn2 = sqrt(1 + l->k2 * sq(ssig2));
484 if (flags & GEOD_ARCMODE || fabs(l->
f) > 0.01)
485 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
486 AB1 = (1 + l->A1m1) * (B12 - l->B11);
489 sbet2 = l->calp0 * ssig2;
491 cbet2 = hypotx(l->salp0, l->calp0 * csig2);
494 cbet2 = csig2 = tiny;
496 salp2 = l->salp0; calp2 = l->calp0 * csig2;
499 s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
502 int E = l->salp0 < 0 ? -1 : 1;
504 somg2 = l->salp0 * ssig2; comg2 = csig2;
508 - (atan2( ssig2, csig2) - atan2( l->ssig1, l->csig1))
509 + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))
510 : atan2(somg2 * l->comg1 - comg2 * l->somg1,
511 comg2 * l->comg1 + somg2 * l->somg1);
512 lam12 = omg12 + l->A3c *
513 ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
515 lon12 = lam12 / degree;
517 AngNormalize(AngNormalize(l->
lon1) + AngNormalize(lon12));
521 lat2 = atan2dx(sbet2, l->f1 * cbet2);
524 azi2 = atan2dx(salp2, calp2);
528 B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),
529 AB2 = (1 + l->A2m1) * (B22 - l->B21),
530 J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);
534 m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
535 - l->csig1 * csig2 * J12);
537 real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / (l->dn1 + dn2);
538 M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1;
539 M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2;
545 B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
547 if (l->calp0 == 0 || l->salp0 == 0) {
549 salp12 = salp2 * l->calp1 - calp2 * l->salp1;
550 calp12 = calp2 * l->calp1 + salp2 * l->salp1;
555 if (salp12 == 0 && calp12 < 0) {
556 salp12 = tiny * l->calp1;
568 salp12 = l->calp0 * l->salp0 *
569 (csig12 <= 0 ? l->csig1 * (1 - csig12) + ssig12 * l->ssig1 :
570 ssig12 * (l->csig1 * ssig12 / (1 + csig12) + l->ssig1));
571 calp12 = sq(l->salp0) + sq(l->calp0) * l->csig1 * csig2;
573 S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
576 if (outmask & GEOD_LATITUDE)
578 if (outmask & GEOD_LONGITUDE)
580 if (outmask & GEOD_AZIMUTH)
582 if (outmask & GEOD_DISTANCE)
584 if (outmask & GEOD_REDUCEDLENGTH)
586 if (outmask & GEOD_GEODESICSCALE) {
587 if (pM12) *pM12 = M12;
588 if (pM21) *pM21 = M21;
590 if (outmask & GEOD_AREA)
593 return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree;
598 geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0);
603 unsigned flags,
real s12_a12,
609 (plat2 ? GEOD_LATITUDE : 0U) |
610 (plon2 ? GEOD_LONGITUDE : 0U) |
611 (pazi2 ? GEOD_AZIMUTH : 0U) |
612 (ps12 ? GEOD_DISTANCE : 0U) |
613 (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
614 (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
615 (pS12 ? GEOD_AREA : 0U);
622 plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
637 real s12 = 0, azi1 = 0, azi2 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
639 int latsign, lonsign, swapp;
640 real sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0;
641 real dn1, dn2, lam12, slam12, clam12;
642 real a12 = 0, sig12, calp1 = 0, salp1 = 0, calp2 = 0, salp2 = 0;
648 (ps12 ? GEOD_DISTANCE : 0U) |
649 (pazi1 || pazi2 ? GEOD_AZIMUTH : 0U) |
650 (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
651 (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
652 (pS12 ? GEOD_AREA : 0U);
659 lon12 = AngRound(AngDiff(lon1, lon2));
661 lonsign = lon12 >= 0 ? 1 : -1;
664 lat1 = AngRound(LatFix(lat1));
665 lat2 = AngRound(LatFix(lat2));
667 swapp = fabs(lat1) >= fabs(lat2) ? 1 : -1;
673 latsign = lat1 < 0 ? 1 : -1;
688 sincosdx(lat1, &sbet1, &cbet1); sbet1 *= g->f1;
690 norm2(&sbet1, &cbet1); cbet1 = maxx(tiny, cbet1);
692 sincosdx(lat2, &sbet2, &cbet2); sbet2 *= g->f1;
694 norm2(&sbet2, &cbet2); cbet2 = maxx(tiny, cbet2);
704 if (cbet1 < -sbet1) {
706 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
708 if (fabs(sbet2) == -sbet1)
712 dn1 = sqrt(1 + g->ep2 * sq(sbet1));
713 dn2 = sqrt(1 + g->ep2 * sq(sbet2));
715 lam12 = lon12 * degree;
716 sincosdx(lon12, &slam12, &clam12);
718 meridian = lat1 == -90 || slam12 == 0;
725 real ssig1, csig1, ssig2, csig2;
726 calp1 = clam12; salp1 = slam12;
727 calp2 = 1; salp2 = 0;
730 ssig1 = sbet1; csig1 = calp1 * cbet1;
731 ssig2 = sbet2; csig2 = calp2 * cbet2;
734 sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (
real)(0)),
735 csig1 * csig2 + ssig1 * ssig2);
736 Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
737 cbet1, cbet2, &s12x, &m12x, 0,
738 outmask & GEOD_GEODESICSCALE ? &M12 : 0,
739 outmask & GEOD_GEODESICSCALE ? &M21 : 0,
748 if (sig12 < 1 || m12x >= 0) {
750 if (sig12 < 3 * tiny)
751 sig12 = m12x = s12x = 0;
754 a12 = sig12 / degree;
763 (g->
f <= 0 || lam12 <= pi - g->f * pi)) {
766 calp1 = calp2 = 0; salp1 = salp2 = 1;
768 sig12 = omg12 = lam12 / g->f1;
769 m12x = g->b * sin(sig12);
770 if (outmask & GEOD_GEODESICSCALE)
771 M12 = M21 = cos(sig12);
774 }
else if (!meridian) {
781 sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
783 &salp1, &calp1, &salp2, &calp2, &dnm,
788 s12x = sig12 * g->b * dnm;
789 m12x = sq(dnm) * g->b * sin(sig12 / dnm);
790 if (outmask & GEOD_GEODESICSCALE)
791 M12 = M21 = cos(sig12 / dnm);
792 a12 = sig12 / degree;
793 omg12 = lam12 / (g->f1 * dnm);
807 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0;
810 real salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
812 for (tripn = FALSE, tripb = FALSE; numit < maxit2; ++numit) {
816 v = (Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
817 &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
818 &eps, &omg12, numit < maxit1, &dv, Ca)
822 if (tripb || !(fabs(v) >= (tripn ? 8 : 2) * tol0))
break;
824 if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))
825 { salp1b = salp1; calp1b = calp1; }
826 else if (v < 0 && (numit > maxit1 || calp1/salp1 < calp1a/salp1a))
827 { salp1a = salp1; calp1a = calp1; }
828 if (numit < maxit1 && dv > 0) {
832 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
833 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
834 if (nsalp1 > 0 && fabs(dalp1) < pi) {
835 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
837 norm2(&salp1, &calp1);
841 tripn = fabs(v) <= 16 * tol0;
853 salp1 = (salp1a + salp1b)/2;
854 calp1 = (calp1a + calp1b)/2;
855 norm2(&salp1, &calp1);
857 tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb ||
858 fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb);
860 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
861 cbet1, cbet2, &s12x, &m12x, 0,
862 outmask & GEOD_GEODESICSCALE ? &M12 : 0,
863 outmask & GEOD_GEODESICSCALE ? &M21 : 0, Ca);
866 a12 = sig12 / degree;
867 omg12 = lam12 - omg12;
871 if (outmask & GEOD_DISTANCE)
874 if (outmask & GEOD_REDUCEDLENGTH)
877 if (outmask & GEOD_AREA) {
880 salp0 = salp1 * cbet1,
881 calp0 = hypotx(calp1, salp1 * sbet1);
883 if (calp0 != 0 && salp0 != 0) {
886 ssig1 = sbet1, csig1 = calp1 * cbet1,
887 ssig2 = sbet2, csig2 = calp2 * cbet2,
888 k2 = sq(calp0) * g->ep2,
889 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
891 A4 = sq(g->
a) * calp0 * salp0 * g->e2;
893 norm2(&ssig1, &csig1);
894 norm2(&ssig2, &csig2);
896 B41 = SinCosSeries(FALSE, ssig1, csig1, Ca, nC4);
897 B42 = SinCosSeries(FALSE, ssig2, csig2, Ca, nC4);
898 S12 = A4 * (B42 - B41);
904 omg12 < (
real)(0.75) * pi &&
905 sbet2 - sbet1 < (
real)(1.75)) {
910 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
911 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
912 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
913 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
917 salp12 = salp2 * calp1 - calp2 * salp1,
918 calp12 = calp2 * calp1 + salp2 * salp1;
923 if (salp12 == 0 && calp12 < 0) {
924 salp12 = tiny * calp1;
927 alp12 = atan2(salp12, calp12);
929 S12 += g->c2 * alp12;
930 S12 *= swapp * lonsign * latsign;
937 swapx(&salp1, &salp2);
938 swapx(&calp1, &calp2);
939 if (outmask & GEOD_GEODESICSCALE)
943 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
944 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
946 if (outmask & GEOD_AZIMUTH) {
948 azi1 = atan2dx(salp1, calp1);
949 azi2 = atan2dx(salp2, calp2);
952 if (outmask & GEOD_DISTANCE)
954 if (outmask & GEOD_AZIMUTH) {
955 if (pazi1) *pazi1 =
azi1;
956 if (pazi2) *pazi2 = azi2;
958 if (outmask & GEOD_REDUCEDLENGTH)
960 if (outmask & GEOD_GEODESICSCALE) {
961 if (pM12) *pM12 = M12;
962 if (pM21) *pM21 = M21;
964 if (outmask & GEOD_AREA)
974 geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
977 real SinCosSeries(boolx sinp,
real sinx,
real cosx,
const real c[],
int n) {
985 ar = 2 * (cosx - sinx) * (cosx + sinx);
986 y0 = n & 1 ? *--c : 0; y1 = 0;
991 y1 = ar * y0 - y1 + *--c;
992 y0 = ar * y1 - y0 + *--c;
995 ? 2 * sinx * cosx * y0
1008 real m0 = 0, J12 = 0, A1 = 0, A2 = 0;
1013 boolx redlp = pm12b || pm0 || pM12 || pM21;
1014 if (ps12b || redlp) {
1026 real B1 = SinCosSeries(TRUE, ssig2, csig2, Ca, nC1) -
1027 SinCosSeries(TRUE, ssig1, csig1, Ca, nC1);
1029 *ps12b = A1 * (sig12 + B1);
1031 real B2 = SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1032 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2);
1033 J12 = m0 * sig12 + (A1 * B1 - A2 * B2);
1038 for (l = 1; l <= nC2; ++l)
1039 Cb[l] = A1 * Ca[l] - A2 * Cb[l];
1040 J12 = m0 * sig12 + (SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1041 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2));
1048 *pm12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1049 csig1 * csig2 * J12;
1051 real csig12 = csig1 * csig2 + ssig1 * ssig2;
1052 real t = g->ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
1054 *pM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
1056 *pM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
1067 r = (p + q - 1) / 6;
1068 if ( !(q == 0 && r <= 0) ) {
1077 disc = S * (S + 2 * r3);
1081 real T3 = S + r3, T;
1085 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
1089 u += T + (T != 0 ? r2 / T : 0);
1092 real ang = atan2(sqrt(-disc), -(S + r3));
1095 u += 2 * r * cos(ang / 3);
1097 v = sqrt(sq(u) + q);
1099 uv = u < 0 ? q / (v - u) : u + v;
1100 w = (uv - q) / (2 * v);
1103 k = uv / (sqrt(uv + sq(w)) + w);
1123 real salp1 = 0, calp1 = 0, salp2 = 0, calp2 = 0, dnm = 0;
1131 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
1132 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
1134 boolx shortline = cbet12 >= 0 && sbet12 < (
real)(0.5) &&
1135 cbet2 * lam12 < (
real)(0.5);
1136 real omg12 = lam12, somg12, comg12, ssig12, csig12;
1137 #if defined(__GNUC__) && __GNUC__ == 4 && \
1138 (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
1146 volatile real xx1 = sbet2 * cbet1;
1147 volatile real xx2 = cbet2 * sbet1;
1148 sbet12a = xx1 + xx2;
1151 sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
1154 real sbetm2 = sq(sbet1 + sbet2);
1157 sbetm2 /= sbetm2 + sq(cbet1 + cbet2);
1158 dnm = sqrt(1 + g->ep2 * sbetm2);
1159 omg12 /= g->f1 * dnm;
1161 somg12 = sin(omg12); comg12 = cos(omg12);
1163 salp1 = cbet2 * somg12;
1164 calp1 = comg12 >= 0 ?
1165 sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) :
1166 sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1168 ssig12 = hypotx(salp1, calp1);
1169 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
1171 if (shortline && ssig12 < g->etol2) {
1173 salp2 = cbet1 * somg12;
1174 calp2 = sbet12 - cbet1 * sbet2 *
1175 (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);
1176 norm2(&salp2, &calp2);
1178 sig12 = atan2(ssig12, csig12);
1179 }
else if (fabs(g->n) > (
real)(0.1) ||
1181 ssig12 >= 6 * fabs(g->n) * pi * sq(cbet1)) {
1186 real y, lamscale, betscale;
1195 k2 = sq(sbet1) * g->ep2,
1196 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1197 lamscale = g->
f * cbet1 * A3f(g, eps) * pi;
1199 betscale = lamscale * cbet1;
1201 x = (lam12 - pi) / lamscale;
1202 y = sbet12a / betscale;
1206 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1207 bet12a = atan2(sbet12a, cbet12a);
1211 Lengths(g, g->n, pi + bet12a,
1212 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1213 cbet1, cbet2, 0, &m12b, &m0, 0, 0, Ca);
1214 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi);
1215 betscale = x < -(
real)(0.01) ? sbet12a / x :
1216 -g->
f * sq(cbet1) * pi;
1217 lamscale = betscale / cbet1;
1218 y = (lam12 - pi) / lamscale;
1221 if (y > -tol1 && x > -1 - xthresh) {
1224 salp1 = minx((
real)(1), -(
real)(x)); calp1 = - sqrt(1 - sq(salp1));
1226 calp1 = maxx((
real)(x > -tol1 ? 0 : -1), (
real)(x));
1227 salp1 = sqrt(1 - sq(calp1));
1264 real k = Astroid(x, y);
1266 omg12a = lamscale * ( g->
f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
1267 somg12 = sin(omg12a); comg12 = -cos(omg12a);
1269 salp1 = cbet2 * somg12;
1270 calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1275 norm2(&salp1, &calp1);
1277 salp1 = 1; calp1 = 0;
1300 boolx diffp,
real* pdlam12,
1303 real salp2 = 0, calp2 = 0, sig12 = 0,
1304 ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0, dlam12 = 0;
1306 real somg1, comg1, somg2, comg2, omg12, lam12;
1309 if (sbet1 == 0 && calp1 == 0)
1315 salp0 = salp1 * cbet1;
1316 calp0 = hypotx(calp1, salp1 * sbet1);
1320 ssig1 = sbet1; somg1 = salp0 * sbet1;
1321 csig1 = comg1 = calp1 * cbet1;
1322 norm2(&ssig1, &csig1);
1329 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
1334 calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ?
1335 sqrt(sq(calp1 * cbet1) +
1337 (cbet2 - cbet1) * (cbet1 + cbet2) :
1338 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1342 ssig2 = sbet2; somg2 = salp0 * sbet2;
1343 csig2 = comg2 = calp2 * cbet2;
1344 norm2(&ssig2, &csig2);
1348 sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (
real)(0)),
1349 csig1 * csig2 + ssig1 * ssig2);
1352 omg12 = atan2(maxx(comg1 * somg2 - somg1 * comg2, (
real)(0)),
1353 comg1 * comg2 + somg1 * somg2);
1354 k2 = sq(calp0) * g->ep2;
1355 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1357 B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) -
1358 SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1));
1359 h0 = -g->
f * A3f(g, eps);
1360 domg12 = salp0 * h0 * (sig12 + B312);
1361 lam12 = omg12 + domg12;
1365 dlam12 = - 2 * g->f1 * dn1 / sbet1;
1367 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1368 cbet1, cbet2, 0, &dlam12, 0, 0, 0, Ca);
1369 dlam12 *= g->f1 / (calp2 * cbet2);
1390 return polyval(nA3 - 1, g->A3x, eps);
1398 for (l = 1; l < nC3; ++l) {
1399 int m = nC3 - l - 1;
1401 c[l] = mult * polyval(m, g->C3x + o, eps);
1411 for (l = 0; l < nC4; ++l) {
1412 int m = nC4 - l - 1;
1413 c[l] = mult * polyval(m, g->C4x + o, eps);
1421 static const real coeff[] = {
1426 real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
1427 return (t + eps) / (1 - eps);
1432 static const real coeff[] = {
1450 for (l = 1; l <= nC1; ++l) {
1451 int m = (nC1 - l) / 2;
1452 c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1460 static const real coeff[] = {
1462 205, -432, 768, 1536,
1464 4005, -4736, 3840, 12288,
1478 for (l = 1; l <= nC1p; ++l) {
1479 int m = (nC1p - l) / 2;
1480 c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1488 static const real coeff[] = {
1490 -11, -28, -192, 0, 256,
1493 real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
1494 return (t - eps) / (1 + eps);
1499 static const real coeff[] = {
1517 for (l = 1; l <= nC2; ++l) {
1518 int m = (nC2 - l) / 2;
1519 c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1527 static const real coeff[] = {
1541 int o = 0, k = 0, j;
1542 for (j = nA3 - 1; j >= 0; --j) {
1543 int m = nA3 - j - 1 < j ? nA3 - j - 1 : j;
1544 g->A3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1551 static const real coeff[] = {
1583 int o = 0, k = 0, l, j;
1584 for (l = 1; l < nC3; ++l) {
1585 for (j = nC3 - 1; j >= l; --j) {
1586 int m = nC3 - j - 1 < j ? nC3 - j - 1 : j;
1587 g->C3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1595 static const real coeff[] = {
1601 -224, -4784, 1573, 45045,
1603 -10656, 14144, -4576, -858, 45045,
1605 64, 624, -4576, 6864, -3003, 15015,
1607 100, 208, 572, 3432, -12012, 30030, 45045,
1613 5792, 1040, -1287, 135135,
1615 5952, -11648, 9152, -2574, 135135,
1617 -64, -624, 4576, -6864, 3003, 135135,
1623 -8448, 4992, -1144, 225225,
1625 -1440, 4160, -4576, 1716, 225225,
1631 3584, -3328, 1144, 315315,
1639 int o = 0, k = 0, l, j;
1640 for (l = 0; l < nC4; ++l) {
1641 for (j = nC4 - 1; j >= l; --j) {
1642 int m = nC4 - j - 1;
1643 g->C4x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1649 int transit(
real lon1,
real lon2) {
1654 lon1 = AngNormalize(lon1);
1655 lon2 = AngNormalize(lon2);
1656 lon12 = AngDiff(lon1, lon2);
1657 return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
1658 (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
1661 int transitdirect(
real lon1,
real lon2) {
1662 lon1 = fmod(lon1, (
real)(720));
1663 lon2 = fmod(lon2, (
real)(720));
1664 return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
1665 ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
1668 void accini(
real s[]) {
1673 void acccopy(
const real s[],
real t[]) {
1675 t[0] = s[0]; t[1] = s[1];
1680 real u, z = sumx(y, s[1], &u);
1681 s[0] = sumx(z, s[0], &s[1]);
1696 void accneg(
real s[]) {
1698 s[0] = -s[0]; s[1] = -s[1];
1702 p->lat0 = p->lon0 = p->
lat = p->
lon = NaN;
1703 p->polyline = (polylinep != 0);
1706 p->
num = p->crossings = 0;
1712 lon = AngNormalize(lon);
1714 p->lat0 = p->
lat = lat;
1715 p->lon0 = p->
lon = lon;
1719 &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
1723 p->crossings += transit(p->
lon, lon);
1725 p->
lat = lat; p->
lon = lon;
1737 0, 0, 0, 0, p->polyline ? 0 : &S12);
1741 p->crossings += transitdirect(p->
lon, lon);
1743 p->
lat = lat; p->
lon = lon;
1750 boolx reverse, boolx sign,
1752 real s12, S12, t[2], area0;
1756 if (!p->polyline && pA) *pA = 0;
1760 if (pP) *pP = p->P[0];
1764 &s12, 0, 0, 0, 0, 0, &S12);
1765 if (pP) *pP = accsum(p->P, s12);
1768 crossings = p->crossings + transit(p->
lon, p->lon0);
1769 area0 = 4 * pi * g->c2;
1771 accadd(t, (t[0] < 0 ? 1 : -1) * area0/2);
1780 else if (t[0] <= -area0/2)
1788 if (pA) *pA = 0 + t[0];
1795 boolx reverse, boolx sign,
1797 real perimeter, tempsum, area0;
1799 unsigned num = p->
num + 1;
1802 if (!p->polyline && pA) *pA = 0;
1805 perimeter = p->P[0];
1806 tempsum = p->polyline ? 0 : p->A[0];
1807 crossings = p->crossings;
1808 for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
1811 i == 0 ? p->
lat : lat, i == 0 ? p->
lon : lon,
1812 i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
1813 &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
1817 crossings += transit(i == 0 ? p->
lon : lon,
1818 i != 0 ? p->lon0 : lon);
1822 if (pP) *pP = perimeter;
1826 area0 = 4 * pi * g->c2;
1828 tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
1835 if (tempsum > area0/2)
1837 else if (tempsum <= -area0/2)
1840 if (tempsum >= area0)
1842 else if (tempsum < 0)
1845 if (pA) *pA = 0 + tempsum;
1852 boolx reverse, boolx sign,
1854 real perimeter, tempsum, area0;
1856 unsigned num = p->
num + 1;
1859 if (!p->polyline && pA) *pA = NaN;
1862 perimeter = p->P[0] + s;
1864 if (pP) *pP = perimeter;
1869 crossings = p->crossings;
1871 real lat, lon, s12, S12;
1876 crossings += transitdirect(p->
lon, lon);
1878 &s12, 0, 0, 0, 0, 0, &S12);
1881 crossings += transit(lon, p->lon0);
1884 area0 = 4 * pi * g->c2;
1886 tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
1893 if (tempsum > area0/2)
1895 else if (tempsum <= -area0/2)
1898 if (tempsum >= area0)
1900 else if (tempsum < 0)
1903 if (pP) *pP = perimeter;
1904 if (pA) *pA = 0 + tempsum;
1914 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.