00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <GeographicLib/Geodesic.hpp>
00030 #include <GeographicLib/GeodesicLine.hpp>
00031
00032 #if defined(_MSC_VER)
00033
00034
00035 # pragma warning (disable: 4701 4127)
00036 #endif
00037
00038 namespace GeographicLib {
00039
00040 using namespace std;
00041
00042 Geodesic::Geodesic(real a, real f)
00043 : maxit2_(maxit1_ + Math::digits() + 10)
00044
00045
00046
00047 , tiny_(sqrt(numeric_limits<real>::min()))
00048 , tol0_(numeric_limits<real>::epsilon())
00049
00050
00051
00052 , tol1_(200 * tol0_)
00053 , tol2_(sqrt(tol0_))
00054
00055 , tolb_(tol0_ * tol2_)
00056 , xthresh_(1000 * tol2_)
00057 , _a(a)
00058 , _f(f <= 1 ? f : 1/f)
00059 , _f1(1 - _f)
00060 , _e2(_f * (2 - _f))
00061 , _ep2(_e2 / Math::sq(_f1))
00062 , _n(_f / ( 2 - _f))
00063 , _b(_a * _f1)
00064 , _c2((Math::sq(_a) + Math::sq(_b) *
00065 (_e2 == 0 ? 1 :
00066 (_e2 > 0 ? Math::atanh(sqrt(_e2)) : atan(sqrt(-_e2))) /
00067 sqrt(abs(_e2))))/2)
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078 , _etol2(0.1 * tol2_ /
00079 sqrt( max(real(0.001), abs(_f)) * min(real(1), 1 - _f/2) / 2 ))
00080 {
00081 if (!(Math::isfinite(_a) && _a > 0))
00082 throw GeographicErr("Major radius is not positive");
00083 if (!(Math::isfinite(_b) && _b > 0))
00084 throw GeographicErr("Minor radius is not positive");
00085 A3coeff();
00086 C3coeff();
00087 C4coeff();
00088 }
00089
00090 const Geodesic& Geodesic::WGS84() {
00091 static const Geodesic wgs84(Constants::WGS84_a(), Constants::WGS84_f());
00092 return wgs84;
00093 }
00094
00095 Math::real Geodesic::SinCosSeries(bool sinp,
00096 real sinx, real cosx,
00097 const real c[], int n) {
00098
00099
00100
00101
00102
00103 c += (n + sinp);
00104 real
00105 ar = 2 * (cosx - sinx) * (cosx + sinx),
00106 y0 = n & 1 ? *--c : 0, y1 = 0;
00107
00108 n /= 2;
00109 while (n--) {
00110
00111 y1 = ar * y0 - y1 + *--c;
00112 y0 = ar * y1 - y0 + *--c;
00113 }
00114 return sinp
00115 ? 2 * sinx * cosx * y0
00116 : cosx * (y0 - y1);
00117 }
00118
00119 GeodesicLine Geodesic::Line(real lat1, real lon1, real azi1, unsigned caps)
00120 const {
00121 return GeodesicLine(*this, lat1, lon1, azi1, caps);
00122 }
00123
00124 Math::real Geodesic::GenDirect(real lat1, real lon1, real azi1,
00125 bool arcmode, real s12_a12, unsigned outmask,
00126 real& lat2, real& lon2, real& azi2,
00127 real& s12, real& m12, real& M12, real& M21,
00128 real& S12) const {
00129 return GeodesicLine(*this, lat1, lon1, azi1,
00130
00131 outmask | (arcmode ? NONE : DISTANCE_IN))
00132 .
00133 GenPosition(arcmode, s12_a12, outmask,
00134 lat2, lon2, azi2, s12, m12, M12, M21, S12);
00135 }
00136
00137 Math::real Geodesic::GenInverse(real lat1, real lon1, real lat2, real lon2,
00138 unsigned outmask,
00139 real& s12, real& azi1, real& azi2,
00140 real& m12, real& M12, real& M21, real& S12)
00141 const {
00142 outmask &= OUT_ALL;
00143
00144
00145
00146 real lon12 = Math::AngDiff(Math::AngNormalize(lon1),
00147 Math::AngNormalize(lon2));
00148
00149 lon12 = AngRound(lon12);
00150
00151 int lonsign = lon12 >= 0 ? 1 : -1;
00152 lon12 *= lonsign;
00153
00154 lat1 = AngRound(lat1);
00155 lat2 = AngRound(lat2);
00156
00157 int swapp = abs(lat1) >= abs(lat2) ? 1 : -1;
00158 if (swapp < 0) {
00159 lonsign *= -1;
00160 swap(lat1, lat2);
00161 }
00162
00163 int latsign = lat1 < 0 ? 1 : -1;
00164 lat1 *= latsign;
00165 lat2 *= latsign;
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x;
00179
00180 phi = lat1 * Math::degree();
00181
00182 sbet1 = _f1 * sin(phi);
00183 cbet1 = lat1 == -90 ? tiny_ : cos(phi);
00184 SinCosNorm(sbet1, cbet1);
00185
00186 phi = lat2 * Math::degree();
00187
00188 sbet2 = _f1 * sin(phi);
00189 cbet2 = abs(lat2) == 90 ? tiny_ : cos(phi);
00190 SinCosNorm(sbet2, cbet2);
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200 if (cbet1 < -sbet1) {
00201 if (cbet2 == cbet1)
00202 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
00203 } else {
00204 if (abs(sbet2) == -sbet1)
00205 cbet2 = cbet1;
00206 }
00207
00208 real
00209 dn1 = sqrt(1 + _ep2 * Math::sq(sbet1)),
00210 dn2 = sqrt(1 + _ep2 * Math::sq(sbet2));
00211
00212 real
00213 lam12 = lon12 * Math::degree(),
00214 slam12 = abs(lon12) == 180 ? 0 : sin(lam12),
00215 clam12 = cos(lam12);
00216
00217
00218 real a12, sig12, calp1, salp1, calp2 = 0, salp2 = 0;
00219
00220 real C1a[nC1_ + 1], C2a[nC2_ + 1], C3a[nC3_];
00221
00222 bool meridian = lat1 == -90 || slam12 == 0;
00223
00224 if (meridian) {
00225
00226
00227
00228
00229 calp1 = clam12; salp1 = slam12;
00230 calp2 = 1; salp2 = 0;
00231
00232 real
00233
00234 ssig1 = sbet1, csig1 = calp1 * cbet1,
00235 ssig2 = sbet2, csig2 = calp2 * cbet2;
00236
00237
00238 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
00239 csig1 * csig2 + ssig1 * ssig2);
00240 {
00241 real dummy;
00242 Lengths(_n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
00243 cbet1, cbet2, s12x, m12x, dummy,
00244 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a);
00245 }
00246
00247
00248
00249
00250
00251
00252
00253 if (sig12 < 1 || m12x >= 0) {
00254 m12x *= _b;
00255 s12x *= _b;
00256 a12 = sig12 / Math::degree();
00257 } else
00258
00259 meridian = false;
00260 }
00261
00262 real omg12 = 0;
00263 if (!meridian &&
00264 sbet1 == 0 &&
00265
00266 (_f <= 0 || lam12 <= Math::pi() - _f * Math::pi())) {
00267
00268
00269 calp1 = calp2 = 0; salp1 = salp2 = 1;
00270 s12x = _a * lam12;
00271 sig12 = omg12 = lam12 / _f1;
00272 m12x = _b * sin(sig12);
00273 if (outmask & GEODESICSCALE)
00274 M12 = M21 = cos(sig12);
00275 a12 = lon12 / _f1;
00276
00277 } else if (!meridian) {
00278
00279
00280
00281
00282
00283 real dnm;
00284 sig12 = InverseStart(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
00285 lam12,
00286 salp1, calp1, salp2, calp2, dnm,
00287 C1a, C2a);
00288
00289 if (sig12 >= 0) {
00290
00291 s12x = sig12 * _b * dnm;
00292 m12x = Math::sq(dnm) * _b * sin(sig12 / dnm);
00293 if (outmask & GEODESICSCALE)
00294 M12 = M21 = cos(sig12 / dnm);
00295 a12 = sig12 / Math::degree();
00296 omg12 = lam12 / (_f1 * dnm);
00297 } else {
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0;
00313 unsigned numit = 0;
00314
00315 real salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
00316 for (bool tripn = false, tripb = false;
00317 numit < maxit2_ || GEOGRAPHICLIB_PANIC;
00318 ++numit) {
00319
00320
00321 real dv;
00322 real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
00323 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
00324 eps, omg12, numit < maxit1_, dv, C1a, C2a, C3a)
00325 - lam12;
00326
00327
00328 if (tripb || !(abs(v) >= (tripn ? 8 : 2) * tol0_)) break;
00329
00330 if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
00331 { salp1b = salp1; calp1b = calp1; }
00332 else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
00333 { salp1a = salp1; calp1a = calp1; }
00334 if (numit < maxit1_ && dv > 0) {
00335 real
00336 dalp1 = -v/dv;
00337 real
00338 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
00339 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
00340 if (nsalp1 > 0 && abs(dalp1) < Math::pi()) {
00341 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
00342 salp1 = nsalp1;
00343 SinCosNorm(salp1, calp1);
00344
00345
00346
00347 tripn = abs(v) <= 16 * tol0_;
00348 continue;
00349 }
00350 }
00351
00352
00353
00354
00355
00356
00357
00358
00359 salp1 = (salp1a + salp1b)/2;
00360 calp1 = (calp1a + calp1b)/2;
00361 SinCosNorm(salp1, calp1);
00362 tripn = false;
00363 tripb = (abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
00364 abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
00365 }
00366 {
00367 real dummy;
00368 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
00369 cbet1, cbet2, s12x, m12x, dummy,
00370 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a);
00371 }
00372 m12x *= _b;
00373 s12x *= _b;
00374 a12 = sig12 / Math::degree();
00375 omg12 = lam12 - omg12;
00376 }
00377 }
00378
00379 if (outmask & DISTANCE)
00380 s12 = 0 + s12x;
00381
00382 if (outmask & REDUCEDLENGTH)
00383 m12 = 0 + m12x;
00384
00385 if (outmask & AREA) {
00386 real
00387
00388 salp0 = salp1 * cbet1,
00389 calp0 = Math::hypot(calp1, salp1 * sbet1);
00390 real alp12;
00391 if (calp0 != 0 && salp0 != 0) {
00392 real
00393
00394 ssig1 = sbet1, csig1 = calp1 * cbet1,
00395 ssig2 = sbet2, csig2 = calp2 * cbet2,
00396 k2 = Math::sq(calp0) * _ep2,
00397 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
00398
00399 A4 = Math::sq(_a) * calp0 * salp0 * _e2;
00400 SinCosNorm(ssig1, csig1);
00401 SinCosNorm(ssig2, csig2);
00402 real C4a[nC4_];
00403 C4f(eps, C4a);
00404 real
00405 B41 = SinCosSeries(false, ssig1, csig1, C4a, nC4_),
00406 B42 = SinCosSeries(false, ssig2, csig2, C4a, nC4_);
00407 S12 = A4 * (B42 - B41);
00408 } else
00409
00410 S12 = 0;
00411
00412 if (!meridian &&
00413 omg12 < real(0.75) * Math::pi() &&
00414 sbet2 - sbet1 < real(1.75)) {
00415
00416
00417
00418 real
00419 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
00420 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
00421 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
00422 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
00423 } else {
00424
00425 real
00426 salp12 = salp2 * calp1 - calp2 * salp1,
00427 calp12 = calp2 * calp1 + salp2 * salp1;
00428
00429
00430
00431
00432 if (salp12 == 0 && calp12 < 0) {
00433 salp12 = tiny_ * calp1;
00434 calp12 = -1;
00435 }
00436 alp12 = atan2(salp12, calp12);
00437 }
00438 S12 += _c2 * alp12;
00439 S12 *= swapp * lonsign * latsign;
00440
00441 S12 += 0;
00442 }
00443
00444
00445 if (swapp < 0) {
00446 swap(salp1, salp2);
00447 swap(calp1, calp2);
00448 if (outmask & GEODESICSCALE)
00449 swap(M12, M21);
00450 }
00451
00452 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
00453 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
00454
00455 if (outmask & AZIMUTH) {
00456
00457 azi1 = 0 - atan2(-salp1, calp1) / Math::degree();
00458 azi2 = 0 - atan2(-salp2, calp2) / Math::degree();
00459 }
00460
00461
00462 return a12;
00463 }
00464
00465 void Geodesic::Lengths(real eps, real sig12,
00466 real ssig1, real csig1, real dn1,
00467 real ssig2, real csig2, real dn2,
00468 real cbet1, real cbet2,
00469 real& s12b, real& m12b, real& m0,
00470 bool scalep, real& M12, real& M21,
00471
00472 real C1a[], real C2a[]) const {
00473
00474
00475 C1f(eps, C1a);
00476 C2f(eps, C2a);
00477 real
00478 A1m1 = A1m1f(eps),
00479 AB1 = (1 + A1m1) * (SinCosSeries(true, ssig2, csig2, C1a, nC1_) -
00480 SinCosSeries(true, ssig1, csig1, C1a, nC1_)),
00481 A2m1 = A2m1f(eps),
00482 AB2 = (1 + A2m1) * (SinCosSeries(true, ssig2, csig2, C2a, nC2_) -
00483 SinCosSeries(true, ssig1, csig1, C2a, nC2_));
00484 m0 = A1m1 - A2m1;
00485 real J12 = m0 * sig12 + (AB1 - AB2);
00486
00487
00488
00489 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12;
00490
00491 s12b = (1 + A1m1) * sig12 + AB1;
00492 if (scalep) {
00493 real csig12 = csig1 * csig2 + ssig1 * ssig2;
00494 real t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
00495 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
00496 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
00497 }
00498 }
00499
00500 Math::real Geodesic::Astroid(real x, real y) {
00501
00502
00503 real k;
00504 real
00505 p = Math::sq(x),
00506 q = Math::sq(y),
00507 r = (p + q - 1) / 6;
00508 if ( !(q == 0 && r <= 0) ) {
00509 real
00510
00511
00512 S = p * q / 4,
00513 r2 = Math::sq(r),
00514 r3 = r * r2,
00515
00516
00517 disc = S * (S + 2 * r3);
00518 real u = r;
00519 if (disc >= 0) {
00520 real T3 = S + r3;
00521
00522
00523
00524 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
00525
00526 real T = Math::cbrt(T3);
00527
00528 u += T + (T ? r2 / T : 0);
00529 } else {
00530
00531 real ang = atan2(sqrt(-disc), -(S + r3));
00532
00533
00534 u += 2 * r * cos(ang / 3);
00535 }
00536 real
00537 v = sqrt(Math::sq(u) + q),
00538
00539 uv = u < 0 ? q / (v - u) : u + v,
00540 w = (uv - q) / (2 * v);
00541
00542
00543 k = uv / (sqrt(uv + Math::sq(w)) + w);
00544 } else {
00545
00546
00547 k = 0;
00548 }
00549 return k;
00550 }
00551
00552 Math::real Geodesic::InverseStart(real sbet1, real cbet1, real dn1,
00553 real sbet2, real cbet2, real dn2,
00554 real lam12,
00555 real& salp1, real& calp1,
00556
00557 real& salp2, real& calp2,
00558
00559 real& dnm,
00560
00561 real C1a[], real C2a[]) const {
00562
00563
00564
00565 real
00566 sig12 = -1,
00567
00568 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
00569 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
00570 #if defined(__GNUC__) && __GNUC__ == 4 && \
00571 (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
00572
00573
00574
00575
00576
00577
00578 real sbet12a;
00579 {
00580 GEOGRAPHICLIB_VOLATILE real xx1 = sbet2 * cbet1;
00581 GEOGRAPHICLIB_VOLATILE real xx2 = cbet2 * sbet1;
00582 sbet12a = xx1 + xx2;
00583 }
00584 #else
00585 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
00586 #endif
00587 bool shortline = cbet12 >= 0 && sbet12 < real(0.5) &&
00588 cbet2 * lam12 < real(0.5);
00589 real omg12 = lam12;
00590 if (shortline) {
00591 real sbetm2 = Math::sq(sbet1 + sbet2);
00592
00593
00594 sbetm2 /= sbetm2 + Math::sq(cbet1 + cbet2);
00595 dnm = sqrt(1 + _ep2 * sbetm2);
00596 omg12 /= _f1 * dnm;
00597 }
00598 real somg12 = sin(omg12), comg12 = cos(omg12);
00599
00600 salp1 = cbet2 * somg12;
00601 calp1 = comg12 >= 0 ?
00602 sbet12 + cbet2 * sbet1 * Math::sq(somg12) / (1 + comg12) :
00603 sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
00604
00605 real
00606 ssig12 = Math::hypot(salp1, calp1),
00607 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
00608
00609 if (shortline && ssig12 < _etol2) {
00610
00611 salp2 = cbet1 * somg12;
00612 calp2 = sbet12 - cbet1 * sbet2 *
00613 (comg12 >= 0 ? Math::sq(somg12) / (1 + comg12) : 1 - comg12);
00614 SinCosNorm(salp2, calp2);
00615
00616 sig12 = atan2(ssig12, csig12);
00617 } else if (abs(_n) > real(0.1) ||
00618 csig12 >= 0 ||
00619 ssig12 >= 6 * abs(_n) * Math::pi() * Math::sq(cbet1)) {
00620
00621 } else {
00622
00623
00624 real y, lamscale, betscale;
00625
00626
00627
00628 GEOGRAPHICLIB_VOLATILE real x;
00629 if (_f >= 0) {
00630
00631 {
00632 real
00633 k2 = Math::sq(sbet1) * _ep2,
00634 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
00635 lamscale = _f * cbet1 * A3f(eps) * Math::pi();
00636 }
00637 betscale = lamscale * cbet1;
00638
00639 x = (lam12 - Math::pi()) / lamscale;
00640 y = sbet12a / betscale;
00641 } else {
00642
00643 real
00644 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
00645 bet12a = atan2(sbet12a, cbet12a);
00646 real m12b, m0, dummy;
00647
00648
00649 Lengths(_n, Math::pi() + bet12a,
00650 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
00651 cbet1, cbet2, dummy, m12b, m0, false,
00652 dummy, dummy, C1a, C2a);
00653 x = -1 + m12b / (cbet1 * cbet2 * m0 * Math::pi());
00654 betscale = x < -real(0.01) ? sbet12a / x :
00655 -_f * Math::sq(cbet1) * Math::pi();
00656 lamscale = betscale / cbet1;
00657 y = (lam12 - Math::pi()) / lamscale;
00658 }
00659
00660 if (y > -tol1_ && x > -1 - xthresh_) {
00661
00662
00663 if (_f >= 0) {
00664 salp1 = min(real(1), -real(x)); calp1 = - sqrt(1 - Math::sq(salp1));
00665 } else {
00666 calp1 = max(real(x > -tol1_ ? 0 : -1), real(x));
00667 salp1 = sqrt(1 - Math::sq(calp1));
00668 }
00669 } else {
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704 real k = Astroid(x, y);
00705 real
00706 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
00707 somg12 = sin(omg12a); comg12 = -cos(omg12a);
00708
00709 salp1 = cbet2 * somg12;
00710 calp1 = sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
00711 }
00712 }
00713 if (salp1 > 0)
00714 SinCosNorm(salp1, calp1);
00715 else {
00716 salp1 = 1; calp1 = 0;
00717 }
00718 return sig12;
00719 }
00720
00721 Math::real Geodesic::Lambda12(real sbet1, real cbet1, real dn1,
00722 real sbet2, real cbet2, real dn2,
00723 real salp1, real calp1,
00724 real& salp2, real& calp2,
00725 real& sig12,
00726 real& ssig1, real& csig1,
00727 real& ssig2, real& csig2,
00728 real& eps, real& domg12,
00729 bool diffp, real& dlam12,
00730
00731 real C1a[], real C2a[], real C3a[]) const {
00732
00733 if (sbet1 == 0 && calp1 == 0)
00734
00735
00736 calp1 = -tiny_;
00737
00738 real
00739
00740 salp0 = salp1 * cbet1,
00741 calp0 = Math::hypot(calp1, salp1 * sbet1);
00742
00743 real somg1, comg1, somg2, comg2, omg12, lam12;
00744
00745
00746 ssig1 = sbet1; somg1 = salp0 * sbet1;
00747 csig1 = comg1 = calp1 * cbet1;
00748 SinCosNorm(ssig1, csig1);
00749
00750
00751
00752
00753
00754
00755 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
00756
00757
00758
00759
00760 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ?
00761 sqrt(Math::sq(calp1 * cbet1) +
00762 (cbet1 < -sbet1 ?
00763 (cbet2 - cbet1) * (cbet1 + cbet2) :
00764 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
00765 abs(calp1);
00766
00767
00768 ssig2 = sbet2; somg2 = salp0 * sbet2;
00769 csig2 = comg2 = calp2 * cbet2;
00770 SinCosNorm(ssig2, csig2);
00771
00772
00773
00774 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
00775 csig1 * csig2 + ssig1 * ssig2);
00776
00777
00778 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, real(0)),
00779 comg1 * comg2 + somg1 * somg2);
00780 real B312, h0;
00781 real k2 = Math::sq(calp0) * _ep2;
00782 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
00783 C3f(eps, C3a);
00784 B312 = (SinCosSeries(true, ssig2, csig2, C3a, nC3_-1) -
00785 SinCosSeries(true, ssig1, csig1, C3a, nC3_-1));
00786 h0 = -_f * A3f(eps);
00787 domg12 = salp0 * h0 * (sig12 + B312);
00788 lam12 = omg12 + domg12;
00789
00790 if (diffp) {
00791 if (calp2 == 0)
00792 dlam12 = - 2 * _f1 * dn1 / sbet1;
00793 else {
00794 real dummy;
00795 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
00796 cbet1, cbet2, dummy, dlam12, dummy,
00797 false, dummy, dummy, C1a, C2a);
00798 dlam12 *= _f1 / (calp2 * cbet2);
00799 }
00800 }
00801
00802 return lam12;
00803 }
00804
00805 Math::real Geodesic::A3f(real eps) const {
00806
00807 real v = 0;
00808 for (int i = nA3x_; i > 0; )
00809 v = eps * v + _A3x[--i];
00810 return v;
00811 }
00812
00813 void Geodesic::C3f(real eps, real c[]) const {
00814
00815
00816 for (int j = nC3x_, k = nC3_ - 1; k > 0; ) {
00817 real t = 0;
00818 for (int i = nC3_ - k; i > 0; --i) {
00819 t = eps * t + _C3x[--j];
00820 }
00821 c[k--] = t;
00822 }
00823
00824 real mult = 1;
00825 for (int k = 1; k < nC3_; ) {
00826 mult *= eps;
00827 c[k++] *= mult;
00828 }
00829 }
00830
00831 void Geodesic::C4f(real eps, real c[]) const {
00832
00833
00834 for (int j = nC4x_, k = nC4_; k > 0; ) {
00835 real t = 0;
00836 for (int i = nC4_ - k + 1; i > 0; --i)
00837 t = eps * t + _C4x[--j];
00838 c[--k] = t;
00839 }
00840
00841 real mult = 1;
00842 for (int k = 1; k < nC4_; ) {
00843 mult *= eps;
00844 c[k++] *= mult;
00845 }
00846 }
00847
00848
00849
00850
00851 Math::real Geodesic::A1m1f(real eps) {
00852 real
00853 eps2 = Math::sq(eps),
00854 t;
00855 switch (nA1_/2) {
00856 case 0:
00857 t = 0;
00858 break;
00859 case 1:
00860 t = eps2/4;
00861 break;
00862 case 2:
00863 t = eps2*(eps2+16)/64;
00864 break;
00865 case 3:
00866 t = eps2*(eps2*(eps2+4)+64)/256;
00867 break;
00868 case 4:
00869 t = eps2*(eps2*(eps2*(25*eps2+64)+256)+4096)/16384;
00870 break;
00871 default:
00872 GEOGRAPHICLIB_STATIC_ASSERT(nA1_ >= 0 && nA1_ <= 8, "Bad value of nA1_");
00873 t = 0;
00874 }
00875 return (t + eps) / (1 - eps);
00876 }
00877
00878
00879 void Geodesic::C1f(real eps, real c[]) {
00880 real
00881 eps2 = Math::sq(eps),
00882 d = eps;
00883 switch (nC1_) {
00884 case 0:
00885 break;
00886 case 1:
00887 c[1] = -d/2;
00888 break;
00889 case 2:
00890 c[1] = -d/2;
00891 d *= eps;
00892 c[2] = -d/16;
00893 break;
00894 case 3:
00895 c[1] = d*(3*eps2-8)/16;
00896 d *= eps;
00897 c[2] = -d/16;
00898 d *= eps;
00899 c[3] = -d/48;
00900 break;
00901 case 4:
00902 c[1] = d*(3*eps2-8)/16;
00903 d *= eps;
00904 c[2] = d*(eps2-2)/32;
00905 d *= eps;
00906 c[3] = -d/48;
00907 d *= eps;
00908 c[4] = -5*d/512;
00909 break;
00910 case 5:
00911 c[1] = d*((6-eps2)*eps2-16)/32;
00912 d *= eps;
00913 c[2] = d*(eps2-2)/32;
00914 d *= eps;
00915 c[3] = d*(9*eps2-16)/768;
00916 d *= eps;
00917 c[4] = -5*d/512;
00918 d *= eps;
00919 c[5] = -7*d/1280;
00920 break;
00921 case 6:
00922 c[1] = d*((6-eps2)*eps2-16)/32;
00923 d *= eps;
00924 c[2] = d*((64-9*eps2)*eps2-128)/2048;
00925 d *= eps;
00926 c[3] = d*(9*eps2-16)/768;
00927 d *= eps;
00928 c[4] = d*(3*eps2-5)/512;
00929 d *= eps;
00930 c[5] = -7*d/1280;
00931 d *= eps;
00932 c[6] = -7*d/2048;
00933 break;
00934 case 7:
00935 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048;
00936 d *= eps;
00937 c[2] = d*((64-9*eps2)*eps2-128)/2048;
00938 d *= eps;
00939 c[3] = d*((72-9*eps2)*eps2-128)/6144;
00940 d *= eps;
00941 c[4] = d*(3*eps2-5)/512;
00942 d *= eps;
00943 c[5] = d*(35*eps2-56)/10240;
00944 d *= eps;
00945 c[6] = -7*d/2048;
00946 d *= eps;
00947 c[7] = -33*d/14336;
00948 break;
00949 case 8:
00950 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048;
00951 d *= eps;
00952 c[2] = d*(eps2*(eps2*(7*eps2-18)+128)-256)/4096;
00953 d *= eps;
00954 c[3] = d*((72-9*eps2)*eps2-128)/6144;
00955 d *= eps;
00956 c[4] = d*((96-11*eps2)*eps2-160)/16384;
00957 d *= eps;
00958 c[5] = d*(35*eps2-56)/10240;
00959 d *= eps;
00960 c[6] = d*(9*eps2-14)/4096;
00961 d *= eps;
00962 c[7] = -33*d/14336;
00963 d *= eps;
00964 c[8] = -429*d/262144;
00965 break;
00966 default:
00967 GEOGRAPHICLIB_STATIC_ASSERT(nC1_ >= 0 && nC1_ <= 8, "Bad value of nC1_");
00968 }
00969 }
00970
00971
00972 void Geodesic::C1pf(real eps, real c[]) {
00973 real
00974 eps2 = Math::sq(eps),
00975 d = eps;
00976 switch (nC1p_) {
00977 case 0:
00978 break;
00979 case 1:
00980 c[1] = d/2;
00981 break;
00982 case 2:
00983 c[1] = d/2;
00984 d *= eps;
00985 c[2] = 5*d/16;
00986 break;
00987 case 3:
00988 c[1] = d*(16-9*eps2)/32;
00989 d *= eps;
00990 c[2] = 5*d/16;
00991 d *= eps;
00992 c[3] = 29*d/96;
00993 break;
00994 case 4:
00995 c[1] = d*(16-9*eps2)/32;
00996 d *= eps;
00997 c[2] = d*(30-37*eps2)/96;
00998 d *= eps;
00999 c[3] = 29*d/96;
01000 d *= eps;
01001 c[4] = 539*d/1536;
01002 break;
01003 case 5:
01004 c[1] = d*(eps2*(205*eps2-432)+768)/1536;
01005 d *= eps;
01006 c[2] = d*(30-37*eps2)/96;
01007 d *= eps;
01008 c[3] = d*(116-225*eps2)/384;
01009 d *= eps;
01010 c[4] = 539*d/1536;
01011 d *= eps;
01012 c[5] = 3467*d/7680;
01013 break;
01014 case 6:
01015 c[1] = d*(eps2*(205*eps2-432)+768)/1536;
01016 d *= eps;
01017 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
01018 d *= eps;
01019 c[3] = d*(116-225*eps2)/384;
01020 d *= eps;
01021 c[4] = d*(2695-7173*eps2)/7680;
01022 d *= eps;
01023 c[5] = 3467*d/7680;
01024 d *= eps;
01025 c[6] = 38081*d/61440;
01026 break;
01027 case 7:
01028 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728;
01029 d *= eps;
01030 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
01031 d *= eps;
01032 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288;
01033 d *= eps;
01034 c[4] = d*(2695-7173*eps2)/7680;
01035 d *= eps;
01036 c[5] = d*(41604-141115*eps2)/92160;
01037 d *= eps;
01038 c[6] = 38081*d/61440;
01039 d *= eps;
01040 c[7] = 459485*d/516096;
01041 break;
01042 case 8:
01043 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728;
01044 d *= eps;
01045 c[2] = d*(eps2*((120150-86171*eps2)*eps2-142080)+115200)/368640;
01046 d *= eps;
01047 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288;
01048 d *= eps;
01049 c[4] = d*(eps2*(1082857*eps2-688608)+258720)/737280;
01050 d *= eps;
01051 c[5] = d*(41604-141115*eps2)/92160;
01052 d *= eps;
01053 c[6] = d*(533134-2200311*eps2)/860160;
01054 d *= eps;
01055 c[7] = 459485*d/516096;
01056 d *= eps;
01057 c[8] = 109167851*d/82575360;
01058 break;
01059 default:
01060 GEOGRAPHICLIB_STATIC_ASSERT(nC1p_ >= 0 && nC1p_ <= 8,
01061 "Bad value of nC1p_");
01062 }
01063 }
01064
01065
01066 Math::real Geodesic::A2m1f(real eps) {
01067 real
01068 eps2 = Math::sq(eps),
01069 t;
01070 switch (nA2_/2) {
01071 case 0:
01072 t = 0;
01073 break;
01074 case 1:
01075 t = eps2/4;
01076 break;
01077 case 2:
01078 t = eps2*(9*eps2+16)/64;
01079 break;
01080 case 3:
01081 t = eps2*(eps2*(25*eps2+36)+64)/256;
01082 break;
01083 case 4:
01084 t = eps2*(eps2*(eps2*(1225*eps2+1600)+2304)+4096)/16384;
01085 break;
01086 default:
01087 GEOGRAPHICLIB_STATIC_ASSERT(nA2_ >= 0 && nA2_ <= 8, "Bad value of nA2_");
01088 t = 0;
01089 }
01090 return t * (1 - eps) - eps;
01091 }
01092
01093
01094 void Geodesic::C2f(real eps, real c[]) {
01095 real
01096 eps2 = Math::sq(eps),
01097 d = eps;
01098 switch (nC2_) {
01099 case 0:
01100 break;
01101 case 1:
01102 c[1] = d/2;
01103 break;
01104 case 2:
01105 c[1] = d/2;
01106 d *= eps;
01107 c[2] = 3*d/16;
01108 break;
01109 case 3:
01110 c[1] = d*(eps2+8)/16;
01111 d *= eps;
01112 c[2] = 3*d/16;
01113 d *= eps;
01114 c[3] = 5*d/48;
01115 break;
01116 case 4:
01117 c[1] = d*(eps2+8)/16;
01118 d *= eps;
01119 c[2] = d*(eps2+6)/32;
01120 d *= eps;
01121 c[3] = 5*d/48;
01122 d *= eps;
01123 c[4] = 35*d/512;
01124 break;
01125 case 5:
01126 c[1] = d*(eps2*(eps2+2)+16)/32;
01127 d *= eps;
01128 c[2] = d*(eps2+6)/32;
01129 d *= eps;
01130 c[3] = d*(15*eps2+80)/768;
01131 d *= eps;
01132 c[4] = 35*d/512;
01133 d *= eps;
01134 c[5] = 63*d/1280;
01135 break;
01136 case 6:
01137 c[1] = d*(eps2*(eps2+2)+16)/32;
01138 d *= eps;
01139 c[2] = d*(eps2*(35*eps2+64)+384)/2048;
01140 d *= eps;
01141 c[3] = d*(15*eps2+80)/768;
01142 d *= eps;
01143 c[4] = d*(7*eps2+35)/512;
01144 d *= eps;
01145 c[5] = 63*d/1280;
01146 d *= eps;
01147 c[6] = 77*d/2048;
01148 break;
01149 case 7:
01150 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048;
01151 d *= eps;
01152 c[2] = d*(eps2*(35*eps2+64)+384)/2048;
01153 d *= eps;
01154 c[3] = d*(eps2*(69*eps2+120)+640)/6144;
01155 d *= eps;
01156 c[4] = d*(7*eps2+35)/512;
01157 d *= eps;
01158 c[5] = d*(105*eps2+504)/10240;
01159 d *= eps;
01160 c[6] = 77*d/2048;
01161 d *= eps;
01162 c[7] = 429*d/14336;
01163 break;
01164 case 8:
01165 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048;
01166 d *= eps;
01167 c[2] = d*(eps2*(eps2*(47*eps2+70)+128)+768)/4096;
01168 d *= eps;
01169 c[3] = d*(eps2*(69*eps2+120)+640)/6144;
01170 d *= eps;
01171 c[4] = d*(eps2*(133*eps2+224)+1120)/16384;
01172 d *= eps;
01173 c[5] = d*(105*eps2+504)/10240;
01174 d *= eps;
01175 c[6] = d*(33*eps2+154)/4096;
01176 d *= eps;
01177 c[7] = 429*d/14336;
01178 d *= eps;
01179 c[8] = 6435*d/262144;
01180 break;
01181 default:
01182 GEOGRAPHICLIB_STATIC_ASSERT(nC2_ >= 0 && nC2_ <= 8, "Bad value of nC2_");
01183 }
01184 }
01185
01186
01187 void Geodesic::A3coeff() {
01188 switch (nA3_) {
01189 case 0:
01190 break;
01191 case 1:
01192 _A3x[0] = 1;
01193 break;
01194 case 2:
01195 _A3x[0] = 1;
01196 _A3x[1] = -1/real(2);
01197 break;
01198 case 3:
01199 _A3x[0] = 1;
01200 _A3x[1] = (_n-1)/2;
01201 _A3x[2] = -1/real(4);
01202 break;
01203 case 4:
01204 _A3x[0] = 1;
01205 _A3x[1] = (_n-1)/2;
01206 _A3x[2] = (-_n-2)/8;
01207 _A3x[3] = -1/real(16);
01208 break;
01209 case 5:
01210 _A3x[0] = 1;
01211 _A3x[1] = (_n-1)/2;
01212 _A3x[2] = (_n*(3*_n-1)-2)/8;
01213 _A3x[3] = (-3*_n-1)/16;
01214 _A3x[4] = -3/real(64);
01215 break;
01216 case 6:
01217 _A3x[0] = 1;
01218 _A3x[1] = (_n-1)/2;
01219 _A3x[2] = (_n*(3*_n-1)-2)/8;
01220 _A3x[3] = ((-_n-3)*_n-1)/16;
01221 _A3x[4] = (-2*_n-3)/64;
01222 _A3x[5] = -3/real(128);
01223 break;
01224 case 7:
01225 _A3x[0] = 1;
01226 _A3x[1] = (_n-1)/2;
01227 _A3x[2] = (_n*(3*_n-1)-2)/8;
01228 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16;
01229 _A3x[4] = ((-10*_n-2)*_n-3)/64;
01230 _A3x[5] = (-5*_n-3)/128;
01231 _A3x[6] = -5/real(256);
01232 break;
01233 case 8:
01234 _A3x[0] = 1;
01235 _A3x[1] = (_n-1)/2;
01236 _A3x[2] = (_n*(3*_n-1)-2)/8;
01237 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16;
01238 _A3x[4] = (_n*((-5*_n-20)*_n-4)-6)/128;
01239 _A3x[5] = ((-5*_n-10)*_n-6)/256;
01240 _A3x[6] = (-15*_n-20)/1024;
01241 _A3x[7] = -25/real(2048);
01242 break;
01243 default:
01244 GEOGRAPHICLIB_STATIC_ASSERT(nA3_ >= 0 && nA3_ <= 8, "Bad value of nA3_");
01245 }
01246 }
01247
01248
01249 void Geodesic::C3coeff() {
01250 switch (nC3_) {
01251 case 0:
01252 break;
01253 case 1:
01254 break;
01255 case 2:
01256 _C3x[0] = 1/real(4);
01257 break;
01258 case 3:
01259 _C3x[0] = (1-_n)/4;
01260 _C3x[1] = 1/real(8);
01261 _C3x[2] = 1/real(16);
01262 break;
01263 case 4:
01264 _C3x[0] = (1-_n)/4;
01265 _C3x[1] = 1/real(8);
01266 _C3x[2] = 3/real(64);
01267 _C3x[3] = (2-3*_n)/32;
01268 _C3x[4] = 3/real(64);
01269 _C3x[5] = 5/real(192);
01270 break;
01271 case 5:
01272 _C3x[0] = (1-_n)/4;
01273 _C3x[1] = (1-_n*_n)/8;
01274 _C3x[2] = (3*_n+3)/64;
01275 _C3x[3] = 5/real(128);
01276 _C3x[4] = ((_n-3)*_n+2)/32;
01277 _C3x[5] = (3-2*_n)/64;
01278 _C3x[6] = 3/real(128);
01279 _C3x[7] = (5-9*_n)/192;
01280 _C3x[8] = 3/real(128);
01281 _C3x[9] = 7/real(512);
01282 break;
01283 case 6:
01284 _C3x[0] = (1-_n)/4;
01285 _C3x[1] = (1-_n*_n)/8;
01286 _C3x[2] = ((3-_n)*_n+3)/64;
01287 _C3x[3] = (2*_n+5)/128;
01288 _C3x[4] = 3/real(128);
01289 _C3x[5] = ((_n-3)*_n+2)/32;
01290 _C3x[6] = ((-3*_n-2)*_n+3)/64;
01291 _C3x[7] = (_n+3)/128;
01292 _C3x[8] = 5/real(256);
01293 _C3x[9] = (_n*(5*_n-9)+5)/192;
01294 _C3x[10] = (9-10*_n)/384;
01295 _C3x[11] = 7/real(512);
01296 _C3x[12] = (7-14*_n)/512;
01297 _C3x[13] = 7/real(512);
01298 _C3x[14] = 21/real(2560);
01299 break;
01300 case 7:
01301 _C3x[0] = (1-_n)/4;
01302 _C3x[1] = (1-_n*_n)/8;
01303 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64;
01304 _C3x[3] = (_n*(2*_n+2)+5)/128;
01305 _C3x[4] = (11*_n+12)/512;
01306 _C3x[5] = 21/real(1024);
01307 _C3x[6] = ((_n-3)*_n+2)/32;
01308 _C3x[7] = (_n*(_n*(2*_n-3)-2)+3)/64;
01309 _C3x[8] = ((2-9*_n)*_n+6)/256;
01310 _C3x[9] = (_n+5)/256;
01311 _C3x[10] = 27/real(2048);
01312 _C3x[11] = (_n*((5-_n)*_n-9)+5)/192;
01313 _C3x[12] = ((-6*_n-10)*_n+9)/384;
01314 _C3x[13] = (21-4*_n)/1536;
01315 _C3x[14] = 3/real(256);
01316 _C3x[15] = (_n*(10*_n-14)+7)/512;
01317 _C3x[16] = (7-10*_n)/512;
01318 _C3x[17] = 9/real(1024);
01319 _C3x[18] = (21-45*_n)/2560;
01320 _C3x[19] = 9/real(1024);
01321 _C3x[20] = 11/real(2048);
01322 break;
01323 case 8:
01324 _C3x[0] = (1-_n)/4;
01325 _C3x[1] = (1-_n*_n)/8;
01326 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64;
01327 _C3x[3] = (_n*((2-2*_n)*_n+2)+5)/128;
01328 _C3x[4] = (_n*(3*_n+11)+12)/512;
01329 _C3x[5] = (10*_n+21)/1024;
01330 _C3x[6] = 243/real(16384);
01331 _C3x[7] = ((_n-3)*_n+2)/32;
01332 _C3x[8] = (_n*(_n*(2*_n-3)-2)+3)/64;
01333 _C3x[9] = (_n*((-6*_n-9)*_n+2)+6)/256;
01334 _C3x[10] = ((1-2*_n)*_n+5)/256;
01335 _C3x[11] = (69*_n+108)/8192;
01336 _C3x[12] = 187/real(16384);
01337 _C3x[13] = (_n*((5-_n)*_n-9)+5)/192;
01338 _C3x[14] = (_n*(_n*(10*_n-6)-10)+9)/384;
01339 _C3x[15] = ((-77*_n-8)*_n+42)/3072;
01340 _C3x[16] = (12-_n)/1024;
01341 _C3x[17] = 139/real(16384);
01342 _C3x[18] = (_n*((20-7*_n)*_n-28)+14)/1024;
01343 _C3x[19] = ((-7*_n-40)*_n+28)/2048;
01344 _C3x[20] = (72-43*_n)/8192;
01345 _C3x[21] = 127/real(16384);
01346 _C3x[22] = (_n*(75*_n-90)+42)/5120;
01347 _C3x[23] = (9-15*_n)/1024;
01348 _C3x[24] = 99/real(16384);
01349 _C3x[25] = (44-99*_n)/8192;
01350 _C3x[26] = 99/real(16384);
01351 _C3x[27] = 429/real(114688);
01352 break;
01353 default:
01354 GEOGRAPHICLIB_STATIC_ASSERT(nC3_ >= 0 && nC3_ <= 8, "Bad value of nC3_");
01355 }
01356 }
01357
01358
01359
01360
01361 void Geodesic::C4coeff() {
01362 switch (nC4_) {
01363 case 0:
01364 break;
01365 case 1:
01366 _C4x[0] = 2/real(3);
01367 break;
01368 case 2:
01369 _C4x[0] = (10-4*_n)/15;
01370 _C4x[1] = -1/real(5);
01371 _C4x[2] = 1/real(45);
01372 break;
01373 case 3:
01374 _C4x[0] = (_n*(8*_n-28)+70)/105;
01375 _C4x[1] = (16*_n-7)/35;
01376 _C4x[2] = -2/real(105);
01377 _C4x[3] = (7-16*_n)/315;
01378 _C4x[4] = -2/real(105);
01379 _C4x[5] = 4/real(525);
01380 break;
01381 case 4:
01382 _C4x[0] = (_n*(_n*(4*_n+24)-84)+210)/315;
01383 _C4x[1] = ((48-32*_n)*_n-21)/105;
01384 _C4x[2] = (-32*_n-6)/315;
01385 _C4x[3] = 11/real(315);
01386 _C4x[4] = (_n*(32*_n-48)+21)/945;
01387 _C4x[5] = (64*_n-18)/945;
01388 _C4x[6] = -1/real(105);
01389 _C4x[7] = (12-32*_n)/1575;
01390 _C4x[8] = -8/real(1575);
01391 _C4x[9] = 8/real(2205);
01392 break;
01393 case 5:
01394 _C4x[0] = (_n*(_n*(_n*(16*_n+44)+264)-924)+2310)/3465;
01395 _C4x[1] = (_n*(_n*(48*_n-352)+528)-231)/1155;
01396 _C4x[2] = (_n*(1088*_n-352)-66)/3465;
01397 _C4x[3] = (121-368*_n)/3465;
01398 _C4x[4] = 4/real(1155);
01399 _C4x[5] = (_n*((352-48*_n)*_n-528)+231)/10395;
01400 _C4x[6] = ((704-896*_n)*_n-198)/10395;
01401 _C4x[7] = (80*_n-99)/10395;
01402 _C4x[8] = 4/real(1155);
01403 _C4x[9] = (_n*(320*_n-352)+132)/17325;
01404 _C4x[10] = (384*_n-88)/17325;
01405 _C4x[11] = -8/real(1925);
01406 _C4x[12] = (88-256*_n)/24255;
01407 _C4x[13] = -16/real(8085);
01408 _C4x[14] = 64/real(31185);
01409 break;
01410 case 6:
01411 _C4x[0] = (_n*(_n*(_n*(_n*(100*_n+208)+572)+3432)-12012)+30030)/45045;
01412 _C4x[1] = (_n*(_n*(_n*(64*_n+624)-4576)+6864)-3003)/15015;
01413 _C4x[2] = (_n*((14144-10656*_n)*_n-4576)-858)/45045;
01414 _C4x[3] = ((-224*_n-4784)*_n+1573)/45045;
01415 _C4x[4] = (1088*_n+156)/45045;
01416 _C4x[5] = 97/real(15015);
01417 _C4x[6] = (_n*(_n*((-64*_n-624)*_n+4576)-6864)+3003)/135135;
01418 _C4x[7] = (_n*(_n*(5952*_n-11648)+9152)-2574)/135135;
01419 _C4x[8] = (_n*(5792*_n+1040)-1287)/135135;
01420 _C4x[9] = (468-2944*_n)/135135;
01421 _C4x[10] = 1/real(9009);
01422 _C4x[11] = (_n*((4160-1440*_n)*_n-4576)+1716)/225225;
01423 _C4x[12] = ((4992-8448*_n)*_n-1144)/225225;
01424 _C4x[13] = (1856*_n-936)/225225;
01425 _C4x[14] = 8/real(10725);
01426 _C4x[15] = (_n*(3584*_n-3328)+1144)/315315;
01427 _C4x[16] = (1024*_n-208)/105105;
01428 _C4x[17] = -136/real(63063);
01429 _C4x[18] = (832-2560*_n)/405405;
01430 _C4x[19] = -128/real(135135);
01431 _C4x[20] = 128/real(99099);
01432 break;
01433 case 7:
01434 _C4x[0] = (_n*(_n*(_n*(_n*(_n*(56*_n+100)+208)+572)+3432)-12012)+30030)/
01435 45045;
01436 _C4x[1] = (_n*(_n*(_n*(_n*(16*_n+64)+624)-4576)+6864)-3003)/15015;
01437 _C4x[2] = (_n*(_n*(_n*(1664*_n-10656)+14144)-4576)-858)/45045;
01438 _C4x[3] = (_n*(_n*(10736*_n-224)-4784)+1573)/45045;
01439 _C4x[4] = ((1088-4480*_n)*_n+156)/45045;
01440 _C4x[5] = (291-464*_n)/45045;
01441 _C4x[6] = 10/real(9009);
01442 _C4x[7] = (_n*(_n*(_n*((-16*_n-64)*_n-624)+4576)-6864)+3003)/135135;
01443 _C4x[8] = (_n*(_n*((5952-768*_n)*_n-11648)+9152)-2574)/135135;
01444 _C4x[9] = (_n*((5792-10704*_n)*_n+1040)-1287)/135135;
01445 _C4x[10] = (_n*(3840*_n-2944)+468)/135135;
01446 _C4x[11] = (112*_n+15)/135135;
01447 _C4x[12] = 10/real(9009);
01448 _C4x[13] = (_n*(_n*(_n*(128*_n-1440)+4160)-4576)+1716)/225225;
01449 _C4x[14] = (_n*(_n*(6784*_n-8448)+4992)-1144)/225225;
01450 _C4x[15] = (_n*(1664*_n+1856)-936)/225225;
01451 _C4x[16] = (168-1664*_n)/225225;
01452 _C4x[17] = -4/real(25025);
01453 _C4x[18] = (_n*((3584-1792*_n)*_n-3328)+1144)/315315;
01454 _C4x[19] = ((1024-2048*_n)*_n-208)/105105;
01455 _C4x[20] = (1792*_n-680)/315315;
01456 _C4x[21] = 64/real(315315);
01457 _C4x[22] = (_n*(3072*_n-2560)+832)/405405;
01458 _C4x[23] = (2048*_n-384)/405405;
01459 _C4x[24] = -512/real(405405);
01460 _C4x[25] = (640-2048*_n)/495495;
01461 _C4x[26] = -256/real(495495);
01462 _C4x[27] = 512/real(585585);
01463 break;
01464 case 8:
01465 _C4x[0] = (_n*(_n*(_n*(_n*(_n*(_n*(588*_n+952)+1700)+3536)+9724)+58344)-
01466 204204)+510510)/765765;
01467 _C4x[1] = (_n*(_n*(_n*(_n*(_n*(96*_n+272)+1088)+10608)-77792)+116688)-
01468 51051)/255255;
01469 _C4x[2] = (_n*(_n*(_n*(_n*(3232*_n+28288)-181152)+240448)-77792)-14586)/
01470 765765;
01471 _C4x[3] = (_n*(_n*((182512-154048*_n)*_n-3808)-81328)+26741)/765765;
01472 _C4x[4] = (_n*(_n*(12480*_n-76160)+18496)+2652)/765765;
01473 _C4x[5] = (_n*(20960*_n-7888)+4947)/765765;
01474 _C4x[6] = (4192*_n+850)/765765;
01475 _C4x[7] = 193/real(85085);
01476 _C4x[8] = (_n*(_n*(_n*(_n*((-96*_n-272)*_n-1088)-10608)+77792)-116688)+
01477 51051)/2297295;
01478 _C4x[9] = (_n*(_n*(_n*((-1344*_n-13056)*_n+101184)-198016)+155584)-43758)/
01479 2297295;
01480 _C4x[10] = (_n*(_n*(_n*(103744*_n-181968)+98464)+17680)-21879)/2297295;
01481 _C4x[11] = (_n*(_n*(52608*_n+65280)-50048)+7956)/2297295;
01482 _C4x[12] = ((1904-39840*_n)*_n+255)/2297295;
01483 _C4x[13] = (510-1472*_n)/459459;
01484 _C4x[14] = 349/real(2297295);
01485 _C4x[15] = (_n*(_n*(_n*(_n*(160*_n+2176)-24480)+70720)-77792)+29172)/
01486 3828825;
01487 _C4x[16] = (_n*(_n*((115328-41472*_n)*_n-143616)+84864)-19448)/3828825;
01488 _C4x[17] = (_n*((28288-126528*_n)*_n+31552)-15912)/3828825;
01489 _C4x[18] = (_n*(64256*_n-28288)+2856)/3828825;
01490 _C4x[19] = (-928*_n-612)/3828825;
01491 _C4x[20] = 464/real(1276275);
01492 _C4x[21] = (_n*(_n*(_n*(7168*_n-30464)+60928)-56576)+19448)/5360355;
01493 _C4x[22] = (_n*(_n*(35840*_n-34816)+17408)-3536)/1786785;
01494 _C4x[23] = ((30464-2560*_n)*_n-11560)/5360355;
01495 _C4x[24] = (1088-16384*_n)/5360355;
01496 _C4x[25] = -16/real(97461);
01497 _C4x[26] = (_n*((52224-32256*_n)*_n-43520)+14144)/6891885;
01498 _C4x[27] = ((34816-77824*_n)*_n-6528)/6891885;
01499 _C4x[28] = (26624*_n-8704)/6891885;
01500 _C4x[29] = 128/real(2297295);
01501 _C4x[30] = (_n*(45056*_n-34816)+10880)/8423415;
01502 _C4x[31] = (24576*_n-4352)/8423415;
01503 _C4x[32] = -6784/real(8423415);
01504 _C4x[33] = (8704-28672*_n)/9954945;
01505 _C4x[34] = -1024/real(3318315);
01506 _C4x[35] = 1024/real(1640925);
01507 break;
01508 default:
01509 GEOGRAPHICLIB_STATIC_ASSERT(nC4_ >= 0 && nC4_ <= 8, "Bad value of nC4_");
01510 }
01511 }
01512
01513 }