31 #if !defined(HAVE_C99_MATH) 34 #define HAVE_C99_MATH 1 36 #define HAVE_C99_MATH 0 40 #if !defined(__cplusplus) 44 #define GEOGRAPHICLIB_GEODESIC_ORDER 6 45 #define nA1 GEOGRAPHICLIB_GEODESIC_ORDER 46 #define nC1 GEOGRAPHICLIB_GEODESIC_ORDER 47 #define nC1p GEOGRAPHICLIB_GEODESIC_ORDER 48 #define nA2 GEOGRAPHICLIB_GEODESIC_ORDER 49 #define nC2 GEOGRAPHICLIB_GEODESIC_ORDER 50 #define nA3 GEOGRAPHICLIB_GEODESIC_ORDER 52 #define nC3 GEOGRAPHICLIB_GEODESIC_ORDER 53 #define nC3x ((nC3 * (nC3 - 1)) / 2) 54 #define nC4 GEOGRAPHICLIB_GEODESIC_ORDER 55 #define nC4x ((nC4 * (nC4 + 1)) / 2) 56 #define nC (GEOGRAPHICLIB_GEODESIC_ORDER + 1) 61 static unsigned init = 0;
62 static const int FALSE = 0;
63 static const int TRUE = 1;
64 static unsigned digits, maxit1, maxit2;
65 static real epsilon, realmin, pi, degree, NaN,
66 tiny, tol0, tol1, tol2, tolb, xthresh;
70 digits = DBL_MANT_DIG;
71 epsilon = DBL_EPSILON;
76 pi = atan2(0.0, -1.0);
79 maxit2 = maxit1 + digits + 10;
89 xthresh = 1000 * tol2;
119 #define copysignx copysign 121 #define remainderx remainder 122 #define remquox remquo 127 x = fabs(x); y = fabs(y);
130 return y * sqrt(1 + x * x);
132 y /= (x != 0 ? x : 1);
133 return x * sqrt(1 + y * y);
145 return z == 0 ? x : x * log(y) / z;
150 y = log1px(2 * y/(1 - y))/2;
151 return x > 0 ? y : (x < 0 ? -y : x);
156 return fabs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
161 return x > 0 ? y : (x < 0 ? -y : x);
174 else if (2 * fabs(z) == y)
175 z -= fmod(x, 2 * y) - z;
176 else if (2 * fabs(z) > y)
177 z += (z < 0 ? y : -y);
182 real z = remainderx(x, y);
185 a = remainderx(x, 2 * y),
186 b = remainderx(x, 4 * y),
187 c = remainderx(x, 8 * y);
188 *n = (
a > z ? 1 : (
a < z ? -1 : 0));
189 *n += (b >
a ? 2 : (b <
a ? -2 : 0));
190 *n += (c > b ? 4 : (c < b ? -4 : 0));
193 if (x/y > 0 && *n <= 0)
195 else if (x/y < 0 && *n >= 0)
204 static real sq(
real x) {
return x * x; }
207 volatile real s = u + v;
208 volatile real up = s - v;
209 volatile real vpp = s - up;
212 if (t) *t = -(up + vpp);
220 real y = N < 0 ? 0 : *p++;
221 while (--N >= 0) y = y * x + *p++;
227 {
return (b <
a) ? b :
a; }
230 {
return (
a < b) ? b :
a; }
233 {
real t = *x; *x = *y; *y = t; }
235 static void norm2(
real* sinx,
real* cosx) {
236 real r = hypotx(*sinx, *cosx);
242 x = remainderx(x, (
real)(360));
243 return x != -180 ? x : 180;
247 {
return fabs(x) > 90 ? NaN : x; }
250 real t, d = AngNormalize(sumx(AngNormalize(-x), AngNormalize(y), &t));
257 return sumx(d == 180 && t > 0 ? -180 : d, t, e);
263 if (x == 0)
return 0;
266 y = y < z ? z - (z - y) : y;
267 return x < 0 ? -y : y;
274 r = remquox(x, (
real)(90), &q);
278 s = sin(r); c = cos(r);
279 #if defined(_MSC_VER) && _MSC_VER < 1900 290 switch ((
unsigned)q & 3U) {
291 case 0U: *sinx = s; *cosx = c;
break;
292 case 1U: *sinx = c; *cosx = -s;
break;
293 case 2U: *sinx = -s; *cosx = -c;
break;
294 default: *sinx = -c; *cosx = s;
break;
296 if (x != 0) { *sinx += (
real)(0); *cosx += (
real)(0); }
305 if (fabs(y) > fabs(x)) { swapx(&x, &y); q = 2; }
306 if (x < 0) { x = -x; ++q; }
308 ang = atan2(y, x) / degree;
315 case 1: ang = (y >= 0 ? 180 : -180) - ang;
break;
316 case 2: ang = 90 - ang;
break;
317 case 3: ang = -90 + ang;
break;
325 static real SinCosSeries(boolx sinp,
327 const real c[],
int n);
360 boolx diffp,
real* pdlam12,
367 static void C1f(
real eps,
real c[]);
368 static void C1pf(
real eps,
real c[]);
370 static void C2f(
real eps,
real c[]);
371 static int transit(
real lon1,
real lon2);
372 static int transitdirect(
real lon1,
real lon2);
373 static void accini(
real s[]);
374 static void acccopy(
const real s[],
real t[]);
375 static void accadd(
real s[],
real y);
377 static void accneg(
real s[]);
378 static void accrem(
real s[],
real y);
380 int crossings, boolx reverse, boolx sign);
382 int crossings, boolx reverse, boolx sign);
389 g->e2 = g->
f * (2 - g->
f);
390 g->ep2 = g->e2 / sq(g->f1);
391 g->n = g->
f / ( 2 - g->
f);
393 g->c2 = (sq(g->
a) + sq(g->b) *
395 (g->e2 > 0 ? atanhx(sqrt(g->e2)) : atan(sqrt(-g->e2))) /
396 sqrt(fabs(g->e2))))/2;
406 g->etol2 = 0.1 * tol2 /
407 sqrt( maxx((
real)(0.001), fabs(g->
f)) * minx((
real)(1), 1 - g->
f/2) / 2 );
419 real cbet1, sbet1, eps;
430 l->
lat1 = LatFix(lat1);
436 sincosdx(AngRound(l->
lat1), &sbet1, &cbet1); sbet1 *= l->f1;
438 norm2(&sbet1, &cbet1); cbet1 = maxx(tiny, cbet1);
439 l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));
442 l->salp0 = l->
salp1 * cbet1;
445 l->calp0 = hypotx(l->
calp1, l->
salp1 * sbet1);
455 l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;
456 l->csig1 = l->comg1 = sbet1 != 0 || l->
calp1 != 0 ? cbet1 * l->
calp1 : 1;
457 norm2(&l->ssig1, &l->csig1);
460 l->k2 = sq(l->calp0) * g->ep2;
461 eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);
463 if (l->
caps & CAP_C1) {
465 l->A1m1 = A1m1f(eps);
467 l->B11 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C1a, nC1);
468 s = sin(l->B11); c = cos(l->B11);
470 l->stau1 = l->ssig1 * c + l->csig1 * s;
471 l->ctau1 = l->csig1 * c - l->ssig1 * s;
476 if (l->
caps & CAP_C1p)
479 if (l->
caps & CAP_C2) {
480 l->A2m1 = A2m1f(eps);
482 l->B21 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C2a, nC2);
485 if (l->
caps & CAP_C3) {
487 l->A3c = -l->
f * l->salp0 * A3f(g, eps);
488 l->B31 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C3a, nC3-1);
491 if (l->
caps & CAP_C4) {
494 l->A4 = sq(l->
a) * l->calp0 * l->salp0 * g->e2;
495 l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4);
505 azi1 = AngNormalize(azi1);
507 sincosdx(AngRound(azi1), &salp1, &calp1);
508 geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
514 unsigned flags,
real s12_a12,
523 real s12,
unsigned caps) {
528 unsigned flags,
real s12_a12,
533 real lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,
534 m12 = 0, M12 = 0, M21 = 0, S12 = 0;
536 real sig12, ssig12, csig12, B12 = 0, AB1 = 0;
537 real omg12, lam12, lon12;
538 real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
548 outmask &= l->
caps & OUT_ALL;
556 sig12 = s12_a12 * degree;
557 sincosdx(s12_a12, &ssig12, &csig12);
561 tau12 = s12_a12 / (l->b * (1 + l->A1m1)),
565 B12 = - SinCosSeries(TRUE,
566 l->stau1 * c + l->ctau1 * s,
567 l->ctau1 * c - l->stau1 * s,
569 sig12 = tau12 - (B12 - l->B11);
570 ssig12 = sin(sig12); csig12 = cos(sig12);
571 if (fabs(l->
f) > 0.01) {
594 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
595 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
596 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
597 serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b;
598 sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2));
599 ssig12 = sin(sig12); csig12 = cos(sig12);
605 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
606 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
607 dn2 = sqrt(1 + l->k2 * sq(ssig2));
610 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
611 AB1 = (1 + l->A1m1) * (B12 - l->B11);
614 sbet2 = l->calp0 * ssig2;
616 cbet2 = hypotx(l->salp0, l->calp0 * csig2);
619 cbet2 = csig2 = tiny;
621 salp2 = l->salp0; calp2 = l->calp0 * csig2;
625 l->b * ((1 + l->A1m1) * sig12 + AB1) :
629 real E = copysignx(1, l->salp0);
631 somg2 = l->salp0 * ssig2; comg2 = csig2;
635 - (atan2( ssig2, csig2) - atan2( l->ssig1, l->csig1))
636 + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))
637 : atan2(somg2 * l->comg1 - comg2 * l->somg1,
638 comg2 * l->comg1 + somg2 * l->somg1);
639 lam12 = omg12 + l->A3c *
640 ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
642 lon12 = lam12 / degree;
644 AngNormalize(AngNormalize(l->
lon1) + AngNormalize(lon12));
648 lat2 = atan2dx(sbet2, l->f1 * cbet2);
651 azi2 = atan2dx(salp2, calp2);
655 B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),
656 AB2 = (1 + l->A2m1) * (B22 - l->B21),
657 J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);
661 m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
662 - l->csig1 * csig2 * J12);
664 real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) /
666 M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1;
667 M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2;
673 B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
675 if (l->calp0 == 0 || l->salp0 == 0) {
688 salp12 = l->calp0 * l->salp0 *
689 (csig12 <= 0 ? l->csig1 * (1 - csig12) + ssig12 * l->ssig1 :
690 ssig12 * (l->csig1 * ssig12 / (1 + csig12) + l->ssig1));
691 calp12 = sq(l->salp0) + sq(l->calp0) * l->csig1 * csig2;
693 S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
715 if (pM12) *pM12 = M12;
716 if (pM21) *pM21 = M21;
721 return (flags &
GEOD_ARCMODE) ? s12_a12 : sig12 / degree;
727 nullptr,
nullptr,
nullptr,
nullptr,
nullptr);
731 l->
a13 = a13; l->
s13 = NaN;
733 nullptr,
nullptr,
nullptr,
nullptr);
737 unsigned flags,
real s13_a13) {
739 geod_setarc(l, s13_a13) :
746 nullptr,
nullptr,
nullptr,
nullptr,
nullptr);
751 unsigned flags,
real s12_a12,
770 plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
778 nullptr,
nullptr,
nullptr,
nullptr,
nullptr);
788 real s12 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
790 int latsign, lonsign, swapp;
791 real sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0;
792 real dn1, dn2, lam12, slam12, clam12;
797 real omg12 = 0, somg12 = 2, comg12 = 0;
809 lon12 = AngDiff(
lon1, lon2, &lon12s);
811 lonsign = lon12 >= 0 ? 1 : -1;
813 lon12 = lonsign * AngRound(lon12);
814 lon12s = AngRound((180 - lon12) - lonsign * lon12s);
815 lam12 = lon12 * degree;
817 sincosdx(lon12s, &slam12, &clam12);
820 sincosdx(lon12, &slam12, &clam12);
824 lat2 = AngRound(LatFix(lat2));
827 swapp = fabs(
lat1) < fabs(lat2) ? -1 : 1;
833 latsign =
lat1 < 0 ? 1 : -1;
848 sincosdx(
lat1, &sbet1, &cbet1); sbet1 *= g->f1;
850 norm2(&sbet1, &cbet1); cbet1 = maxx(tiny, cbet1);
852 sincosdx(lat2, &sbet2, &cbet2); sbet2 *= g->f1;
854 norm2(&sbet2, &cbet2); cbet2 = maxx(tiny, cbet2);
864 if (cbet1 < -sbet1) {
866 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
868 if (fabs(sbet2) == -sbet1)
872 dn1 = sqrt(1 + g->ep2 * sq(sbet1));
873 dn2 = sqrt(1 + g->ep2 * sq(sbet2));
875 meridian =
lat1 == -90 || slam12 == 0;
882 real ssig1, csig1, ssig2, csig2;
884 calp2 = 1; salp2 = 0;
887 ssig1 = sbet1; csig1 =
calp1 * cbet1;
888 ssig2 = sbet2; csig2 = calp2 * cbet2;
891 sig12 = atan2(maxx((
real)(0), csig1 * ssig2 - ssig1 * csig2),
892 csig1 * csig2 + ssig1 * ssig2);
893 Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
894 cbet1, cbet2, &s12x, &m12x,
nullptr,
905 if (sig12 < 1 || m12x >= 0) {
907 if (sig12 < 3 * tiny)
908 sig12 = m12x = s12x = 0;
911 a12 = sig12 / degree;
920 (g->
f <= 0 || lon12s >= g->
f * 180)) {
925 sig12 = omg12 = lam12 / g->f1;
926 m12x = g->b * sin(sig12);
928 M12 = M21 = cos(sig12);
931 }
else if (!meridian) {
938 sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
939 lam12, slam12, clam12,
945 s12x = sig12 * g->b * dnm;
946 m12x = sq(dnm) * g->b * sin(sig12 / dnm);
948 M12 = M21 = cos(sig12 / dnm);
949 a12 = sig12 / degree;
950 omg12 = lam12 / (g->f1 * dnm);
964 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0;
967 real salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
970 for (; numit < maxit2; ++numit) {
974 v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
salp1,
calp1,
976 &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
977 &eps, &domg12, numit < maxit1, &dv, Ca);
980 if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0))
break;
982 if (v > 0 && (numit > maxit1 ||
calp1/
salp1 > calp1b/salp1b))
984 else if (v < 0 && (numit > maxit1 ||
calp1/
salp1 < calp1a/salp1a))
986 if (numit < maxit1 && dv > 0) {
990 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
992 if (nsalp1 > 0 && fabs(dalp1) < pi) {
999 tripn = fabs(v) <= 16 * tol0;
1011 salp1 = (salp1a + salp1b)/2;
1012 calp1 = (calp1a + calp1b)/2;
1015 tripb = (fabs(salp1a -
salp1) + (calp1a -
calp1) < tolb ||
1016 fabs(
salp1 - salp1b) + (
calp1 - calp1b) < tolb);
1018 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1019 cbet1, cbet2, &s12x, &m12x,
nullptr,
1024 a12 = sig12 / degree;
1027 real sdomg12 = sin(domg12), cdomg12 = cos(domg12);
1028 somg12 = slam12 * cdomg12 - clam12 * sdomg12;
1029 comg12 = clam12 * cdomg12 + slam12 * sdomg12;
1043 salp0 =
salp1 * cbet1,
1046 if (calp0 != 0 && salp0 != 0) {
1049 ssig1 = sbet1, csig1 =
calp1 * cbet1,
1050 ssig2 = sbet2, csig2 = calp2 * cbet2,
1051 k2 = sq(calp0) * g->ep2,
1052 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
1054 A4 = sq(g->
a) * calp0 * salp0 * g->e2;
1056 norm2(&ssig1, &csig1);
1057 norm2(&ssig2, &csig2);
1059 B41 = SinCosSeries(FALSE, ssig1, csig1, Ca, nC4);
1060 B42 = SinCosSeries(FALSE, ssig2, csig2, Ca, nC4);
1061 S12 = A4 * (B42 - B41);
1066 if (!meridian && somg12 > 1) {
1067 somg12 = sin(omg12); comg12 = cos(omg12);
1072 comg12 > -(
real)(0.7071) &&
1073 sbet2 - sbet1 < (
real)(1.75)) {
1078 domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
1079 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
1080 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
1090 if (salp12 == 0 && calp12 < 0) {
1091 salp12 = tiny *
calp1;
1094 alp12 = atan2(salp12, calp12);
1096 S12 += g->c2 * alp12;
1097 S12 *= swapp * lonsign * latsign;
1104 swapx(&
salp1, &salp2);
1105 swapx(&
calp1, &calp2);
1110 salp1 *= swapp * lonsign;
calp1 *= swapp * latsign;
1111 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
1113 if (psalp1) *psalp1 =
salp1;
1114 if (pcalp1) *pcalp1 =
calp1;
1115 if (psalp2) *psalp2 = salp2;
1116 if (pcalp2) *pcalp2 = calp2;
1123 if (pM12) *pM12 = M12;
1124 if (pM21) *pM21 = M21;
1138 a12 = geod_geninverse_int(g,
lat1,
lon1, lat2, lon2, ps12,
1140 pm12, pM12, pM21, pS12);
1142 if (pazi2) *pazi2 = atan2dx(salp2, calp2);
1151 a12 = geod_geninverse_int(g,
lat1,
lon1, lat2, lon2,
nullptr,
1153 nullptr,
nullptr,
nullptr,
nullptr),
1159 geod_setarc(l, a12);
1166 nullptr,
nullptr,
nullptr,
nullptr);
1169 real SinCosSeries(boolx sinp,
real sinx,
real cosx,
const real c[],
int n) {
1177 ar = 2 * (cosx - sinx) * (cosx + sinx);
1178 y0 = (n & 1) ? *--c : 0; y1 = 0;
1183 y1 = ar * y0 - y1 + *--c;
1184 y0 = ar * y1 - y0 + *--c;
1187 ? 2 * sinx * cosx * y0
1200 real m0 = 0, J12 = 0, A1 = 0, A2 = 0;
1205 boolx redlp = pm12b || pm0 || pM12 || pM21;
1206 if (ps12b || redlp) {
1218 real B1 = SinCosSeries(TRUE, ssig2, csig2, Ca, nC1) -
1219 SinCosSeries(TRUE, ssig1, csig1, Ca, nC1);
1221 *ps12b = A1 * (sig12 + B1);
1223 real B2 = SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1224 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2);
1225 J12 = m0 * sig12 + (A1 * B1 - A2 * B2);
1230 for (l = 1; l <= nC2; ++l)
1231 Cb[l] = A1 * Ca[l] - A2 * Cb[l];
1232 J12 = m0 * sig12 + (SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1233 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2));
1240 *pm12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1241 csig1 * csig2 * J12;
1243 real csig12 = csig1 * csig2 + ssig1 * ssig2;
1244 real t = g->ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
1246 *pM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
1248 *pM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
1259 r = (p + q - 1) / 6;
1260 if ( !(q == 0 && r <= 0) ) {
1269 disc = S * (S + 2 * r3);
1273 real T3 = S + r3, T;
1277 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
1281 u += T + (T != 0 ? r2 / T : 0);
1284 real ang = atan2(sqrt(-disc), -(S + r3));
1287 u += 2 * r * cos(ang / 3);
1289 v = sqrt(sq(u) + q);
1291 uv = u < 0 ? q / (v - u) : u + v;
1292 w = (uv - q) / (2 * v);
1295 k = uv / (sqrt(uv + sq(w)) + w);
1323 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
1324 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
1326 boolx shortline = cbet12 >= 0 && sbet12 < (
real)(0.5) &&
1327 cbet2 * lam12 < (
real)(0.5);
1328 real somg12, comg12, ssig12, csig12;
1329 sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
1331 real sbetm2 = sq(sbet1 + sbet2), omg12;
1334 sbetm2 /= sbetm2 + sq(cbet1 + cbet2);
1335 dnm = sqrt(1 + g->ep2 * sbetm2);
1336 omg12 = lam12 / (g->f1 * dnm);
1337 somg12 = sin(omg12); comg12 = cos(omg12);
1339 somg12 = slam12; comg12 = clam12;
1342 salp1 = cbet2 * somg12;
1343 calp1 = comg12 >= 0 ?
1344 sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) :
1345 sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1348 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
1350 if (shortline && ssig12 < g->etol2) {
1352 salp2 = cbet1 * somg12;
1353 calp2 = sbet12 - cbet1 * sbet2 *
1354 (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);
1355 norm2(&salp2, &calp2);
1357 sig12 = atan2(ssig12, csig12);
1358 }
else if (fabs(g->n) > (
real)(0.1) ||
1360 ssig12 >= 6 * fabs(g->n) * pi * sq(cbet1)) {
1365 real y, lamscale, betscale;
1370 real lam12x = atan2(-slam12, -clam12);
1375 k2 = sq(sbet1) * g->ep2,
1376 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1377 lamscale = g->
f * cbet1 * A3f(g, eps) * pi;
1379 betscale = lamscale * cbet1;
1381 x = lam12x / lamscale;
1382 y = sbet12a / betscale;
1386 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1387 bet12a = atan2(sbet12a, cbet12a);
1391 Lengths(g, g->n, pi + bet12a,
1392 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1393 cbet1, cbet2,
nullptr, &m12b, &m0,
nullptr,
nullptr, Ca);
1394 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi);
1395 betscale = x < -(
real)(0.01) ? sbet12a / x :
1396 -g->
f * sq(cbet1) * pi;
1397 lamscale = betscale / cbet1;
1398 y = lam12x / lamscale;
1401 if (y > -tol1 && x > -1 - xthresh) {
1444 real k = Astroid(x, y);
1446 omg12a = lamscale * ( g->
f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
1447 somg12 = sin(omg12a); comg12 = -cos(omg12a);
1449 salp1 = cbet2 * somg12;
1450 calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1482 boolx diffp,
real* pdlam12,
1485 real salp2 = 0, calp2 = 0, sig12 = 0,
1486 ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0,
1487 domg12 = 0, dlam12 = 0;
1489 real somg1, comg1, somg2, comg2, somg12, comg12, lam12;
1492 if (sbet1 == 0 &&
calp1 == 0)
1498 salp0 =
salp1 * cbet1;
1503 ssig1 = sbet1; somg1 = salp0 * sbet1;
1504 csig1 = comg1 =
calp1 * cbet1;
1505 norm2(&ssig1, &csig1);
1512 salp2 = cbet2 != cbet1 ? salp0 / cbet2 :
salp1;
1517 calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ?
1518 sqrt(sq(
calp1 * cbet1) +
1520 (cbet2 - cbet1) * (cbet1 + cbet2) :
1521 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1525 ssig2 = sbet2; somg2 = salp0 * sbet2;
1526 csig2 = comg2 = calp2 * cbet2;
1527 norm2(&ssig2, &csig2);
1531 sig12 = atan2(maxx((
real)(0), csig1 * ssig2 - ssig1 * csig2),
1532 csig1 * csig2 + ssig1 * ssig2);
1535 somg12 = maxx((
real)(0), comg1 * somg2 - somg1 * comg2);
1536 comg12 = comg1 * comg2 + somg1 * somg2;
1538 eta = atan2(somg12 * clam120 - comg12 * slam120,
1539 comg12 * clam120 + somg12 * slam120);
1540 k2 = sq(calp0) * g->ep2;
1541 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1543 B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) -
1544 SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1));
1545 domg12 = -g->
f * A3f(g, eps) * salp0 * (sig12 + B312);
1546 lam12 = eta + domg12;
1550 dlam12 = - 2 * g->f1 * dn1 / sbet1;
1552 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1553 cbet1, cbet2,
nullptr, &dlam12,
nullptr,
nullptr,
nullptr, Ca);
1554 dlam12 *= g->f1 / (calp2 * cbet2);
1575 return polyval(nA3 - 1, g->A3x, eps);
1583 for (l = 1; l < nC3; ++l) {
1584 int m = nC3 - l - 1;
1586 c[l] = mult * polyval(m, g->C3x + o, eps);
1596 for (l = 0; l < nC4; ++l) {
1597 int m = nC4 - l - 1;
1598 c[l] = mult * polyval(m, g->C4x + o, eps);
1606 static const real coeff[] = {
1611 real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
1612 return (t + eps) / (1 - eps);
1617 static const real coeff[] = {
1635 for (l = 1; l <= nC1; ++l) {
1636 int m = (nC1 - l) / 2;
1637 c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1645 static const real coeff[] = {
1647 205, -432, 768, 1536,
1649 4005, -4736, 3840, 12288,
1663 for (l = 1; l <= nC1p; ++l) {
1664 int m = (nC1p - l) / 2;
1665 c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1673 static const real coeff[] = {
1675 -11, -28, -192, 0, 256,
1678 real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
1679 return (t - eps) / (1 + eps);
1684 static const real coeff[] = {
1702 for (l = 1; l <= nC2; ++l) {
1703 int m = (nC2 - l) / 2;
1704 c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1712 static const real coeff[] = {
1726 int o = 0, k = 0, j;
1727 for (j = nA3 - 1; j >= 0; --j) {
1728 int m = nA3 - j - 1 < j ? nA3 - j - 1 : j;
1729 g->A3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1736 static const real coeff[] = {
1768 int o = 0, k = 0, l, j;
1769 for (l = 1; l < nC3; ++l) {
1770 for (j = nC3 - 1; j >= l; --j) {
1771 int m = nC3 - j - 1 < j ? nC3 - j - 1 : j;
1772 g->C3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1780 static const real coeff[] = {
1786 -224, -4784, 1573, 45045,
1788 -10656, 14144, -4576, -858, 45045,
1790 64, 624, -4576, 6864, -3003, 15015,
1792 100, 208, 572, 3432, -12012, 30030, 45045,
1798 5792, 1040, -1287, 135135,
1800 5952, -11648, 9152, -2574, 135135,
1802 -64, -624, 4576, -6864, 3003, 135135,
1808 -8448, 4992, -1144, 225225,
1810 -1440, 4160, -4576, 1716, 225225,
1816 3584, -3328, 1144, 315315,
1824 int o = 0, k = 0, l, j;
1825 for (l = 0; l < nC4; ++l) {
1826 for (j = nC4 - 1; j >= l; --j) {
1827 int m = nC4 - j - 1;
1828 g->C4x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1840 lon2 = AngNormalize(lon2);
1841 lon12 = AngDiff(
lon1, lon2,
nullptr);
1842 return lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
1843 (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
1850 lon2 = remainderx(lon2, (
real)(720));
1851 return ( (lon2 <= 0 && lon2 > -360 ? 1 : 0) -
1852 (lon1 <= 0 && lon1 > -360 ? 1 : 0) );
1855 void accini(
real s[]) {
1860 void acccopy(
const real s[],
real t[]) {
1862 t[0] = s[0]; t[1] = s[1];
1867 real u, z = sumx(y, s[1], &u);
1868 s[0] = sumx(z, s[0], &s[1]);
1883 void accneg(
real s[]) {
1885 s[0] = -s[0]; s[1] = -s[1];
1890 s[0] = remainderx(s[0], y);
1891 accadd(s, (
real)(0));
1895 p->polyline = (polylinep != 0);
1900 p->lat0 = p->lon0 = p->
lat = p->
lon = NaN;
1903 p->
num = p->crossings = 0;
1909 lon = AngNormalize(lon);
1911 p->lat0 = p->
lat = lat;
1912 p->lon0 = p->
lon = lon;
1916 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
1917 p->polyline ?
nullptr : &S12);
1921 p->crossings += transit(p->
lon, lon);
1923 p->
lat = lat; p->
lon = lon;
1934 real lat = 0, lon = 0, S12 = 0;
1936 &lat, &lon,
nullptr,
1937 nullptr,
nullptr,
nullptr,
nullptr,
1938 p->polyline ?
nullptr : &S12);
1942 p->crossings += transitdirect(p->
lon, lon);
1944 p->
lat = lat; p->
lon = lon;
1951 boolx reverse, boolx sign,
1953 real s12, S12, t[2];
1956 if (!p->polyline && pA) *pA = 0;
1960 if (pP) *pP = p->P[0];
1964 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr, &S12);
1965 if (pP) *pP = accsum(p->P, s12);
1968 if (pA) *pA = areareduceA(t, 4 * pi * g->c2,
1969 p->crossings + transit(p->
lon, p->lon0),
1977 boolx reverse, boolx sign,
1979 real perimeter, tempsum;
1981 unsigned num = p->
num + 1;
1984 if (!p->polyline && pA) *pA = 0;
1987 perimeter = p->P[0];
1988 tempsum = p->polyline ? 0 : p->A[0];
1989 crossings = p->crossings;
1990 for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
1993 i == 0 ? p->
lat : lat, i == 0 ? p->
lon : lon,
1994 i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
1995 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
1996 p->polyline ?
nullptr : &S12);
2000 crossings += transit(i == 0 ? p->
lon : lon,
2001 i != 0 ? p->lon0 : lon);
2005 if (pP) *pP = perimeter;
2009 if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
2016 boolx reverse, boolx sign,
2018 real perimeter, tempsum;
2020 unsigned num = p->
num + 1;
2023 if (!p->polyline && pA) *pA = NaN;
2026 perimeter = p->P[0] + s;
2028 if (pP) *pP = perimeter;
2033 crossings = p->crossings;
2037 real lat = 0, lon = 0, s12, S12 = 0;
2039 &lat, &lon,
nullptr,
2040 nullptr,
nullptr,
nullptr,
nullptr, &S12);
2042 crossings += transitdirect(p->
lon, lon);
2044 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr, &S12);
2047 crossings += transit(lon, p->lon0);
2050 if (pP) *pP = perimeter;
2051 if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
2061 for (i = 0; i < n; ++i)
2067 int crossings, boolx reverse, boolx sign) {
2068 accrem(area, area0);
2070 accadd(area, (area[0] < 0 ? 1 : -1) * area0/2);
2077 if (area[0] > area0/2)
2078 accadd(area, -area0);
2079 else if (area[0] <= -area0/2)
2080 accadd(area, +area0);
2082 if (area[0] >= area0)
2083 accadd(area, -area0);
2084 else if (area[0] < 0)
2085 accadd(area, +area0);
2091 int crossings, boolx reverse, boolx sign) {
2092 area = remainderx(area, area0);
2094 area += (area < 0 ? 1 : -1) * area0/2;
2103 else if (area <= -area0/2)
void GEOD_DLL geod_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
GeographicLib::Math::real real
void GEOD_DLL geod_gendirectline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned flags, double s12_a12, unsigned caps)
void GEOD_DLL geod_polygon_addedge(const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)
unsigned GEOD_DLL geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
void GEOD_DLL geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
void GEOD_DLL geod_directline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, unsigned caps)
void GEOD_DLL geod_inverseline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, unsigned caps)
void GEOD_DLL geod_polygon_clear(struct geod_polygon *p)
double GEOD_DLL 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)
double GEOD_DLL 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)
unsigned GEOD_DLL 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)
double GEOD_DLL 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_DLL geod_gensetdistance(struct geod_geodesicline *l, unsigned flags, double s13_a13)
void GEOD_DLL geod_setdistance(struct geod_geodesicline *l, double s13)
void GEOD_DLL geod_polygonarea(const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)
void GEOD_DLL geod_polygon_init(struct geod_polygon *p, int polylinep)
void GEOD_DLL geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
void GEOD_DLL geod_init(struct geod_geodesic *g, double a, double f)
void GEOD_DLL geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
unsigned GEOD_DLL 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)
void GEOD_DLL geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
API for the geodesic routines in C.