35 # pragma warning (disable: 4701 4127) 43 : maxit2_(maxit1_ +
Math::digits() + 10)
47 , tiny_(sqrt(numeric_limits<real>::min()))
48 , tol0_(numeric_limits<real>::epsilon())
55 , tolb_(tol0_ * tol2_)
56 , xthresh_(1000 * tol2_)
58 , _f(f <= 1 ? f : 1/f)
61 , _ep2(_e2 /
Math::sq(_f1))
65 Math::eatanhe(real(1), (_f < 0 ? -1 : 1) * sqrt(abs(_e2))) / _e2)
77 , _etol2(0.1 * tol2_ /
78 sqrt( max(real(0.001), abs(_f)) * min(real(1), 1 - _f/2) / 2 ))
96 const real c[],
int n) {
104 ar = 2 * (cosx - sinx) * (cosx + sinx),
105 y0 = n & 1 ? *--c : 0, y1 = 0;
110 y1 = ar * y0 - y1 + *--c;
111 y0 = ar * y1 - y0 + *--c;
114 ? 2 * sinx * cosx * y0
124 bool arcmode, real s12_a12,
unsigned outmask,
125 real& lat2, real& lon2, real& azi2,
126 real& s12, real& m12, real& M12, real& M21,
132 GenPosition(arcmode, s12_a12, outmask,
133 lat2, lon2, azi2, s12, m12, M12, M21, S12);
138 real& s12, real& azi1, real& azi2,
139 real& m12, real& M12, real& M21, real& S12)
150 int lonsign = lon12 >= 0 ? 1 : -1;
156 int swapp = abs(lat1) >= abs(lat2) ? 1 : -1;
162 int latsign = lat1 < 0 ? 1 : -1;
177 real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x;
181 sbet1 = _f1 * sin(phi);
182 cbet1 = lat1 == -90 ? tiny_ : cos(phi);
187 sbet2 = _f1 * sin(phi);
188 cbet2 = abs(lat2) == 90 ? tiny_ : cos(phi);
199 if (cbet1 < -sbet1) {
201 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
203 if (abs(sbet2) == -sbet1)
208 dn1 = sqrt(1 + _ep2 *
Math::sq(sbet1)),
209 dn2 = sqrt(1 + _ep2 *
Math::sq(sbet2));
213 slam12 = abs(lon12) == 180 ? 0 : sin(lam12),
217 real a12, sig12, calp1, salp1, calp2 = 0, salp2 = 0;
219 real C1a[nC1_ + 1], C2a[nC2_ + 1], C3a[nC3_];
221 bool meridian = lat1 == -90 || slam12 == 0;
228 calp1 = clam12; salp1 = slam12;
229 calp2 = 1; salp2 = 0;
233 ssig1 = sbet1, csig1 = calp1 * cbet1,
234 ssig2 = sbet2, csig2 = calp2 * cbet2;
237 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
238 csig1 * csig2 + ssig1 * ssig2);
241 Lengths(_n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
242 cbet1, cbet2, s12x, m12x, dummy,
252 if (sig12 < 1 || m12x >= 0) {
268 calp1 = calp2 = 0; salp1 = salp2 = 1;
270 sig12 = omg12 = lam12 / _f1;
271 m12x = _b * sin(sig12);
273 M12 = M21 = cos(sig12);
276 }
else if (!meridian) {
283 sig12 = InverseStart(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
285 salp1, calp1, salp2, calp2, dnm,
290 s12x = sig12 * _b * dnm;
291 m12x =
Math::sq(dnm) * _b * sin(sig12 / dnm);
293 M12 = M21 = cos(sig12 / dnm);
295 omg12 = lam12 / (_f1 * dnm);
311 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0;
314 real salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
315 for (
bool tripn =
false, tripb =
false;
321 real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
322 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
323 eps, omg12, numit < maxit1_, dv, C1a, C2a, C3a)
327 if (tripb || !(abs(v) >= (tripn ? 8 : 2) * tol0_))
break;
329 if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
330 { salp1b = salp1; calp1b = calp1; }
331 else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
332 { salp1a = salp1; calp1a = calp1; }
333 if (numit < maxit1_ && dv > 0) {
337 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
338 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
339 if (nsalp1 > 0 && abs(dalp1) <
Math::pi()) {
340 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
346 tripn = abs(v) <= 16 * tol0_;
358 salp1 = (salp1a + salp1b)/2;
359 calp1 = (calp1a + calp1b)/2;
362 tripb = (abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
363 abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
367 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
368 cbet1, cbet2, s12x, m12x, dummy,
374 omg12 = lam12 - omg12;
384 if (outmask &
AREA) {
387 salp0 = salp1 * cbet1,
390 if (calp0 != 0 && salp0 != 0) {
393 ssig1 = sbet1, csig1 = calp1 * cbet1,
394 ssig2 = sbet2, csig2 = calp2 * cbet2,
396 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
398 A4 =
Math::sq(_a) * calp0 * salp0 * _e2;
404 B41 = SinCosSeries(
false, ssig1, csig1, C4a, nC4_),
405 B42 = SinCosSeries(
false, ssig2, csig2, C4a, nC4_);
406 S12 = A4 * (B42 - B41);
413 sbet2 - sbet1 < real(1.75)) {
418 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
419 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
420 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
421 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
425 salp12 = salp2 * calp1 - calp2 * salp1,
426 calp12 = calp2 * calp1 + salp2 * salp1;
431 if (salp12 == 0 && calp12 < 0) {
432 salp12 = tiny_ * calp1;
435 alp12 = atan2(salp12, calp12);
438 S12 *= swapp * lonsign * latsign;
451 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
452 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
463 void Geodesic::Lengths(real eps, real sig12,
464 real ssig1, real csig1, real dn1,
465 real ssig2, real csig2, real dn2,
466 real cbet1, real cbet2,
467 real& s12b, real& m12b, real& m0,
468 bool scalep, real& M12, real& M21,
470 real C1a[], real C2a[])
const {
477 AB1 = (1 + A1m1) * (SinCosSeries(
true, ssig2, csig2, C1a, nC1_) -
478 SinCosSeries(
true, ssig1, csig1, C1a, nC1_)),
480 AB2 = (1 + A2m1) * (SinCosSeries(
true, ssig2, csig2, C2a, nC2_) -
481 SinCosSeries(
true, ssig1, csig1, C2a, nC2_));
483 real J12 = m0 * sig12 + (AB1 - AB2);
487 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12;
489 s12b = (1 + A1m1) * sig12 + AB1;
491 real csig12 = csig1 * csig2 + ssig1 * ssig2;
492 real t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
493 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
494 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
498 Math::real Geodesic::Astroid(real x, real y) {
506 if ( !(q == 0 && r <= 0) ) {
515 disc = S * (S + 2 * r3);
522 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
526 u += T + (T ? r2 / T : 0);
529 real ang = atan2(sqrt(-disc), -(S + r3));
532 u += 2 * r * cos(ang / 3);
537 uv = u < 0 ? q / (v - u) : u + v,
538 w = (uv - q) / (2 * v);
541 k = uv / (sqrt(uv +
Math::sq(w)) + w);
550 Math::real Geodesic::InverseStart(real sbet1, real cbet1, real dn1,
551 real sbet2, real cbet2, real dn2,
553 real& salp1, real& calp1,
555 real& salp2, real& calp2,
559 real C1a[], real C2a[])
const {
566 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
567 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
568 #if defined(__GNUC__) && __GNUC__ == 4 && \ 569 (__GNUC_MINOR__ < 6 || defined(__MINGW32__)) 583 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
585 bool shortline = cbet12 >= 0 && sbet12 < real(0.5) &&
586 cbet2 * lam12 < real(0.5);
589 real sbetm2 =
Math::sq(sbet1 + sbet2);
592 sbetm2 /= sbetm2 +
Math::sq(cbet1 + cbet2);
593 dnm = sqrt(1 + _ep2 * sbetm2);
596 real somg12 = sin(omg12), comg12 = cos(omg12);
598 salp1 = cbet2 * somg12;
599 calp1 = comg12 >= 0 ?
600 sbet12 + cbet2 * sbet1 *
Math::sq(somg12) / (1 + comg12) :
601 sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
605 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
607 if (shortline && ssig12 < _etol2) {
609 salp2 = cbet1 * somg12;
610 calp2 = sbet12 - cbet1 * sbet2 *
611 (comg12 >= 0 ?
Math::sq(somg12) / (1 + comg12) : 1 - comg12);
614 sig12 = atan2(ssig12, csig12);
615 }
else if (abs(_n) > real(0.1) ||
622 real y, lamscale, betscale;
632 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
633 lamscale = _f * cbet1 * A3f(eps) *
Math::pi();
635 betscale = lamscale * cbet1;
637 x = (lam12 -
Math::pi()) / lamscale;
638 y = sbet12a / betscale;
642 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
643 bet12a = atan2(sbet12a, cbet12a);
644 real m12b, m0, dummy;
648 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
649 cbet1, cbet2, dummy, m12b, m0,
false,
650 dummy, dummy, C1a, C2a);
651 x = -1 + m12b / (cbet1 * cbet2 * m0 *
Math::pi());
652 betscale = x < -real(0.01) ? sbet12a / x :
654 lamscale = betscale / cbet1;
655 y = (lam12 -
Math::pi()) / lamscale;
658 if (y > -tol1_ && x > -1 - xthresh_) {
662 salp1 = min(real(1), -real(x)); calp1 = - sqrt(1 -
Math::sq(salp1));
664 calp1 = max(real(x > -tol1_ ? 0 : -1), real(x));
702 real k = Astroid(x, y);
704 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
705 somg12 = sin(omg12a); comg12 = -cos(omg12a);
707 salp1 = cbet2 * somg12;
708 calp1 = sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
715 salp1 = 1; calp1 = 0;
720 Math::real Geodesic::Lambda12(real sbet1, real cbet1, real dn1,
721 real sbet2, real cbet2, real dn2,
722 real salp1, real calp1,
723 real& salp2, real& calp2,
725 real& ssig1, real& csig1,
726 real& ssig2, real& csig2,
727 real& eps, real& domg12,
728 bool diffp, real& dlam12,
730 real C1a[], real C2a[], real C3a[])
const {
732 if (sbet1 == 0 && calp1 == 0)
739 salp0 = salp1 * cbet1,
742 real somg1, comg1, somg2, comg2, omg12, lam12;
745 ssig1 = sbet1; somg1 = salp0 * sbet1;
746 csig1 = comg1 = calp1 * cbet1;
754 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
759 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ?
762 (cbet2 - cbet1) * (cbet1 + cbet2) :
763 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
767 ssig2 = sbet2; somg2 = salp0 * sbet2;
768 csig2 = comg2 = calp2 * cbet2;
773 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
774 csig1 * csig2 + ssig1 * ssig2);
777 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, real(0)),
778 comg1 * comg2 + somg1 * somg2);
781 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
783 B312 = (SinCosSeries(
true, ssig2, csig2, C3a, nC3_-1) -
784 SinCosSeries(
true, ssig1, csig1, C3a, nC3_-1));
786 domg12 = salp0 * h0 * (sig12 + B312);
787 lam12 = omg12 + domg12;
791 dlam12 = - 2 * _f1 * dn1 / sbet1;
794 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
795 cbet1, cbet2, dummy, dlam12, dummy,
796 false, dummy, dummy, C1a, C2a);
797 dlam12 *= _f1 / (calp2 * cbet2);
809 void Geodesic::C3f(real eps, real c[])
const {
814 for (
int l = 1; l < nC3_; ++l) {
815 int m = nC3_ - l - 1;
823 void Geodesic::C4f(real eps, real c[])
const {
828 for (
int l = 0; l < nC4_; ++l) {
829 int m = nC4_ - l - 1;
863 #if GEOGRAPHICLIB_GEODESIC_ORDER/2 == 1 864 static const real coeff[] = {
868 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 2 869 static const real coeff[] = {
873 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 3 874 static const real coeff[] = {
878 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 4 879 static const real coeff[] = {
881 25, 64, 256, 4096, 0, 16384,
884 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 886 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) == nA1_/2 + 2,
887 "Coefficient array size mismatch in A1m1f");
890 return (t + eps) / (1 - eps);
894 void Geodesic::C1f(real eps, real c[]) {
896 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3 897 static const real coeff[] = {
905 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4 906 static const real coeff[] = {
916 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5 917 static const real coeff[] = {
929 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6 930 static const real coeff[] = {
944 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7 945 static const real coeff[] = {
947 19, -64, 384, -1024, 2048,
961 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8 962 static const real coeff[] = {
964 19, -64, 384, -1024, 2048,
966 7, -18, 128, -256, 4096,
970 -11, 96, -160, 16384,
981 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 983 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) ==
984 (nC1_*nC1_ + 7*nC1_ - 2*(nC1_/2)) / 4,
985 "Coefficient array size mismatch in C1f");
990 for (
int l = 1; l <= nC1_; ++l) {
991 int m = (nC1_ - l) / 2;
992 c[l] = d *
Math::polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1000 void Geodesic::C1pf(real eps, real c[]) {
1002 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3 1003 static const real coeff[] = {
1011 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4 1012 static const real coeff[] = {
1022 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5 1023 static const real coeff[] = {
1025 205, -432, 768, 1536,
1035 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6 1036 static const real coeff[] = {
1038 205, -432, 768, 1536,
1040 4005, -4736, 3840, 12288,
1050 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7 1051 static const real coeff[] = {
1053 -4879, 9840, -20736, 36864, 73728,
1055 4005, -4736, 3840, 12288,
1057 8703, -7200, 3712, 12288,
1061 -141115, 41604, 92160,
1067 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8 1068 static const real coeff[] = {
1070 -4879, 9840, -20736, 36864, 73728,
1072 -86171, 120150, -142080, 115200, 368640,
1074 8703, -7200, 3712, 12288,
1076 1082857, -688608, 258720, 737280,
1078 -141115, 41604, 92160,
1080 -2200311, 533134, 860160,
1084 109167851, 82575360,
1087 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 1089 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) ==
1090 (nC1p_*nC1p_ + 7*nC1p_ - 2*(nC1p_/2)) / 4,
1091 "Coefficient array size mismatch in C1pf");
1096 for (
int l = 1; l <= nC1p_; ++l) {
1097 int m = (nC1p_ - l) / 2;
1098 c[l] = d *
Math::polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1108 #if GEOGRAPHICLIB_GEODESIC_ORDER/2 == 1 1109 static const real coeff[] = {
1113 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 2 1114 static const real coeff[] = {
1118 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 3 1119 static const real coeff[] = {
1123 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 4 1124 static const real coeff[] = {
1126 1225, 1600, 2304, 4096, 0, 16384,
1129 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 1131 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) == nA2_/2 + 2,
1132 "Coefficient array size mismatch in A2m1f");
1135 return t * (1 - eps) - eps;
1139 void Geodesic::C2f(real eps, real c[]) {
1141 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3 1142 static const real coeff[] = {
1150 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4 1151 static const real coeff[] = {
1161 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5 1162 static const real coeff[] = {
1174 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6 1175 static const real coeff[] = {
1189 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7 1190 static const real coeff[] = {
1192 41, 64, 128, 1024, 2048,
1206 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8 1207 static const real coeff[] = {
1209 41, 64, 128, 1024, 2048,
1211 47, 70, 128, 768, 4096,
1215 133, 224, 1120, 16384,
1226 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 1228 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) ==
1229 (nC2_*nC2_ + 7*nC2_ - 2*(nC2_/2)) / 4,
1230 "Coefficient array size mismatch in C2f");
1235 for (
int l = 1; l <= nC2_; ++l) {
1236 int m = (nC2_ - l) / 2;
1237 c[l] = d *
Math::polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1245 void Geodesic::A3coeff() {
1247 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3 1248 static const real coeff[] = {
1256 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4 1257 static const real coeff[] = {
1267 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5 1268 static const real coeff[] = {
1280 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6 1281 static const real coeff[] = {
1295 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7 1296 static const real coeff[] = {
1312 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8 1313 static const real coeff[] = {
1321 -5, -20, -4, -6, 128,
1332 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 1334 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) ==
1335 (nA3_*nA3_ + 7*nA3_ - 2*(nA3_/2)) / 4,
1336 "Coefficient array size mismatch in A3f");
1338 for (
int j = nA3_ - 1; j >= 0; --j) {
1339 int m = min(nA3_ - j - 1, j);
1340 _A3x[k++] =
Math::polyval(m, coeff + o, _n) / coeff[o + m + 1];
1347 void Geodesic::C3coeff() {
1349 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3 1350 static const real coeff[] = {
1358 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4 1359 static const real coeff[] = {
1375 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5 1376 static const real coeff[] = {
1398 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6 1399 static const real coeff[] = {
1431 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7 1432 static const real coeff[] = {
1476 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8 1477 static const real coeff[] = {
1511 10, -6, -10, 9, 384,
1521 -7, 20, -28, 14, 1024,
1536 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 1538 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) ==
1539 ((nC3_-1)*(nC3_*nC3_ + 7*nC3_ - 2*(nC3_/2)))/8,
1540 "Coefficient array size mismatch in C3coeff");
1542 for (
int l = 1; l < nC3_; ++l) {
1543 for (
int j = nC3_ - 1; j >= l; --j) {
1544 int m = min(nC3_ - j - 1, j);
1545 _C3x[k++] =
Math::polyval(m, coeff + o, _n) / coeff[o + m + 1];
1552 void Geodesic::C4coeff() {
1554 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3 1555 static const real coeff[] = {
1569 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4 1570 static const real coeff[] = {
1578 4, 24, -84, 210, 315,
1592 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5 1593 static const real coeff[] = {
1599 1088, -352, -66, 3465,
1601 48, -352, 528, -231, 1155,
1603 16, 44, 264, -924, 2310, 3465,
1609 -896, 704, -198, 10395,
1611 -48, 352, -528, 231, 10395,
1617 320, -352, 132, 17325,
1625 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6 1626 static const real coeff[] = {
1632 -224, -4784, 1573, 45045,
1634 -10656, 14144, -4576, -858, 45045,
1636 64, 624, -4576, 6864, -3003, 15015,
1638 100, 208, 572, 3432, -12012, 30030, 45045,
1644 5792, 1040, -1287, 135135,
1646 5952, -11648, 9152, -2574, 135135,
1648 -64, -624, 4576, -6864, 3003, 135135,
1654 -8448, 4992, -1144, 225225,
1656 -1440, 4160, -4576, 1716, 225225,
1662 3584, -3328, 1144, 315315,
1670 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7 1671 static const real coeff[] = {
1677 -4480, 1088, 156, 45045,
1679 10736, -224, -4784, 1573, 45045,
1681 1664, -10656, 14144, -4576, -858, 45045,
1683 16, 64, 624, -4576, 6864, -3003, 15015,
1685 56, 100, 208, 572, 3432, -12012, 30030, 45045,
1691 3840, -2944, 468, 135135,
1693 -10704, 5792, 1040, -1287, 135135,
1695 -768, 5952, -11648, 9152, -2574, 135135,
1697 -16, -64, -624, 4576, -6864, 3003, 135135,
1703 1664, 1856, -936, 225225,
1705 6784, -8448, 4992, -1144, 225225,
1707 128, -1440, 4160, -4576, 1716, 225225,
1713 -2048, 1024, -208, 105105,
1715 -1792, 3584, -3328, 1144, 315315,
1721 3072, -2560, 832, 405405,
1729 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8 1730 static const real coeff[] = {
1736 20960, -7888, 4947, 765765,
1738 12480, -76160, 18496, 2652, 765765,
1740 -154048, 182512, -3808, -81328, 26741, 765765,
1742 3232, 28288, -181152, 240448, -77792, -14586, 765765,
1744 96, 272, 1088, 10608, -77792, 116688, -51051, 255255,
1746 588, 952, 1700, 3536, 9724, 58344, -204204, 510510, 765765,
1752 -39840, 1904, 255, 2297295,
1754 52608, 65280, -50048, 7956, 2297295,
1756 103744, -181968, 98464, 17680, -21879, 2297295,
1758 -1344, -13056, 101184, -198016, 155584, -43758, 2297295,
1760 -96, -272, -1088, -10608, 77792, -116688, 51051, 2297295,
1764 -928, -612, 3828825,
1766 64256, -28288, 2856, 3828825,
1768 -126528, 28288, 31552, -15912, 3828825,
1770 -41472, 115328, -143616, 84864, -19448, 3828825,
1772 160, 2176, -24480, 70720, -77792, 29172, 3828825,
1776 -16384, 1088, 5360355,
1778 -2560, 30464, -11560, 5360355,
1780 35840, -34816, 17408, -3536, 1786785,
1782 7168, -30464, 60928, -56576, 19448, 5360355,
1786 26624, -8704, 6891885,
1788 -77824, 34816, -6528, 6891885,
1790 -32256, 52224, -43520, 14144, 6891885,
1794 24576, -4352, 8423415,
1796 45056, -34816, 10880, 8423415,
1800 -28672, 8704, 9954945,
1805 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER" 1807 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(real) ==
1808 (nC4_ * (nC4_ + 1) * (nC4_ + 5)) / 6,
1809 "Coefficient array size mismatch in C4coeff");
1811 for (
int l = 0; l < nC4_; ++l) {
1812 for (
int j = nC4_ - 1; j >= l; --j) {
1813 int m = nC4_ - j - 1;
1814 _C4x[k++] =
Math::polyval(m, coeff + o, _n) / coeff[o + m + 1];
static T AngNormalize(T x)
Header for GeographicLib::GeodesicLine class.
Math::real GenDirect(real lat1, real lon1, real azi1, bool arcmode, real s12_a12, unsigned outmask, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
static const Geodesic & WGS84()
static bool isfinite(T x)
Mathematical functions needed by GeographicLib.
static void norm(T &x, T &y)
#define GEOGRAPHICLIB_VOLATILE
Header for GeographicLib::Geodesic class.
friend class GeodesicLine
static T atan2d(T y, T x)
static T polyval(int N, const T p[], T x)
Namespace for GeographicLib.
static T AngDiff(T x, T y)
Math::real GenInverse(real lat1, real lon1, real lat2, real lon2, unsigned outmask, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21, real &S12) const
Exception handling for GeographicLib.
#define GEOGRAPHICLIB_PANIC