30 #if !defined(__cplusplus)
34 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
35 #define nA1 GEOGRAPHICLIB_GEODESIC_ORDER
36 #define nC1 GEOGRAPHICLIB_GEODESIC_ORDER
37 #define nC1p GEOGRAPHICLIB_GEODESIC_ORDER
38 #define nA2 GEOGRAPHICLIB_GEODESIC_ORDER
39 #define nC2 GEOGRAPHICLIB_GEODESIC_ORDER
40 #define nA3 GEOGRAPHICLIB_GEODESIC_ORDER
42 #define nC3 GEOGRAPHICLIB_GEODESIC_ORDER
43 #define nC3x ((nC3 * (nC3 - 1)) / 2)
44 #define nC4 GEOGRAPHICLIB_GEODESIC_ORDER
45 #define nC4x ((nC4 * (nC4 + 1)) / 2)
46 #define nC (GEOGRAPHICLIB_GEODESIC_ORDER + 1)
49 enum booly { FALSE = 0, TRUE = 1 };
53 enum dms { qd = 90, hd = 2 * qd, td = 2 * hd };
55 static unsigned init = 0;
56 static unsigned digits, maxit1, maxit2;
57 static double epsilon, realmin, pi, degree, NaN,
58 tiny, tol0, tol1, tol2, tolb, xthresh;
60 static void Init(
void) {
62 digits = DBL_MANT_DIG;
63 epsilon = DBL_EPSILON;
68 pi = atan2(0.0, -1.0);
71 maxit2 = maxit1 + digits + 10;
81 xthresh = 1000 * tol2;
99 static double sq(
double x) {
return x * x; }
101 static double sumx(
double u,
double v,
double* t) {
102 volatile double s = u + v;
103 volatile double up = s - v;
104 volatile double vpp = s - up;
107 if (t) *t = s != 0 ? 0 - (up + vpp) : s;
114 static double polyvalx(
int N,
const double p[],
double x) {
115 double y = N < 0 ? 0 : *p++;
116 while (--N >= 0) y = y * x + *p++;
120 static void swapx(
double* x,
double* y)
121 {
double t = *x; *x = *y; *y = t; }
123 static void norm2(
double* sinx,
double* cosx) {
124 #if defined(_MSC_VER) && defined(_M_IX86)
133 double r = sqrt(*sinx * *sinx + *cosx * *cosx);
135 double r = hypot(*sinx, *cosx);
141 static double AngNormalize(
double x) {
142 double y = remainder(x, (
double)td);
143 return fabs(y) == hd ? copysign((
double)hd, x) : y;
146 static double LatFix(
double x)
147 {
return fabs(x) > qd ? NaN : x; }
149 static double AngDiff(
double x,
double y,
double* e) {
152 double t, d = sumx(remainder(-x, (
double)td), remainder( y, (
double)td), &t);
155 d = sumx(remainder(d, (
double)td), t, &t);
157 if (d == 0 || fabs(d) == hd)
160 d = copysign(d, t == 0 ? y - x : -t);
165 static double AngRound(
double x) {
167 const double z = 1.0/16.0;
168 volatile double y = fabs(x);
169 volatile double w = z - y;
171 y = w > 0 ? z - w : y;
172 return copysign(y, x);
175 static void sincosdx(
double x,
double* sinx,
double* cosx) {
178 double r, s, c;
int q = 0;
179 r = remquo(x, (
double)qd, &q);
183 s = sin(r); c = cos(r);
184 switch ((
unsigned)q & 3U) {
185 case 0U: *sinx = s; *cosx = c;
break;
186 case 1U: *sinx = c; *cosx = -s;
break;
187 case 2U: *sinx = -s; *cosx = -c;
break;
188 default: *sinx = -c; *cosx = s;
break;
193 if (*sinx == 0) *sinx = copysign(*sinx, x);
196 static void sincosde(
double x,
double t,
double* sinx,
double* cosx) {
199 double r, s, c;
int q = 0;
200 r = AngRound(remquo(x, (
double)qd, &q) + t);
204 s = sin(r); c = cos(r);
205 switch ((
unsigned)q & 3U) {
206 case 0U: *sinx = s; *cosx = c;
break;
207 case 1U: *sinx = c; *cosx = -s;
break;
208 case 2U: *sinx = -s; *cosx = -c;
break;
209 default: *sinx = -c; *cosx = s;
break;
214 if (*sinx == 0) *sinx = copysign(*sinx, x);
217 static double atan2dx(
double y,
double x) {
222 int q = 0;
double ang;
223 if (fabs(y) > fabs(x)) { swapx(&x, &y); q = 2; }
224 if (signbit(x)) { x = -x; ++q; }
226 ang = atan2(y, x) / degree;
228 case 1: ang = copysign((
double)hd, y) - ang;
break;
229 case 2: ang = qd - ang;
break;
230 case 3: ang = -qd + ang;
break;
239 static double SinCosSeries(boolx sinp,
240 double sinx,
double cosx,
241 const double c[],
int n);
243 double eps,
double sig12,
244 double ssig1,
double csig1,
double dn1,
245 double ssig2,
double csig2,
double dn2,
246 double cbet1,
double cbet2,
247 double* ps12b,
double* pm12b,
double* pm0,
248 double* pM12,
double* pM21,
251 static double Astroid(
double x,
double y);
253 double sbet1,
double cbet1,
double dn1,
254 double sbet2,
double cbet2,
double dn2,
255 double lam12,
double slam12,
double clam12,
256 double* psalp1,
double* pcalp1,
258 double* psalp2,
double* pcalp2,
264 double sbet1,
double cbet1,
double dn1,
265 double sbet2,
double cbet2,
double dn2,
266 double salp1,
double calp1,
267 double slam120,
double clam120,
268 double* psalp2,
double* pcalp2,
270 double* pssig1,
double* pcsig1,
271 double* pssig2,
double* pcsig2,
274 boolx diffp,
double* pdlam12,
277 static double A3f(
const struct geod_geodesic* g,
double eps);
278 static void C3f(
const struct geod_geodesic* g,
double eps,
double c[]);
279 static void C4f(
const struct geod_geodesic* g,
double eps,
double c[]);
280 static double A1m1f(
double eps);
281 static void C1f(
double eps,
double c[]);
282 static void C1pf(
double eps,
double c[]);
283 static double A2m1f(
double eps);
284 static void C2f(
double eps,
double c[]);
285 static int transit(
double lon1,
double lon2);
286 static int transitdirect(
double lon1,
double lon2);
287 static void accini(
double s[]);
288 static void acccopy(
const double s[],
double t[]);
289 static void accadd(
double s[],
double y);
290 static double accsum(
const double s[],
double y);
291 static void accneg(
double s[]);
292 static void accrem(
double s[],
double y);
293 static double areareduceA(
double area[],
double area0,
294 int crossings, boolx reverse, boolx sign);
295 static double areareduceB(
double area,
double area0,
296 int crossings, boolx reverse, boolx sign);
303 g->e2 = g->
f * (2 - g->
f);
304 g->ep2 = g->e2 / sq(g->f1);
305 g->n = g->
f / ( 2 - g->
f);
307 g->c2 = (sq(g->
a) + sq(g->b) *
309 (g->e2 > 0 ? atanh(sqrt(g->e2)) : atan(sqrt(-g->e2))) /
310 sqrt(fabs(g->e2))))/2;
320 g->etol2 = 0.1 * tol2 /
321 sqrt( fmax(0.001, fabs(g->
f)) * fmin(1.0, 1 - g->
f/2) / 2 );
330 double lat1,
double lon1,
331 double azi1,
double salp1,
double calp1,
333 double cbet1, sbet1, eps;
344 l->
lat1 = LatFix(lat1);
350 sincosdx(AngRound(l->
lat1), &sbet1, &cbet1); sbet1 *= l->f1;
352 norm2(&sbet1, &cbet1); cbet1 = fmax(tiny, cbet1);
353 l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));
356 l->salp0 = l->
salp1 * cbet1;
369 l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;
370 l->csig1 = l->comg1 = sbet1 != 0 || l->
calp1 != 0 ? cbet1 * l->
calp1 : 1;
371 norm2(&l->ssig1, &l->csig1);
374 l->k2 = sq(l->calp0) * g->ep2;
375 eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);
377 if (l->
caps & CAP_C1) {
379 l->A1m1 = A1m1f(eps);
381 l->B11 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C1a, nC1);
382 s = sin(l->B11); c = cos(l->B11);
384 l->stau1 = l->ssig1 * c + l->csig1 * s;
385 l->ctau1 = l->csig1 * c - l->ssig1 * s;
390 if (l->
caps & CAP_C1p)
393 if (l->
caps & CAP_C2) {
394 l->A2m1 = A2m1f(eps);
396 l->B21 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C2a, nC2);
399 if (l->
caps & CAP_C3) {
401 l->A3c = -l->
f * l->salp0 * A3f(g, eps);
402 l->B31 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C3a, nC3-1);
405 if (l->
caps & CAP_C4) {
408 l->A4 = sq(l->
a) * l->calp0 * l->salp0 * g->e2;
409 l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4);
417 double lat1,
double lon1,
double azi1,
unsigned caps) {
419 azi1 = AngNormalize(azi1);
421 sincosdx(AngRound(azi1), &salp1, &calp1);
422 geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
427 double lat1,
double lon1,
double azi1,
428 unsigned flags,
double s12_a12,
436 double lat1,
double lon1,
double azi1,
437 double s12,
unsigned caps) {
442 unsigned flags,
double s12_a12,
443 double* plat2,
double* plon2,
double* pazi2,
444 double* ps12,
double* pm12,
445 double* pM12,
double* pM21,
447 double lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,
448 m12 = 0, M12 = 0, M21 = 0, S12 = 0;
450 double sig12, ssig12, csig12, B12 = 0, AB1 = 0;
451 double omg12, lam12, lon12;
452 double ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
462 outmask &= l->
caps & OUT_ALL;
469 sig12 = s12_a12 * degree;
470 sincosdx(s12_a12, &ssig12, &csig12);
474 tau12 = s12_a12 / (l->b * (1 + l->A1m1)),
478 B12 = - SinCosSeries(TRUE,
479 l->stau1 * c + l->ctau1 * s,
480 l->ctau1 * c - l->stau1 * s,
482 sig12 = tau12 - (B12 - l->B11);
483 ssig12 = sin(sig12); csig12 = cos(sig12);
484 if (fabs(l->
f) > 0.01) {
507 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
508 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
509 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
510 serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b;
511 sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2));
512 ssig12 = sin(sig12); csig12 = cos(sig12);
518 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
519 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
520 dn2 = sqrt(1 + l->k2 * sq(ssig2));
523 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
524 AB1 = (1 + l->A1m1) * (B12 - l->B11);
527 sbet2 = l->calp0 * ssig2;
529 cbet2 = hypot(l->salp0, l->calp0 * csig2);
532 cbet2 = csig2 = tiny;
534 salp2 = l->salp0; calp2 = l->calp0 * csig2;
538 l->b * ((1 + l->A1m1) * sig12 + AB1) :
542 double E = copysign(1, l->salp0);
544 somg2 = l->salp0 * ssig2; comg2 = csig2;
548 - (atan2( ssig2, csig2) - atan2( l->ssig1, l->csig1))
549 + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))
550 : atan2(somg2 * l->comg1 - comg2 * l->somg1,
551 comg2 * l->comg1 + somg2 * l->somg1);
552 lam12 = omg12 + l->A3c *
553 ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
555 lon12 = lam12 / degree;
557 AngNormalize(AngNormalize(l->
lon1) + AngNormalize(lon12));
561 lat2 = atan2dx(sbet2, l->f1 * cbet2);
564 azi2 = atan2dx(salp2, calp2);
568 B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),
569 AB2 = (1 + l->A2m1) * (B22 - l->B21),
570 J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);
574 m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
575 - l->csig1 * csig2 * J12);
577 double t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) /
579 M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1;
580 M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2;
586 B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
587 double salp12, calp12;
588 if (l->calp0 == 0 || l->salp0 == 0) {
601 salp12 = l->calp0 * l->salp0 *
602 (csig12 <= 0 ? l->csig1 * (1 - csig12) + ssig12 * l->ssig1 :
603 ssig12 * (l->csig1 * ssig12 / (1 + csig12) + l->ssig1));
604 calp12 = sq(l->salp0) + sq(l->calp0) * l->csig1 * csig2;
606 S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
628 if (pM12) *pM12 = M12;
629 if (pM21) *pM21 = M21;
634 return (flags &
GEOD_ARCMODE) ? s12_a12 : sig12 / degree;
640 nullptr,
nullptr,
nullptr,
nullptr,
nullptr);
644 l->
a13 = a13; l->
s13 = NaN;
646 nullptr,
nullptr,
nullptr,
nullptr);
650 unsigned flags,
double s13_a13) {
652 geod_setarc(l, s13_a13) :
657 double* plat2,
double* plon2,
double* pazi2) {
659 nullptr,
nullptr,
nullptr,
nullptr,
nullptr);
663 double lat1,
double lon1,
double azi1,
664 unsigned flags,
double s12_a12,
665 double* plat2,
double* plon2,
double* pazi2,
666 double* ps12,
double* pm12,
double* pM12,
double* pM21,
683 plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
689 double* plat2,
double* plon2,
double* pazi2) {
691 nullptr,
nullptr,
nullptr,
nullptr,
nullptr);
694 static double geod_geninverse_int(
const struct geod_geodesic* g,
696 double lat2,
double lon2,
698 double* psalp1,
double* pcalp1,
699 double* psalp2,
double* pcalp2,
700 double* pm12,
double* pM12,
double* pM21,
702 double s12 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
703 double lon12, lon12s;
704 int latsign, lonsign, swapp;
705 double sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0;
706 double dn1, dn2, lam12, slam12, clam12;
707 double a12 = 0, sig12,
calp1 = 0,
salp1 = 0, calp2 = 0, salp2 = 0;
711 double omg12 = 0, somg12 = 2, comg12 = 0;
723 lon12 = AngDiff(
lon1, lon2, &lon12s);
725 lonsign = signbit(lon12) ? -1 : 1;
726 lon12 *= lonsign; lon12s *= lonsign;
727 lam12 = lon12 * degree;
729 sincosde(lon12, lon12s, &slam12, &clam12);
730 lon12s = (hd - lon12) - lon12s;
734 lat2 = AngRound(LatFix(lat2));
737 swapp = fabs(
lat1) < fabs(lat2) || lat2 != lat2 ? -1 : 1;
743 latsign = signbit(
lat1) ? 1 : -1;
758 sincosdx(
lat1, &sbet1, &cbet1); sbet1 *= g->f1;
760 norm2(&sbet1, &cbet1); cbet1 = fmax(tiny, cbet1);
762 sincosdx(lat2, &sbet2, &cbet2); sbet2 *= g->f1;
764 norm2(&sbet2, &cbet2); cbet2 = fmax(tiny, cbet2);
774 if (cbet1 < -sbet1) {
776 sbet2 = copysign(sbet1, sbet2);
778 if (fabs(sbet2) == -sbet1)
782 dn1 = sqrt(1 + g->ep2 * sq(sbet1));
783 dn2 = sqrt(1 + g->ep2 * sq(sbet2));
785 meridian =
lat1 == -qd || slam12 == 0;
792 double ssig1, csig1, ssig2, csig2;
794 calp2 = 1; salp2 = 0;
797 ssig1 = sbet1; csig1 =
calp1 * cbet1;
798 ssig2 = sbet2; csig2 = calp2 * cbet2;
801 sig12 = atan2(fmax(0.0, csig1 * ssig2 - ssig1 * csig2) + 0,
802 csig1 * csig2 + ssig1 * ssig2);
803 Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
804 cbet1, cbet2, &s12x, &m12x,
nullptr,
815 if (sig12 < 1 || m12x >= 0) {
817 if (sig12 < 3 * tiny ||
819 (sig12 < tol0 && (s12x < 0 || m12x < 0)))
820 sig12 = m12x = s12x = 0;
823 a12 = sig12 / degree;
832 (g->
f <= 0 || lon12s >= g->
f * hd)) {
837 sig12 = omg12 = lam12 / g->f1;
838 m12x = g->b * sin(sig12);
840 M12 = M21 = cos(sig12);
843 }
else if (!meridian) {
850 sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
851 lam12, slam12, clam12,
857 s12x = sig12 * g->b * dnm;
858 m12x = sq(dnm) * g->b * sin(sig12 / dnm);
860 M12 = M21 = cos(sig12 / dnm);
861 a12 = sig12 / degree;
862 omg12 = lam12 / (g->f1 * dnm);
876 double ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0;
879 double salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
886 v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
salp1,
calp1,
888 &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
889 &eps, &domg12, numit < maxit1, &dv, Ca);
892 !(fabs(v) >= (tripn ? 8 : 1) * tol0) ||
897 if (v > 0 && (numit > maxit1 ||
calp1/
salp1 > calp1b/salp1b))
899 else if (v < 0 && (numit > maxit1 ||
calp1/
salp1 < calp1a/salp1a))
901 if (numit < maxit1 && dv > 0) {
904 if (fabs(dalp1) < pi) {
906 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
915 tripn = fabs(v) <= 16 * tol0;
928 salp1 = (salp1a + salp1b)/2;
929 calp1 = (calp1a + calp1b)/2;
932 tripb = (fabs(salp1a -
salp1) + (calp1a -
calp1) < tolb ||
933 fabs(
salp1 - salp1b) + (
calp1 - calp1b) < tolb);
935 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
936 cbet1, cbet2, &s12x, &m12x,
nullptr,
941 a12 = sig12 / degree;
944 double sdomg12 = sin(domg12), cdomg12 = cos(domg12);
945 somg12 = slam12 * cdomg12 - clam12 * sdomg12;
946 comg12 = clam12 * cdomg12 + slam12 * sdomg12;
960 salp0 =
salp1 * cbet1,
963 if (calp0 != 0 && salp0 != 0) {
966 ssig1 = sbet1, csig1 =
calp1 * cbet1,
967 ssig2 = sbet2, csig2 = calp2 * cbet2,
968 k2 = sq(calp0) * g->ep2,
969 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
971 A4 = sq(g->
a) * calp0 * salp0 * g->e2;
973 norm2(&ssig1, &csig1);
974 norm2(&ssig2, &csig2);
976 B41 = SinCosSeries(FALSE, ssig1, csig1, Ca, nC4);
977 B42 = SinCosSeries(FALSE, ssig2, csig2, Ca, nC4);
978 S12 = A4 * (B42 - B41);
983 if (!meridian && somg12 == 2) {
984 somg12 = sin(omg12); comg12 = cos(omg12);
990 sbet2 - sbet1 < 1.75) {
995 domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
996 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
997 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
1007 if (salp12 == 0 && calp12 < 0) {
1008 salp12 = tiny *
calp1;
1011 alp12 = atan2(salp12, calp12);
1013 S12 += g->c2 * alp12;
1014 S12 *= swapp * lonsign * latsign;
1021 swapx(&
salp1, &salp2);
1022 swapx(&
calp1, &calp2);
1027 salp1 *= swapp * lonsign;
calp1 *= swapp * latsign;
1028 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
1030 if (psalp1) *psalp1 =
salp1;
1031 if (pcalp1) *pcalp1 =
calp1;
1032 if (psalp2) *psalp2 = salp2;
1033 if (pcalp2) *pcalp2 = calp2;
1040 if (pM12) *pM12 = M12;
1041 if (pM21) *pM21 = M21;
1051 double lat1,
double lon1,
double lat2,
double lon2,
1052 double* ps12,
double* pazi1,
double* pazi2,
1053 double* pm12,
double* pM12,
double* pM21,
1056 a12 = geod_geninverse_int(g,
lat1,
lon1, lat2, lon2, ps12,
1058 pm12, pM12, pM21, pS12);
1060 if (pazi2) *pazi2 = atan2dx(salp2, calp2);
1066 double lat1,
double lon1,
double lat2,
double lon2,
1069 a12 = geod_geninverse_int(g,
lat1,
lon1, lat2, lon2,
nullptr,
1071 nullptr,
nullptr,
nullptr,
nullptr),
1077 geod_setarc(l, a12);
1081 double lat1,
double lon1,
double lat2,
double lon2,
1082 double* ps12,
double* pazi1,
double* pazi2) {
1084 nullptr,
nullptr,
nullptr,
nullptr);
1087 double SinCosSeries(boolx sinp,
double sinx,
double cosx,
1088 const double c[],
int n) {
1096 ar = 2 * (cosx - sinx) * (cosx + sinx);
1097 y0 = (n & 1) ? *--c : 0; y1 = 0;
1102 y1 = ar * y0 - y1 + *--c;
1103 y0 = ar * y1 - y0 + *--c;
1106 ? 2 * sinx * cosx * y0
1111 double eps,
double sig12,
1112 double ssig1,
double csig1,
double dn1,
1113 double ssig2,
double csig2,
double dn2,
1114 double cbet1,
double cbet2,
1115 double* ps12b,
double* pm12b,
double* pm0,
1116 double* pM12,
double* pM21,
1119 double m0 = 0, J12 = 0, A1 = 0, A2 = 0;
1124 boolx redlp = pm12b || pm0 || pM12 || pM21;
1125 if (ps12b || redlp) {
1137 double B1 = SinCosSeries(TRUE, ssig2, csig2, Ca, nC1) -
1138 SinCosSeries(TRUE, ssig1, csig1, Ca, nC1);
1140 *ps12b = A1 * (sig12 + B1);
1142 double B2 = SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1143 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2);
1144 J12 = m0 * sig12 + (A1 * B1 - A2 * B2);
1149 for (l = 1; l <= nC2; ++l)
1150 Cb[l] = A1 * Ca[l] - A2 * Cb[l];
1151 J12 = m0 * sig12 + (SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1152 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2));
1159 *pm12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1160 csig1 * csig2 * J12;
1162 double csig12 = csig1 * csig2 + ssig1 * ssig2;
1163 double t = g->ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
1165 *pM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
1167 *pM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
1171 double Astroid(
double x,
double y) {
1178 r = (p + q - 1) / 6;
1179 if ( !(q == 0 && r <= 0) ) {
1188 disc = S * (S + 2 * r3);
1192 double T3 = S + r3, T;
1196 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
1200 u += T + (T != 0 ? r2 / T : 0);
1203 double ang = atan2(sqrt(-disc), -(S + r3));
1206 u += 2 * r * cos(ang / 3);
1208 v = sqrt(sq(u) + q);
1210 uv = u < 0 ? q / (v - u) : u + v;
1211 w = (uv - q) / (2 * v);
1214 k = uv / (sqrt(uv + sq(w)) + w);
1224 double sbet1,
double cbet1,
double dn1,
1225 double sbet2,
double cbet2,
double dn2,
1226 double lam12,
double slam12,
double clam12,
1227 double* psalp1,
double* pcalp1,
1229 double* psalp2,
double* pcalp2,
1234 double salp1 = 0,
calp1 = 0, salp2 = 0, calp2 = 0, dnm = 0;
1242 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
1243 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
1245 boolx shortline = cbet12 >= 0 && sbet12 < 0.5 && cbet2 * lam12 < 0.5;
1246 double somg12, comg12, ssig12, csig12;
1247 sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
1249 double sbetm2 = sq(sbet1 + sbet2), omg12;
1252 sbetm2 /= sbetm2 + sq(cbet1 + cbet2);
1253 dnm = sqrt(1 + g->ep2 * sbetm2);
1254 omg12 = lam12 / (g->f1 * dnm);
1255 somg12 = sin(omg12); comg12 = cos(omg12);
1257 somg12 = slam12; comg12 = clam12;
1260 salp1 = cbet2 * somg12;
1261 calp1 = comg12 >= 0 ?
1262 sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) :
1263 sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1266 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
1268 if (shortline && ssig12 < g->etol2) {
1270 salp2 = cbet1 * somg12;
1271 calp2 = sbet12 - cbet1 * sbet2 *
1272 (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);
1273 norm2(&salp2, &calp2);
1275 sig12 = atan2(ssig12, csig12);
1276 }
else if (fabs(g->n) > 0.1 ||
1278 ssig12 >= 6 * fabs(g->n) * pi * sq(cbet1)) {
1283 double x, y, lamscale, betscale;
1284 double lam12x = atan2(-slam12, -clam12);
1289 k2 = sq(sbet1) * g->ep2,
1290 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1291 lamscale = g->
f * cbet1 * A3f(g, eps) * pi;
1293 betscale = lamscale * cbet1;
1295 x = lam12x / lamscale;
1296 y = sbet12a / betscale;
1300 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1301 bet12a = atan2(sbet12a, cbet12a);
1305 Lengths(g, g->n, pi + bet12a,
1306 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1307 cbet1, cbet2,
nullptr, &m12b, &m0,
nullptr,
nullptr, Ca);
1308 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi);
1309 betscale = x < -0.01 ? sbet12a / x :
1310 -g->
f * sq(cbet1) * pi;
1311 lamscale = betscale / cbet1;
1312 y = lam12x / lamscale;
1315 if (y > -tol1 && x > -1 - xthresh) {
1320 calp1 = fmax(x > -tol1 ? 0.0 : -1.0, x);
1358 double k = Astroid(x, y);
1360 omg12a = lamscale * ( g->
f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
1361 somg12 = sin(omg12a); comg12 = -cos(omg12a);
1363 salp1 = cbet2 * somg12;
1364 calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1386 double sbet1,
double cbet1,
double dn1,
1387 double sbet2,
double cbet2,
double dn2,
1389 double slam120,
double clam120,
1390 double* psalp2,
double* pcalp2,
1392 double* pssig1,
double* pcsig1,
1393 double* pssig2,
double* pcsig2,
1396 boolx diffp,
double* pdlam12,
1399 double salp2 = 0, calp2 = 0, sig12 = 0,
1400 ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0,
1401 domg12 = 0, dlam12 = 0;
1402 double salp0, calp0;
1403 double somg1, comg1, somg2, comg2, somg12, comg12, lam12;
1404 double B312, eta, k2;
1406 if (sbet1 == 0 &&
calp1 == 0)
1412 salp0 =
salp1 * cbet1;
1417 ssig1 = sbet1; somg1 = salp0 * sbet1;
1418 csig1 = comg1 =
calp1 * cbet1;
1419 norm2(&ssig1, &csig1);
1426 salp2 = cbet2 != cbet1 ? salp0 / cbet2 :
salp1;
1431 calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ?
1432 sqrt(sq(
calp1 * cbet1) +
1434 (cbet2 - cbet1) * (cbet1 + cbet2) :
1435 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1439 ssig2 = sbet2; somg2 = salp0 * sbet2;
1440 csig2 = comg2 = calp2 * cbet2;
1441 norm2(&ssig2, &csig2);
1445 sig12 = atan2(fmax(0.0, csig1 * ssig2 - ssig1 * csig2) + 0,
1446 csig1 * csig2 + ssig1 * ssig2);
1449 somg12 = fmax(0.0, comg1 * somg2 - somg1 * comg2) + 0;
1450 comg12 = comg1 * comg2 + somg1 * somg2;
1452 eta = atan2(somg12 * clam120 - comg12 * slam120,
1453 comg12 * clam120 + somg12 * slam120);
1454 k2 = sq(calp0) * g->ep2;
1455 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1457 B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) -
1458 SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1));
1459 domg12 = -g->
f * A3f(g, eps) * salp0 * (sig12 + B312);
1460 lam12 = eta + domg12;
1464 dlam12 = - 2 * g->f1 * dn1 / sbet1;
1466 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1467 cbet1, cbet2,
nullptr, &dlam12,
nullptr,
nullptr,
nullptr, Ca);
1468 dlam12 *= g->f1 / (calp2 * cbet2);
1489 return polyvalx(nA3 - 1, g->A3x, eps);
1492 void C3f(
const struct geod_geodesic* g,
double eps,
double c[]) {
1497 for (l = 1; l < nC3; ++l) {
1498 int m = nC3 - l - 1;
1500 c[l] = mult * polyvalx(m, g->C3x + o, eps);
1505 void C4f(
const struct geod_geodesic* g,
double eps,
double c[]) {
1510 for (l = 0; l < nC4; ++l) {
1511 int m = nC4 - l - 1;
1512 c[l] = mult * polyvalx(m, g->C4x + o, eps);
1519 double A1m1f(
double eps) {
1520 static const double coeff[] = {
1525 double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1];
1526 return (t + eps) / (1 - eps);
1530 void C1f(
double eps,
double c[]) {
1531 static const double coeff[] = {
1549 for (l = 1; l <= nC1; ++l) {
1550 int m = (nC1 - l) / 2;
1551 c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];
1558 void C1pf(
double eps,
double c[]) {
1559 static const double coeff[] = {
1561 205, -432, 768, 1536,
1563 4005, -4736, 3840, 12288,
1577 for (l = 1; l <= nC1p; ++l) {
1578 int m = (nC1p - l) / 2;
1579 c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];
1586 double A2m1f(
double eps) {
1587 static const double coeff[] = {
1589 -11, -28, -192, 0, 256,
1592 double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1];
1593 return (t - eps) / (1 + eps);
1597 void C2f(
double eps,
double c[]) {
1598 static const double coeff[] = {
1616 for (l = 1; l <= nC2; ++l) {
1617 int m = (nC2 - l) / 2;
1618 c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];
1626 static const double coeff[] = {
1640 int o = 0, k = 0, j;
1641 for (j = nA3 - 1; j >= 0; --j) {
1642 int m = nA3 - j - 1 < j ? nA3 - j - 1 : j;
1643 g->A3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];
1650 static const double coeff[] = {
1682 int o = 0, k = 0, l, j;
1683 for (l = 1; l < nC3; ++l) {
1684 for (j = nC3 - 1; j >= l; --j) {
1685 int m = nC3 - j - 1 < j ? nC3 - j - 1 : j;
1686 g->C3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];
1694 static const double coeff[] = {
1700 -224, -4784, 1573, 45045,
1702 -10656, 14144, -4576, -858, 45045,
1704 64, 624, -4576, 6864, -3003, 15015,
1706 100, 208, 572, 3432, -12012, 30030, 45045,
1712 5792, 1040, -1287, 135135,
1714 5952, -11648, 9152, -2574, 135135,
1716 -64, -624, 4576, -6864, 3003, 135135,
1722 -8448, 4992, -1144, 225225,
1724 -1440, 4160, -4576, 1716, 225225,
1730 3584, -3328, 1144, 315315,
1738 int o = 0, k = 0, l, j;
1739 for (l = 0; l < nC4; ++l) {
1740 for (j = nC4 - 1; j >= l; --j) {
1741 int m = nC4 - j - 1;
1742 g->C4x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];
1748 int transit(
double lon1,
double lon2) {
1753 lon12 = AngDiff(
lon1, lon2,
nullptr);
1755 lon2 = AngNormalize(lon2);
1757 lon12 > 0 && ((lon1 < 0 && lon2 >= 0) ||
1758 (
lon1 > 0 && lon2 == 0)) ? 1 :
1759 (lon12 < 0 && lon1 >= 0 && lon2 < 0 ? -1 : 0);
1762 int transitdirect(
double lon1,
double lon2) {
1765 lon1 = remainder(
lon1, 2.0 * td); lon2 = remainder(lon2, 2.0 * td);
1766 return ( (lon2 >= 0 && lon2 < td ? 0 : 1) -
1767 (
lon1 >= 0 &&
lon1 < td ? 0 : 1) );
1770 void accini(
double s[]) {
1775 void acccopy(
const double s[],
double t[]) {
1777 t[0] = s[0]; t[1] = s[1];
1780 void accadd(
double s[],
double y) {
1782 double u, z = sumx(y, s[1], &u);
1783 s[0] = sumx(z, s[0], &s[1]);
1790 double accsum(
const double s[],
double y) {
1798 void accneg(
double s[]) {
1800 s[0] = -s[0]; s[1] = -s[1];
1803 void accrem(
double s[],
double y) {
1805 s[0] = remainder(s[0], y);
1810 p->polyline = (polylinep != 0);
1815 p->lat0 = p->lon0 = p->
lat = p->
lon = NaN;
1818 p->
num = p->crossings = 0;
1823 double lat,
double lon) {
1825 p->lat0 = p->
lat = lat;
1826 p->lon0 = p->
lon = lon;
1828 double s12, S12 = 0;
1830 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
1831 p->polyline ?
nullptr : &S12);
1835 p->crossings += transit(p->
lon, lon);
1837 p->
lat = lat; p->
lon = lon;
1844 double azi,
double s) {
1848 double lat = 0, lon = 0, S12 = 0;
1850 &lat, &lon,
nullptr,
1851 nullptr,
nullptr,
nullptr,
nullptr,
1852 p->polyline ?
nullptr : &S12);
1856 p->crossings += transitdirect(p->
lon, lon);
1858 p->
lat = lat; p->
lon = lon;
1865 boolx reverse, boolx sign,
1866 double* pA,
double* pP) {
1867 double s12, S12, t[2];
1870 if (!p->polyline && pA) *pA = 0;
1874 if (pP) *pP = p->P[0];
1878 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr, &S12);
1879 if (pP) *pP = accsum(p->P, s12);
1882 if (pA) *pA = areareduceA(t, 4 * pi * g->c2,
1883 p->crossings + transit(p->
lon, p->lon0),
1890 double lat,
double lon,
1891 boolx reverse, boolx sign,
1892 double* pA,
double* pP) {
1893 double perimeter, tempsum;
1895 unsigned num = p->
num + 1;
1898 if (!p->polyline && pA) *pA = 0;
1901 perimeter = p->P[0];
1902 tempsum = p->polyline ? 0 : p->A[0];
1903 crossings = p->crossings;
1904 for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
1905 double s12, S12 = 0;
1907 i == 0 ? p->
lat : lat, i == 0 ? p->
lon : lon,
1908 i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
1909 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
1910 p->polyline ?
nullptr : &S12);
1914 crossings += transit(i == 0 ? p->
lon : lon,
1915 i != 0 ? p->lon0 : lon);
1919 if (pP) *pP = perimeter;
1923 if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
1929 double azi,
double s,
1930 boolx reverse, boolx sign,
1931 double* pA,
double* pP) {
1932 double perimeter, tempsum;
1934 unsigned num = p->
num + 1;
1937 if (!p->polyline && pA) *pA = NaN;
1940 perimeter = p->P[0] + s;
1942 if (pP) *pP = perimeter;
1947 crossings = p->crossings;
1951 double lat = 0, lon = 0, s12, S12 = 0;
1953 &lat, &lon,
nullptr,
1954 nullptr,
nullptr,
nullptr,
nullptr, &S12);
1956 crossings += transitdirect(p->
lon, lon);
1958 &s12,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr, &S12);
1961 crossings += transit(lon, p->lon0);
1964 if (pP) *pP = perimeter;
1965 if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
1970 double lats[],
double lons[],
int n,
1971 double* pA,
double* pP) {
1975 for (i = 0; i < n; ++i)
1980 double areareduceA(
double area[],
double area0,
1981 int crossings, boolx reverse, boolx sign) {
1982 accrem(area, area0);
1984 accadd(area, (area[0] < 0 ? 1 : -1) * area0/2);
1991 if (area[0] > area0/2)
1992 accadd(area, -area0);
1993 else if (area[0] <= -area0/2)
1994 accadd(area, +area0);
1996 if (area[0] >= area0)
1997 accadd(area, -area0);
1998 else if (area[0] < 0)
1999 accadd(area, +area0);
2004 double areareduceB(
double area,
double area0,
2005 int crossings, boolx reverse, boolx sign) {
2006 area = remainder(area, area0);
2008 area += (area < 0 ? 1 : -1) * area0/2;
2017 else if (area <= -area0/2)
API for the geodesic routines in C.
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)
void GEOD_DLL geod_polygon_addedge(const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)
void GEOD_DLL geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
void GEOD_DLL geod_setdistance(struct geod_geodesicline *l, double s13)
void GEOD_DLL geod_gensetdistance(struct geod_geodesicline *l, unsigned flags, double s13_a13)
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)
void GEOD_DLL geod_polygon_clear(struct geod_polygon *p)
void GEOD_DLL geod_init(struct geod_geodesic *g, double a, double f)
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_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
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_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
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_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
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)
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)
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_polygonarea(const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)
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_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_init(struct geod_polygon *p, int polylinep)