34 # pragma warning (disable: 4701)
37 namespace GeographicLib {
44 const Math::real Geodesic::tiny_ = sqrt(numeric_limits<real>::min());
45 const Math::real Geodesic::tol0_ = numeric_limits<real>::epsilon();
49 const Math::real Geodesic::tol1_ = 200 * tol0_;
50 const Math::real Geodesic::tol2_ = sqrt(tol0_);
52 const Math::real Geodesic::tolb_ = tol0_ * tol2_;
53 const Math::real Geodesic::xthresh_ = 1000 * tol2_;
57 , _f(f <= 1 ? f : 1/f)
60 , _ep2(_e2 /
Math::sq(_f1))
65 (_e2 > 0 ?
Math::atanh(sqrt(_e2)) : atan(sqrt(-_e2))) /
77 , _etol2(0.1 * tol2_ /
78 sqrt( max(real(0.001), abs(_f)) * min(real(1), 1 - _f/2) / 2 ))
90 Constants::WGS84_f<real>());
94 const real c[],
int n)
throw() {
102 ar = 2 * (cosx - sinx) * (cosx + sinx),
103 y0 = n & 1 ? *--c : 0, y1 = 0;
108 y1 = ar * y0 - y1 + *--c;
109 y0 = ar * y1 - y0 + *--c;
112 ? 2 * sinx * cosx * y0
122 bool arcmode, real s12_a12,
unsigned outmask,
123 real& lat2, real& lon2, real& azi2,
124 real& s12, real& m12, real& M12, real& M21,
125 real& S12)
const throw() {
128 outmask | (arcmode ? NONE : DISTANCE_IN))
130 GenPosition(arcmode, s12_a12, outmask,
131 lat2, lon2, azi2, s12, m12, M12, M21, S12);
136 real& s12, real& azi1, real& azi2,
137 real& m12, real& M12, real& M21, real& S12)
146 lon12 = AngRound(lon12);
148 int lonsign = lon12 >= 0 ? 1 : -1;
151 lat1 = AngRound(lat1);
152 lat2 = AngRound(lat2);
154 int swapp = abs(lat1) >= abs(lat2) ? 1 : -1;
160 int latsign = lat1 < 0 ? 1 : -1;
175 real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x;
177 phi = lat1 * Math::degree<real>();
179 sbet1 = _f1 * sin(phi);
180 cbet1 = lat1 == -90 ? tiny_ : cos(phi);
181 SinCosNorm(sbet1, cbet1);
183 phi = lat2 * Math::degree<real>();
185 sbet2 = _f1 * sin(phi);
186 cbet2 = abs(lat2) == 90 ? tiny_ : cos(phi);
187 SinCosNorm(sbet2, cbet2);
197 if (cbet1 < -sbet1) {
199 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
201 if (abs(sbet2) == -sbet1)
206 dn1 = sqrt(1 + _ep2 *
Math::sq(sbet1)),
207 dn2 = sqrt(1 + _ep2 *
Math::sq(sbet2));
210 lam12 = lon12 * Math::degree<real>(),
211 slam12 = abs(lon12) == 180 ? 0 : sin(lam12),
214 real a12, sig12, calp1, salp1, calp2, salp2;
216 real C1a[nC1_ + 1], C2a[nC2_ + 1], C3a[nC3_];
218 bool meridian = lat1 == -90 || slam12 == 0;
225 calp1 = clam12; salp1 = slam12;
226 calp2 = 1; salp2 = 0;
230 ssig1 = sbet1, csig1 = calp1 * cbet1,
231 ssig2 = sbet2, csig2 = calp2 * cbet2;
234 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2,
real(0)),
235 csig1 * csig2 + ssig1 * ssig2);
238 Lengths(_n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
239 cbet1, cbet2, s12x, m12x, dummy,
240 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a);
249 if (sig12 < 1 || m12x >= 0) {
252 a12 = sig12 / Math::degree<real>();
262 (_f <= 0 || lam12 <= Math::pi<real>() - _f * Math::pi<real>())) {
265 calp1 = calp2 = 0; salp1 = salp2 = 1;
267 sig12 = omg12 = lam12 / _f1;
268 m12x = _b * sin(sig12);
269 if (outmask & GEODESICSCALE)
270 M12 = M21 = cos(sig12);
273 }
else if (!meridian) {
280 sig12 = InverseStart(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
282 salp1, calp1, salp2, calp2, dnm,
287 s12x = sig12 * _b * dnm;
288 m12x =
Math::sq(dnm) * _b * sin(sig12 / dnm);
289 if (outmask & GEODESICSCALE)
290 M12 = M21 = cos(sig12 / dnm);
291 a12 = sig12 / Math::degree<real>();
292 omg12 = lam12 / (_f1 * dnm);
306 real ssig1, csig1, ssig2, csig2, eps;
309 real salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
310 for (
bool tripn =
false, tripb =
false; numit < maxit2_; ++numit) {
314 real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
315 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
316 eps, omg12, numit < maxit1_, dv, C1a, C2a, C3a)
320 if (tripb || !(abs(v) >= (tripn ? 8 : 2) * tol0_))
break;
322 if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
323 { salp1b = salp1; calp1b = calp1; }
324 else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
325 { salp1a = salp1; calp1a = calp1; }
326 if (numit < maxit1_ && dv > 0) {
330 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
331 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
332 if (nsalp1 > 0 && abs(dalp1) < Math::pi<real>()) {
333 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
335 SinCosNorm(salp1, calp1);
339 tripn = abs(v) <= 16 * tol0_;
351 salp1 = (salp1a + salp1b)/2;
352 calp1 = (calp1a + calp1b)/2;
353 SinCosNorm(salp1, calp1);
355 tripb = (abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
356 abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
360 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
361 cbet1, cbet2, s12x, m12x, dummy,
362 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a);
366 a12 = sig12 / Math::degree<real>();
367 omg12 = lam12 - omg12;
371 if (outmask & DISTANCE)
374 if (outmask & REDUCEDLENGTH)
377 if (outmask & AREA) {
380 salp0 = salp1 * cbet1,
383 if (calp0 != 0 && salp0 != 0) {
386 ssig1 = sbet1, csig1 = calp1 * cbet1,
387 ssig2 = sbet2, csig2 = calp2 * cbet2,
389 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
391 A4 =
Math::sq(_a) * calp0 * salp0 * _e2;
392 SinCosNorm(ssig1, csig1);
393 SinCosNorm(ssig2, csig2);
397 B41 = SinCosSeries(
false, ssig1, csig1, C4a, nC4_),
398 B42 = SinCosSeries(
false, ssig2, csig2, C4a, nC4_);
399 S12 = A4 * (B42 - B41);
405 omg12 <
real(0.75) * Math::pi<real>() &&
406 sbet2 - sbet1 <
real(1.75)) {
411 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
412 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
413 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
414 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
418 salp12 = salp2 * calp1 - calp2 * salp1,
419 calp12 = calp2 * calp1 + salp2 * salp1;
424 if (salp12 == 0 && calp12 < 0) {
425 salp12 = tiny_ * calp1;
428 alp12 = atan2(salp12, calp12);
431 S12 *= swapp * lonsign * latsign;
440 if (outmask & GEODESICSCALE)
444 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
445 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
447 if (outmask & AZIMUTH) {
449 azi1 = 0 - atan2(-salp1, calp1) / Math::degree<real>();
450 azi2 = 0 - atan2(-salp2, calp2) / Math::degree<real>();
457 void Geodesic::Lengths(
real eps,
real sig12,
464 real C1a[],
real C2a[])
const throw() {
471 AB1 = (1 + A1m1) * (SinCosSeries(
true, ssig2, csig2, C1a, nC1_) -
472 SinCosSeries(
true, ssig1, csig1, C1a, nC1_)),
474 AB2 = (1 + A2m1) * (SinCosSeries(
true, ssig2, csig2, C2a, nC2_) -
475 SinCosSeries(
true, ssig1, csig1, C2a, nC2_));
477 real J12 = m0 * sig12 + (AB1 - AB2);
481 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12;
483 s12b = (1 + A1m1) * sig12 + AB1;
485 real csig12 = csig1 * csig2 + ssig1 * ssig2;
486 real t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
487 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
488 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
500 if ( !(q == 0 && r <= 0) ) {
509 disc = S * (S + 2 * r3);
516 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
520 u += T + (T != 0 ? r2 / T : 0);
523 real ang = atan2(sqrt(-disc), -(S + r3));
526 u += 2 * r * cos(ang / 3);
531 uv = u < 0 ? q / (v - u) : u + v,
532 w = (uv - q) / (2 * v);
535 k = uv / (sqrt(uv +
Math::sq(w)) + w);
553 real C1a[],
real C2a[])
const throw() {
560 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
561 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
562 #if defined(__GNUC__) && __GNUC__ == 4 && \
563 (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
572 volatile real xx1 = sbet2 * cbet1;
573 volatile real xx2 = cbet2 * sbet1;
577 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
579 bool shortline = cbet12 >= 0 && sbet12 <
real(0.5) &&
580 cbet2 * lam12 <
real(0.5);
586 sbetm2 /= sbetm2 +
Math::sq(cbet1 + cbet2);
587 dnm = sqrt(1 + _ep2 * sbetm2);
590 real somg12 = sin(omg12), comg12 = cos(omg12);
592 salp1 = cbet2 * somg12;
593 calp1 = comg12 >= 0 ?
594 sbet12 + cbet2 * sbet1 *
Math::sq(somg12) / (1 + comg12) :
595 sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
599 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
601 if (shortline && ssig12 < _etol2) {
603 salp2 = cbet1 * somg12;
604 calp2 = sbet12 - cbet1 * sbet2 *
605 (comg12 >= 0 ?
Math::sq(somg12) / (1 + comg12) : 1 - comg12);
606 SinCosNorm(salp2, calp2);
608 sig12 = atan2(ssig12, csig12);
609 }
else if (abs(_n) >
real(0.1) ||
611 ssig12 >= 6 * abs(_n) * Math::pi<real>() *
Math::sq(cbet1)) {
616 real y, lamscale, betscale;
626 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
627 lamscale = _f * cbet1 * A3f(eps) * Math::pi<real>();
629 betscale = lamscale * cbet1;
631 x = (lam12 - Math::pi<real>()) / lamscale;
632 y = sbet12a / betscale;
636 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
637 bet12a = atan2(sbet12a, cbet12a);
638 real m12b, m0, dummy;
641 Lengths(_n, Math::pi<real>() + bet12a,
642 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
643 cbet1, cbet2, dummy, m12b, m0,
false,
644 dummy, dummy, C1a, C2a);
645 x = -1 + m12b / (cbet1 * cbet2 * m0 * Math::pi<real>());
646 betscale = x < -
real(0.01) ? sbet12a / x :
647 -_f *
Math::sq(cbet1) * Math::pi<real>();
648 lamscale = betscale / cbet1;
649 y = (lam12 - Math::pi<real>()) / lamscale;
652 if (y > -tol1_ && x > -1 - xthresh_) {
658 calp1 = max(
real(x > -tol1_ ? 0 : -1),
real(x));
696 real k = Astroid(x, y);
698 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
699 somg12 = sin(omg12a); comg12 = -cos(omg12a);
701 salp1 = cbet2 * somg12;
702 calp1 = sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
706 SinCosNorm(salp1, calp1);
708 salp1 = 1; calp1 = 0;
721 bool diffp,
real& dlam12,
726 if (sbet1 == 0 && calp1 == 0)
733 salp0 = salp1 * cbet1,
736 real somg1, comg1, somg2, comg2, omg12, lam12;
739 ssig1 = sbet1; somg1 = salp0 * sbet1;
740 csig1 = comg1 = calp1 * cbet1;
741 SinCosNorm(ssig1, csig1);
748 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
753 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ?
756 (cbet2 - cbet1) * (cbet1 + cbet2) :
757 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
761 ssig2 = sbet2; somg2 = salp0 * sbet2;
762 csig2 = comg2 = calp2 * cbet2;
763 SinCosNorm(ssig2, csig2);
767 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2,
real(0)),
768 csig1 * csig2 + ssig1 * ssig2);
771 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2,
real(0)),
772 comg1 * comg2 + somg1 * somg2);
775 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
777 B312 = (SinCosSeries(
true, ssig2, csig2, C3a, nC3_-1) -
778 SinCosSeries(
true, ssig1, csig1, C3a, nC3_-1));
780 domg12 = salp0 * h0 * (sig12 + B312);
781 lam12 = omg12 + domg12;
785 dlam12 = - 2 * _f1 * dn1 / sbet1;
788 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
789 cbet1, cbet2, dummy, dlam12, dummy,
790 false, dummy, dummy, C1a, C2a);
791 dlam12 *= _f1 / (calp2 * cbet2);
801 for (
int i = nA3x_; i > 0; )
802 v = eps * v + _A3x[--i];
806 void Geodesic::C3f(
real eps,
real c[])
const throw() {
809 for (
int j = nC3x_, k = nC3_ - 1; k > 0; ) {
811 for (
int i = nC3_ - k; i > 0; --i) {
812 t = eps * t + _C3x[--j];
818 for (
int k = 1; k < nC3_; ) {
824 void Geodesic::C4f(
real eps,
real c[])
const throw() {
827 for (
int j = nC4x_, k = nC4_; k > 0; ) {
829 for (
int i = nC4_ - k + 1; i > 0; --i)
830 t = eps * t + _C4x[--j];
835 for (
int k = 1; k < nC4_; ) {
856 t = eps2*(eps2+16)/64;
859 t = eps2*(eps2*(eps2+4)+64)/256;
862 t = eps2*(eps2*(eps2*(25*eps2+64)+256)+4096)/16384;
868 return (t + eps) / (1 - eps);
872 void Geodesic::C1f(
real eps,
real c[])
throw() {
888 c[1] = d*(3*eps2-8)/16;
895 c[1] = d*(3*eps2-8)/16;
897 c[2] = d*(eps2-2)/32;
904 c[1] = d*((6-eps2)*eps2-16)/32;
906 c[2] = d*(eps2-2)/32;
908 c[3] = d*(9*eps2-16)/768;
915 c[1] = d*((6-eps2)*eps2-16)/32;
917 c[2] = d*((64-9*eps2)*eps2-128)/2048;
919 c[3] = d*(9*eps2-16)/768;
921 c[4] = d*(3*eps2-5)/512;
928 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048;
930 c[2] = d*((64-9*eps2)*eps2-128)/2048;
932 c[3] = d*((72-9*eps2)*eps2-128)/6144;
934 c[4] = d*(3*eps2-5)/512;
936 c[5] = d*(35*eps2-56)/10240;
943 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048;
945 c[2] = d*(eps2*(eps2*(7*eps2-18)+128)-256)/4096;
947 c[3] = d*((72-9*eps2)*eps2-128)/6144;
949 c[4] = d*((96-11*eps2)*eps2-160)/16384;
951 c[5] = d*(35*eps2-56)/10240;
953 c[6] = d*(9*eps2-14)/4096;
957 c[8] = -429*d/262144;
965 void Geodesic::C1pf(
real eps,
real c[])
throw() {
981 c[1] = d*(16-9*eps2)/32;
988 c[1] = d*(16-9*eps2)/32;
990 c[2] = d*(30-37*eps2)/96;
997 c[1] = d*(eps2*(205*eps2-432)+768)/1536;
999 c[2] = d*(30-37*eps2)/96;
1001 c[3] = d*(116-225*eps2)/384;
1008 c[1] = d*(eps2*(205*eps2-432)+768)/1536;
1010 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
1012 c[3] = d*(116-225*eps2)/384;
1014 c[4] = d*(2695-7173*eps2)/7680;
1018 c[6] = 38081*d/61440;
1021 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728;
1023 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
1025 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288;
1027 c[4] = d*(2695-7173*eps2)/7680;
1029 c[5] = d*(41604-141115*eps2)/92160;
1031 c[6] = 38081*d/61440;
1033 c[7] = 459485*d/516096;
1036 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728;
1038 c[2] = d*(eps2*((120150-86171*eps2)*eps2-142080)+115200)/368640;
1040 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288;
1042 c[4] = d*(eps2*(1082857*eps2-688608)+258720)/737280;
1044 c[5] = d*(41604-141115*eps2)/92160;
1046 c[6] = d*(533134-2200311*eps2)/860160;
1048 c[7] = 459485*d/516096;
1050 c[8] = 109167851*d/82575360;
1053 STATIC_ASSERT(nC1p_ >= 0 && nC1p_ <= 8,
"Bad value of nC1p_");
1070 t = eps2*(9*eps2+16)/64;
1073 t = eps2*(eps2*(25*eps2+36)+64)/256;
1076 t = eps2*(eps2*(eps2*(1225*eps2+1600)+2304)+4096)/16384;
1082 return t * (1 - eps) - eps;
1086 void Geodesic::C2f(
real eps,
real c[])
throw() {
1102 c[1] = d*(eps2+8)/16;
1109 c[1] = d*(eps2+8)/16;
1111 c[2] = d*(eps2+6)/32;
1118 c[1] = d*(eps2*(eps2+2)+16)/32;
1120 c[2] = d*(eps2+6)/32;
1122 c[3] = d*(15*eps2+80)/768;
1129 c[1] = d*(eps2*(eps2+2)+16)/32;
1131 c[2] = d*(eps2*(35*eps2+64)+384)/2048;
1133 c[3] = d*(15*eps2+80)/768;
1135 c[4] = d*(7*eps2+35)/512;
1142 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048;
1144 c[2] = d*(eps2*(35*eps2+64)+384)/2048;
1146 c[3] = d*(eps2*(69*eps2+120)+640)/6144;
1148 c[4] = d*(7*eps2+35)/512;
1150 c[5] = d*(105*eps2+504)/10240;
1157 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048;
1159 c[2] = d*(eps2*(eps2*(47*eps2+70)+128)+768)/4096;
1161 c[3] = d*(eps2*(69*eps2+120)+640)/6144;
1163 c[4] = d*(eps2*(133*eps2+224)+1120)/16384;
1165 c[5] = d*(105*eps2+504)/10240;
1167 c[6] = d*(33*eps2+154)/4096;
1171 c[8] = 6435*d/262144;
1179 void Geodesic::A3coeff() throw() {
1188 _A3x[1] = -1/
real(2);
1193 _A3x[2] = -1/
real(4);
1198 _A3x[2] = (-_n-2)/8;
1199 _A3x[3] = -1/
real(16);
1204 _A3x[2] = (_n*(3*_n-1)-2)/8;
1205 _A3x[3] = (-3*_n-1)/16;
1206 _A3x[4] = -3/
real(64);
1211 _A3x[2] = (_n*(3*_n-1)-2)/8;
1212 _A3x[3] = ((-_n-3)*_n-1)/16;
1213 _A3x[4] = (-2*_n-3)/64;
1214 _A3x[5] = -3/
real(128);
1219 _A3x[2] = (_n*(3*_n-1)-2)/8;
1220 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16;
1221 _A3x[4] = ((-10*_n-2)*_n-3)/64;
1222 _A3x[5] = (-5*_n-3)/128;
1223 _A3x[6] = -5/
real(256);
1228 _A3x[2] = (_n*(3*_n-1)-2)/8;
1229 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16;
1230 _A3x[4] = (_n*((-5*_n-20)*_n-4)-6)/128;
1231 _A3x[5] = ((-5*_n-10)*_n-6)/256;
1232 _A3x[6] = (-15*_n-20)/1024;
1233 _A3x[7] = -25/
real(2048);
1241 void Geodesic::C3coeff() throw() {
1248 _C3x[0] = 1/
real(4);
1252 _C3x[1] = 1/
real(8);
1253 _C3x[2] = 1/
real(16);
1257 _C3x[1] = 1/
real(8);
1258 _C3x[2] = 3/
real(64);
1259 _C3x[3] = (2-3*_n)/32;
1260 _C3x[4] = 3/
real(64);
1261 _C3x[5] = 5/
real(192);
1265 _C3x[1] = (1-_n*_n)/8;
1266 _C3x[2] = (3*_n+3)/64;
1267 _C3x[3] = 5/
real(128);
1268 _C3x[4] = ((_n-3)*_n+2)/32;
1269 _C3x[5] = (3-2*_n)/64;
1270 _C3x[6] = 3/
real(128);
1271 _C3x[7] = (5-9*_n)/192;
1272 _C3x[8] = 3/
real(128);
1273 _C3x[9] = 7/
real(512);
1277 _C3x[1] = (1-_n*_n)/8;
1278 _C3x[2] = ((3-_n)*_n+3)/64;
1279 _C3x[3] = (2*_n+5)/128;
1280 _C3x[4] = 3/
real(128);
1281 _C3x[5] = ((_n-3)*_n+2)/32;
1282 _C3x[6] = ((-3*_n-2)*_n+3)/64;
1283 _C3x[7] = (_n+3)/128;
1284 _C3x[8] = 5/
real(256);
1285 _C3x[9] = (_n*(5*_n-9)+5)/192;
1286 _C3x[10] = (9-10*_n)/384;
1287 _C3x[11] = 7/
real(512);
1288 _C3x[12] = (7-14*_n)/512;
1289 _C3x[13] = 7/
real(512);
1290 _C3x[14] = 21/
real(2560);
1294 _C3x[1] = (1-_n*_n)/8;
1295 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64;
1296 _C3x[3] = (_n*(2*_n+2)+5)/128;
1297 _C3x[4] = (11*_n+12)/512;
1298 _C3x[5] = 21/
real(1024);
1299 _C3x[6] = ((_n-3)*_n+2)/32;
1300 _C3x[7] = (_n*(_n*(2*_n-3)-2)+3)/64;
1301 _C3x[8] = ((2-9*_n)*_n+6)/256;
1302 _C3x[9] = (_n+5)/256;
1303 _C3x[10] = 27/
real(2048);
1304 _C3x[11] = (_n*((5-_n)*_n-9)+5)/192;
1305 _C3x[12] = ((-6*_n-10)*_n+9)/384;
1306 _C3x[13] = (21-4*_n)/1536;
1307 _C3x[14] = 3/
real(256);
1308 _C3x[15] = (_n*(10*_n-14)+7)/512;
1309 _C3x[16] = (7-10*_n)/512;
1310 _C3x[17] = 9/
real(1024);
1311 _C3x[18] = (21-45*_n)/2560;
1312 _C3x[19] = 9/
real(1024);
1313 _C3x[20] = 11/
real(2048);
1317 _C3x[1] = (1-_n*_n)/8;
1318 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64;
1319 _C3x[3] = (_n*((2-2*_n)*_n+2)+5)/128;
1320 _C3x[4] = (_n*(3*_n+11)+12)/512;
1321 _C3x[5] = (10*_n+21)/1024;
1322 _C3x[6] = 243/
real(16384);
1323 _C3x[7] = ((_n-3)*_n+2)/32;
1324 _C3x[8] = (_n*(_n*(2*_n-3)-2)+3)/64;
1325 _C3x[9] = (_n*((-6*_n-9)*_n+2)+6)/256;
1326 _C3x[10] = ((1-2*_n)*_n+5)/256;
1327 _C3x[11] = (69*_n+108)/8192;
1328 _C3x[12] = 187/
real(16384);
1329 _C3x[13] = (_n*((5-_n)*_n-9)+5)/192;
1330 _C3x[14] = (_n*(_n*(10*_n-6)-10)+9)/384;
1331 _C3x[15] = ((-77*_n-8)*_n+42)/3072;
1332 _C3x[16] = (12-_n)/1024;
1333 _C3x[17] = 139/
real(16384);
1334 _C3x[18] = (_n*((20-7*_n)*_n-28)+14)/1024;
1335 _C3x[19] = ((-7*_n-40)*_n+28)/2048;
1336 _C3x[20] = (72-43*_n)/8192;
1337 _C3x[21] = 127/
real(16384);
1338 _C3x[22] = (_n*(75*_n-90)+42)/5120;
1339 _C3x[23] = (9-15*_n)/1024;
1340 _C3x[24] = 99/
real(16384);
1341 _C3x[25] = (44-99*_n)/8192;
1342 _C3x[26] = 99/
real(16384);
1343 _C3x[27] = 429/
real(114688);
1353 void Geodesic::C4coeff() throw() {
1358 _C4x[0] = 2/
real(3);
1361 _C4x[0] = (10-4*_n)/15;
1362 _C4x[1] = -1/
real(5);
1363 _C4x[2] = 1/
real(45);
1366 _C4x[0] = (_n*(8*_n-28)+70)/105;
1367 _C4x[1] = (16*_n-7)/35;
1368 _C4x[2] = -2/
real(105);
1369 _C4x[3] = (7-16*_n)/315;
1370 _C4x[4] = -2/
real(105);
1371 _C4x[5] = 4/
real(525);
1374 _C4x[0] = (_n*(_n*(4*_n+24)-84)+210)/315;
1375 _C4x[1] = ((48-32*_n)*_n-21)/105;
1376 _C4x[2] = (-32*_n-6)/315;
1377 _C4x[3] = 11/
real(315);
1378 _C4x[4] = (_n*(32*_n-48)+21)/945;
1379 _C4x[5] = (64*_n-18)/945;
1380 _C4x[6] = -1/
real(105);
1381 _C4x[7] = (12-32*_n)/1575;
1382 _C4x[8] = -8/
real(1575);
1383 _C4x[9] = 8/
real(2205);
1386 _C4x[0] = (_n*(_n*(_n*(16*_n+44)+264)-924)+2310)/3465;
1387 _C4x[1] = (_n*(_n*(48*_n-352)+528)-231)/1155;
1388 _C4x[2] = (_n*(1088*_n-352)-66)/3465;
1389 _C4x[3] = (121-368*_n)/3465;
1390 _C4x[4] = 4/
real(1155);
1391 _C4x[5] = (_n*((352-48*_n)*_n-528)+231)/10395;
1392 _C4x[6] = ((704-896*_n)*_n-198)/10395;
1393 _C4x[7] = (80*_n-99)/10395;
1394 _C4x[8] = 4/
real(1155);
1395 _C4x[9] = (_n*(320*_n-352)+132)/17325;
1396 _C4x[10] = (384*_n-88)/17325;
1397 _C4x[11] = -8/
real(1925);
1398 _C4x[12] = (88-256*_n)/24255;
1399 _C4x[13] = -16/
real(8085);
1400 _C4x[14] = 64/
real(31185);
1403 _C4x[0] = (_n*(_n*(_n*(_n*(100*_n+208)+572)+3432)-12012)+30030)/45045;
1404 _C4x[1] = (_n*(_n*(_n*(64*_n+624)-4576)+6864)-3003)/15015;
1405 _C4x[2] = (_n*((14144-10656*_n)*_n-4576)-858)/45045;
1406 _C4x[3] = ((-224*_n-4784)*_n+1573)/45045;
1407 _C4x[4] = (1088*_n+156)/45045;
1408 _C4x[5] = 97/
real(15015);
1409 _C4x[6] = (_n*(_n*((-64*_n-624)*_n+4576)-6864)+3003)/135135;
1410 _C4x[7] = (_n*(_n*(5952*_n-11648)+9152)-2574)/135135;
1411 _C4x[8] = (_n*(5792*_n+1040)-1287)/135135;
1412 _C4x[9] = (468-2944*_n)/135135;
1413 _C4x[10] = 1/
real(9009);
1414 _C4x[11] = (_n*((4160-1440*_n)*_n-4576)+1716)/225225;
1415 _C4x[12] = ((4992-8448*_n)*_n-1144)/225225;
1416 _C4x[13] = (1856*_n-936)/225225;
1417 _C4x[14] = 8/
real(10725);
1418 _C4x[15] = (_n*(3584*_n-3328)+1144)/315315;
1419 _C4x[16] = (1024*_n-208)/105105;
1420 _C4x[17] = -136/
real(63063);
1421 _C4x[18] = (832-2560*_n)/405405;
1422 _C4x[19] = -128/
real(135135);
1423 _C4x[20] = 128/
real(99099);
1426 _C4x[0] = (_n*(_n*(_n*(_n*(_n*(56*_n+100)+208)+572)+3432)-12012)+30030)/
1428 _C4x[1] = (_n*(_n*(_n*(_n*(16*_n+64)+624)-4576)+6864)-3003)/15015;
1429 _C4x[2] = (_n*(_n*(_n*(1664*_n-10656)+14144)-4576)-858)/45045;
1430 _C4x[3] = (_n*(_n*(10736*_n-224)-4784)+1573)/45045;
1431 _C4x[4] = ((1088-4480*_n)*_n+156)/45045;
1432 _C4x[5] = (291-464*_n)/45045;
1433 _C4x[6] = 10/
real(9009);
1434 _C4x[7] = (_n*(_n*(_n*((-16*_n-64)*_n-624)+4576)-6864)+3003)/135135;
1435 _C4x[8] = (_n*(_n*((5952-768*_n)*_n-11648)+9152)-2574)/135135;
1436 _C4x[9] = (_n*((5792-10704*_n)*_n+1040)-1287)/135135;
1437 _C4x[10] = (_n*(3840*_n-2944)+468)/135135;
1438 _C4x[11] = (112*_n+15)/135135;
1439 _C4x[12] = 10/
real(9009);
1440 _C4x[13] = (_n*(_n*(_n*(128*_n-1440)+4160)-4576)+1716)/225225;
1441 _C4x[14] = (_n*(_n*(6784*_n-8448)+4992)-1144)/225225;
1442 _C4x[15] = (_n*(1664*_n+1856)-936)/225225;
1443 _C4x[16] = (168-1664*_n)/225225;
1444 _C4x[17] = -4/
real(25025);
1445 _C4x[18] = (_n*((3584-1792*_n)*_n-3328)+1144)/315315;
1446 _C4x[19] = ((1024-2048*_n)*_n-208)/105105;
1447 _C4x[20] = (1792*_n-680)/315315;
1448 _C4x[21] = 64/
real(315315);
1449 _C4x[22] = (_n*(3072*_n-2560)+832)/405405;
1450 _C4x[23] = (2048*_n-384)/405405;
1451 _C4x[24] = -512/
real(405405);
1452 _C4x[25] = (640-2048*_n)/495495;
1453 _C4x[26] = -256/
real(495495);
1454 _C4x[27] = 512/
real(585585);
1457 _C4x[0] = (_n*(_n*(_n*(_n*(_n*(_n*(588*_n+952)+1700)+3536)+9724)+58344)-
1458 204204)+510510)/765765;
1459 _C4x[1] = (_n*(_n*(_n*(_n*(_n*(96*_n+272)+1088)+10608)-77792)+116688)-
1461 _C4x[2] = (_n*(_n*(_n*(_n*(3232*_n+28288)-181152)+240448)-77792)-14586)/
1463 _C4x[3] = (_n*(_n*((182512-154048*_n)*_n-3808)-81328)+26741)/765765;
1464 _C4x[4] = (_n*(_n*(12480*_n-76160)+18496)+2652)/765765;
1465 _C4x[5] = (_n*(20960*_n-7888)+4947)/765765;
1466 _C4x[6] = (4192*_n+850)/765765;
1467 _C4x[7] = 193/
real(85085);
1468 _C4x[8] = (_n*(_n*(_n*(_n*((-96*_n-272)*_n-1088)-10608)+77792)-116688)+
1470 _C4x[9] = (_n*(_n*(_n*((-1344*_n-13056)*_n+101184)-198016)+155584)-43758)/
1472 _C4x[10] = (_n*(_n*(_n*(103744*_n-181968)+98464)+17680)-21879)/2297295;
1473 _C4x[11] = (_n*(_n*(52608*_n+65280)-50048)+7956)/2297295;
1474 _C4x[12] = ((1904-39840*_n)*_n+255)/2297295;
1475 _C4x[13] = (510-1472*_n)/459459;
1476 _C4x[14] = 349/
real(2297295);
1477 _C4x[15] = (_n*(_n*(_n*(_n*(160*_n+2176)-24480)+70720)-77792)+29172)/
1479 _C4x[16] = (_n*(_n*((115328-41472*_n)*_n-143616)+84864)-19448)/3828825;
1480 _C4x[17] = (_n*((28288-126528*_n)*_n+31552)-15912)/3828825;
1481 _C4x[18] = (_n*(64256*_n-28288)+2856)/3828825;
1482 _C4x[19] = (-928*_n-612)/3828825;
1483 _C4x[20] = 464/
real(1276275);
1484 _C4x[21] = (_n*(_n*(_n*(7168*_n-30464)+60928)-56576)+19448)/5360355;
1485 _C4x[22] = (_n*(_n*(35840*_n-34816)+17408)-3536)/1786785;
1486 _C4x[23] = ((30464-2560*_n)*_n-11560)/5360355;
1487 _C4x[24] = (1088-16384*_n)/5360355;
1488 _C4x[25] = -16/
real(97461);
1489 _C4x[26] = (_n*((52224-32256*_n)*_n-43520)+14144)/6891885;
1490 _C4x[27] = ((34816-77824*_n)*_n-6528)/6891885;
1491 _C4x[28] = (26624*_n-8704)/6891885;
1492 _C4x[29] = 128/
real(2297295);
1493 _C4x[30] = (_n*(45056*_n-34816)+10880)/8423415;
1494 _C4x[31] = (24576*_n-4352)/8423415;
1495 _C4x[32] = -6784/
real(8423415);
1496 _C4x[33] = (8704-28672*_n)/9954945;
1497 _C4x[34] = -1024/
real(3318315);
1498 _C4x[35] = 1024/
real(1640925);
static T AngNormalize(T x)
Header for GeographicLib::GeodesicLine class.
GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
GeographicLib::Math::real real
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
static bool isfinite(T x)
Mathematical functions needed by GeographicLib.
Header for GeographicLib::Geodesic 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
static T AngDiff(T x, T y)
#define STATIC_ASSERT(cond, reason)
static const Geodesic WGS84
Exception handling for GeographicLib.