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 #define GEOGRAPHICLIB_GEODESIC_CPP "$Id: Geodesic.cpp 6937 2011-02-01 20:17:13Z karney $"
00033
00034 RCSID_DECL(GEOGRAPHICLIB_GEODESIC_CPP)
00035 RCSID_DECL(GEOGRAPHICLIB_GEODESIC_HPP)
00036
00037 namespace GeographicLib {
00038
00039 using namespace std;
00040
00041
00042
00043
00044 const Math::real Geodesic::eps2 = sqrt(numeric_limits<real>::min());
00045 const Math::real Geodesic::tol0 = numeric_limits<real>::epsilon();
00046 const Math::real Geodesic::tol1 = 100 * tol0;
00047 const Math::real Geodesic::tol2 = sqrt(numeric_limits<real>::epsilon());
00048 const Math::real Geodesic::xthresh = 1000 * tol2;
00049
00050 Geodesic::Geodesic(real a, real r)
00051 : _a(a)
00052 , _r(r)
00053 , _f(_r != 0 ? 1 / _r : 0)
00054 , _f1(1 - _f)
00055 , _e2(_f * (2 - _f))
00056 , _ep2(_e2 / sq(_f1))
00057 , _n(_f / ( 2 - _f))
00058 , _b(_a * _f1)
00059 , _c2((sq(_a) + sq(_b) *
00060 (_e2 == 0 ? 1 :
00061 (_e2 > 0 ? Math::atanh(sqrt(_e2)) : atan(sqrt(-_e2))) /
00062 sqrt(abs(_e2))))/2)
00063 , _etol2(tol2 / max(real(0.1), sqrt(abs(_e2))))
00064 {
00065 if (!(_a > 0))
00066 throw GeographicErr("Major radius is not positive");
00067 if (!(_f < 1))
00068 throw GeographicErr("Minor radius is not positive");
00069 A3coeff();
00070 C3coeff();
00071 C4coeff();
00072 }
00073
00074 const Geodesic Geodesic::WGS84(Constants::WGS84_a<real>(),
00075 Constants::WGS84_r<real>());
00076
00077 Math::real Geodesic::SinCosSeries(bool sinp,
00078 real sinx, real cosx,
00079 const real c[], int n) throw() {
00080
00081
00082
00083
00084
00085 c += (n + sinp);
00086 real
00087 ar = 2 * (cosx - sinx) * (cosx + sinx),
00088 y0 = n & 1 ? *--c : 0, y1 = 0;
00089
00090 n /= 2;
00091 while (n--) {
00092
00093 y1 = ar * y0 - y1 + *--c;
00094 y0 = ar * y1 - y0 + *--c;
00095 }
00096 return sinp
00097 ? 2 * sinx * cosx * y0
00098 : cosx * (y0 - y1);
00099 }
00100
00101 GeodesicLine Geodesic::Line(real lat1, real lon1, real azi1, unsigned caps)
00102 const throw() {
00103 return GeodesicLine(*this, lat1, lon1, azi1, caps);
00104 }
00105
00106 Math::real Geodesic::GenDirect(real lat1, real lon1, real azi1,
00107 bool arcmode, real s12_a12, unsigned outmask,
00108 real& lat2, real& lon2, real& azi2,
00109 real& s12, real& m12, real& M12, real& M21,
00110 real& S12) const throw() {
00111 return
00112 GeodesicLine(*this, lat1, lon1, azi1,
00113
00114 outmask | (arcmode ? NONE : DISTANCE_IN))
00115 .
00116 GenPosition(arcmode, s12_a12, outmask,
00117 lat2, lon2, azi2, s12, m12, M12, M21, S12);
00118 }
00119
00120 Math::real Geodesic::GenInverse(real lat1, real lon1, real lat2, real lon2,
00121 unsigned outmask,
00122 real& s12, real& azi1, real& azi2,
00123 real& m12, real& M12, real& M21, real& S12)
00124 const throw() {
00125 outmask &= OUT_ALL;
00126 lon1 = AngNormalize(lon1);
00127 real lon12 = AngNormalize(AngNormalize(lon2) - lon1);
00128
00129
00130 lon12 = AngRound(lon12);
00131
00132 int lonsign = lon12 >= 0 ? 1 : -1;
00133 lon12 *= lonsign;
00134 if (lon12 == 180)
00135 lonsign = 1;
00136
00137 lat1 = AngRound(lat1);
00138 lat2 = AngRound(lat2);
00139
00140 int swapp = abs(lat1) >= abs(lat2) ? 1 : -1;
00141 if (swapp < 0) {
00142 lonsign *= -1;
00143 swap(lat1, lat2);
00144 }
00145
00146 int latsign = lat1 < 0 ? 1 : -1;
00147 lat1 *= latsign;
00148 lat2 *= latsign;
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161 real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x;
00162
00163 phi = lat1 * Math::degree<real>();
00164
00165 sbet1 = _f1 * sin(phi);
00166 cbet1 = lat1 == -90 ? eps2 : cos(phi);
00167 SinCosNorm(sbet1, cbet1);
00168
00169 phi = lat2 * Math::degree<real>();
00170
00171 sbet2 = _f1 * sin(phi);
00172 cbet2 = abs(lat2) == 90 ? eps2 : cos(phi);
00173 SinCosNorm(sbet2, cbet2);
00174
00175 real
00176 lam12 = lon12 * Math::degree<real>(),
00177 slam12 = lon12 == 180 ? 0 : sin(lam12),
00178 clam12 = cos(lam12);
00179
00180 real a12, sig12, calp1, salp1, calp2, salp2;
00181
00182 real C1a[nC1 + 1], C2a[nC2 + 1], C3a[nC3];
00183
00184 bool meridian = lat1 == -90 || slam12 == 0;
00185
00186 if (meridian) {
00187
00188
00189
00190
00191 calp1 = clam12; salp1 = slam12;
00192 calp2 = 1; salp2 = 0;
00193
00194 real
00195
00196 ssig1 = sbet1, csig1 = calp1 * cbet1,
00197 ssig2 = sbet2, csig2 = calp2 * cbet2;
00198
00199
00200 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
00201 csig1 * csig2 + ssig1 * ssig2);
00202 {
00203 real dummy;
00204 Lengths(_n, sig12, ssig1, csig1, ssig2, csig2,
00205 cbet1, cbet2, s12x, m12x, dummy,
00206 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a);
00207 }
00208
00209
00210
00211
00212
00213
00214
00215 if (sig12 < 1 || m12x >= 0) {
00216 m12x *= _a;
00217 s12x *= _b;
00218 a12 = sig12 / Math::degree<real>();
00219 } else
00220
00221 meridian = false;
00222 }
00223
00224 real omg12;
00225 if (!meridian &&
00226 sbet1 == 0 &&
00227
00228 (_f <= 0 || lam12 <= Math::pi<real>() - _f * Math::pi<real>())) {
00229
00230
00231 calp1 = calp2 = 0; salp1 = salp2 = 1;
00232 s12x = _a * lam12;
00233 m12x = _b * sin(lam12 / _f1);
00234 if (outmask & GEODESICSCALE)
00235 M12 = M21 = cos(lam12 / _f1);
00236 a12 = lon12 / _f1;
00237 sig12 = lam12 / _f1;
00238
00239 } else if (!meridian) {
00240
00241
00242
00243
00244
00245 sig12 = InverseStart(sbet1, cbet1, sbet2, cbet2,
00246 lam12,
00247 salp1, calp1, salp2, calp2,
00248 C1a, C2a);
00249
00250 if (sig12 >= 0) {
00251
00252 real w1 = sqrt(1 - _e2 * sq(cbet1));
00253 s12x = sig12 * _a * w1;
00254 m12x = sq(w1) * _a / _f1 * sin(sig12 * _f1 / w1);
00255 if (outmask & GEODESICSCALE)
00256 M12 = M21 = cos(sig12 * _f1 / w1);
00257 a12 = sig12 / Math::degree<real>();
00258 omg12 = lam12 / w1;
00259 } else {
00260
00261
00262 real ssig1, csig1, ssig2, csig2, eps;
00263 real ov = 0;
00264 unsigned numit = 0;
00265 for (unsigned trip = 0; numit < maxit; ++numit) {
00266 real dv;
00267 real v = Lambda12(sbet1, cbet1, sbet2, cbet2, salp1, calp1,
00268 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
00269 eps, omg12, trip < 1, dv, C1a, C2a, C3a) - lam12;
00270
00271 if (!(abs(v) > eps2) || !(trip < 1)) {
00272 if (!(abs(v) <= max(tol1, ov)))
00273 numit = maxit;
00274 break;
00275 }
00276 real
00277 dalp1 = -v/dv;
00278 real
00279 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
00280 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
00281 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
00282 salp1 = max(real(0), nsalp1);
00283 SinCosNorm(salp1, calp1);
00284
00285
00286
00287
00288
00289
00290 if (!(abs(v) >= tol1 && sq(v) >= ov * tol0)) ++trip;
00291 ov = abs(v);
00292 }
00293
00294 if (numit >= maxit) {
00295
00296 if (outmask & DISTANCE)
00297 s12 = Math::NaN();
00298 if (outmask & AZIMUTH)
00299 azi1 = azi2 = Math::NaN();
00300 if (outmask & REDUCEDLENGTH)
00301 m12 = Math::NaN();
00302 if (outmask & GEODESICSCALE)
00303 M12 = M21 = Math::NaN();
00304 if (outmask & AREA)
00305 S12 = Math::NaN();
00306 return Math::NaN();
00307 }
00308
00309 {
00310 real dummy;
00311 Lengths(eps, sig12, ssig1, csig1, ssig2, csig2,
00312 cbet1, cbet2, s12x, m12x, dummy,
00313 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a);
00314 }
00315 m12x *= _a;
00316 s12x *= _b;
00317 a12 = sig12 / Math::degree<real>();
00318 omg12 = lam12 - omg12;
00319 }
00320 }
00321
00322 if (outmask & DISTANCE)
00323 s12 = s12x;
00324
00325 if (outmask & REDUCEDLENGTH)
00326 m12 = m12x;
00327
00328 if (outmask & AREA) {
00329 real
00330
00331 salp0 = salp1 * cbet1,
00332 calp0 = Math::hypot(calp1, salp1 * sbet1);
00333 real alp12;
00334 if (calp0 != 0 && salp0 != 0) {
00335 real
00336
00337 ssig1 = sbet1, csig1 = calp1 * cbet1,
00338 ssig2 = sbet2, csig2 = calp2 * cbet2,
00339 k2 = sq(calp0) * _ep2,
00340
00341 A4 = sq(_a) * calp0 * salp0 * _e2;
00342 SinCosNorm(ssig1, csig1);
00343 SinCosNorm(ssig2, csig2);
00344 real C4a[nC4];
00345 C4f(k2, C4a);
00346 real
00347 B41 = Geodesic::SinCosSeries(false, ssig1, csig1, C4a, nC4),
00348 B42 = Geodesic::SinCosSeries(false, ssig2, csig2, C4a, nC4);
00349 S12 = A4 * (B42 - B41);
00350 } else
00351
00352 S12 = 0;
00353 if (!meridian &&
00354 omg12 < real(0.75) * Math::pi<real>() &&
00355 sbet2 - sbet1 < real(1.75)) {
00356
00357
00358
00359 real
00360 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
00361 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
00362 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
00363 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
00364 } else {
00365
00366 real
00367 salp12 = salp2 * calp1 - calp2 * salp1,
00368 calp12 = calp2 * calp1 + salp2 * salp1;
00369
00370
00371
00372
00373 if (salp12 == 0 && calp12 < 0) {
00374 salp12 = Geodesic::eps2 * calp1;
00375 calp12 = -1;
00376 }
00377 alp12 = atan2(salp12, calp12);
00378 }
00379 S12 += _c2 * alp12;
00380 S12 *= swapp * lonsign * latsign;
00381
00382 S12 += 0;
00383 }
00384
00385
00386 if (swapp < 0) {
00387 swap(salp1, salp2);
00388 swap(calp1, calp2);
00389 if (outmask & GEODESICSCALE)
00390 swap(M12, M21);
00391 }
00392
00393 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
00394 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
00395
00396 if (outmask & AZIMUTH) {
00397
00398 azi1 = 0 - atan2(-salp1, calp1) / Math::degree<real>();
00399 azi2 = 0 - atan2(-salp2, calp2) / Math::degree<real>();
00400 }
00401
00402
00403 return a12;
00404 }
00405
00406 void Geodesic::Lengths(real eps, real sig12,
00407 real ssig1, real csig1, real ssig2, real csig2,
00408 real cbet1, real cbet2,
00409 real& s12b, real& m12a, real& m0,
00410 bool scalep, real& M12, real& M21,
00411
00412 real C1a[], real C2a[]) const throw() {
00413
00414
00415 C1f(eps, C1a);
00416 C2f(eps, C2a);
00417 real
00418 A1m1 = A1m1f(eps),
00419 AB1 = (1 + A1m1) * (SinCosSeries(true, ssig2, csig2, C1a, nC1) -
00420 SinCosSeries(true, ssig1, csig1, C1a, nC1)),
00421 A2m1 = A2m1f(eps),
00422 AB2 = (1 + A2m1) * (SinCosSeries(true, ssig2, csig2, C2a, nC2) -
00423 SinCosSeries(true, ssig1, csig1, C2a, nC2)),
00424 cbet1sq = sq(cbet1),
00425 cbet2sq = sq(cbet2),
00426 w1 = sqrt(1 - _e2 * cbet1sq),
00427 w2 = sqrt(1 - _e2 * cbet2sq),
00428
00429 m0x = A1m1 - A2m1,
00430 J12 = m0x * sig12 + (AB1 - AB2);
00431 m0 = m0x;
00432
00433
00434
00435 m12a = (w2 * (csig1 * ssig2) - w1 * (ssig1 * csig2))
00436 - _f1 * csig1 * csig2 * J12;
00437
00438 s12b = (1 + A1m1) * sig12 + AB1;
00439 if (scalep) {
00440 real csig12 = csig1 * csig2 + ssig1 * ssig2;
00441 J12 *= _f1;
00442 M12 = csig12 + (_e2 * (cbet1sq - cbet2sq) * ssig2 / (w1 + w2)
00443 - csig2 * J12) * ssig1 / w1;
00444 M21 = csig12 - (_e2 * (cbet1sq - cbet2sq) * ssig1 / (w1 + w2)
00445 - csig1 * J12) * ssig2 / w2;
00446 }
00447 }
00448
00449 Math::real Geodesic::Astroid(real x, real y) throw() {
00450
00451
00452 real k;
00453 real
00454 p = sq(x),
00455 q = sq(y),
00456 r = (p + q - 1) / 6;
00457 if ( !(q == 0 && r <= 0) ) {
00458 real
00459
00460
00461 S = p * q / 4,
00462 r2 = sq(r),
00463 r3 = r * r2,
00464
00465
00466 disc = S * (S + 2 * r3);
00467 real u = r;
00468 if (disc >= 0) {
00469 real T3 = S + r3;
00470
00471
00472
00473 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
00474
00475 real T = Math::cbrt(T3);
00476
00477 u += T + (T != 0 ? r2 / T : 0);
00478 } else {
00479
00480 real ang = atan2(sqrt(-disc), -(S + r3));
00481
00482
00483 u += 2 * r * cos(ang / 3);
00484 }
00485 real
00486 v = sqrt(sq(u) + q),
00487
00488 uv = u < 0 ? q / (v - u) : u + v,
00489 w = (uv - q) / (2 * v);
00490
00491
00492 k = uv / (sqrt(uv + sq(w)) + w);
00493 } else {
00494
00495
00496 k = 0;
00497 }
00498 return k;
00499 }
00500
00501 Math::real Geodesic::InverseStart(real sbet1, real cbet1,
00502 real sbet2, real cbet2,
00503 real lam12,
00504 real& salp1, real& calp1,
00505
00506 real& salp2, real& calp2,
00507
00508 real C1a[], real C2a[]) const throw() {
00509
00510
00511
00512 real
00513 sig12 = -1,
00514
00515 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
00516 cbet12 = cbet2 * cbet1 + sbet2 * sbet1,
00517 sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
00518
00519 bool shortline = cbet12 >= 0 && sbet12 < real(0.5) &&
00520 lam12 <= Math::pi<real>() / 6;
00521 real
00522 omg12 = shortline ? lam12 / sqrt(1 - _e2 * sq(cbet1)) : lam12,
00523 somg12 = sin(omg12), comg12 = cos(omg12);
00524
00525 salp1 = cbet2 * somg12;
00526 calp1 = comg12 >= 0 ?
00527 sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) :
00528 sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
00529
00530 real
00531 ssig12 = Math::hypot(salp1, calp1),
00532 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
00533
00534 if (shortline && ssig12 < _etol2) {
00535
00536 salp2 = cbet1 * somg12;
00537 calp2 = sbet12 - cbet1 * sbet2 * sq(somg12) / (1 + comg12);
00538 SinCosNorm(salp2, calp2);
00539
00540 sig12 = atan2(ssig12, csig12);
00541 } else if (csig12 >= 0 ||
00542 ssig12 >= 3 * abs(_f) * Math::pi<real>() * sq(cbet1)) {
00543
00544
00545 } else {
00546
00547
00548 real x, y, lamscale, betscale;
00549 if (_f >= 0) {
00550
00551 {
00552 real
00553 k2 = sq(sbet1) * _ep2,
00554 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
00555 lamscale = _f * cbet1 * A3f(eps) * Math::pi<real>();
00556 }
00557 betscale = lamscale * cbet1;
00558
00559 x = (lam12 - Math::pi<real>()) / lamscale;
00560 y = sbet12a / betscale;
00561 } else {
00562
00563 real
00564 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
00565 bet12a = atan2(sbet12a, cbet12a);
00566 real m0, dummy;
00567
00568
00569 Lengths(_n, Math::pi<real>() + bet12a, sbet1, -cbet1, sbet2, cbet2,
00570 cbet1, cbet2, dummy, x, m0, false, dummy, dummy, C1a, C2a);
00571 x = -1 + x/(_f1 * cbet1 * cbet2 * m0 * Math::pi<real>());
00572 betscale = x < -real(0.01) ? sbet12a / x :
00573 -_f * sq(cbet1) * Math::pi<real>();
00574 lamscale = betscale / cbet1;
00575 y = (lam12 - Math::pi<real>()) / lamscale;
00576 }
00577
00578 if (y > -tol1 && x > -1 - xthresh) {
00579
00580 if (_f >= 0) {
00581 salp1 = min(real( 1), -x); calp1 = - sqrt(1 - sq(salp1));
00582 } else {
00583 calp1 = max(real(x > -tol1 ? 0 : -1), x);
00584 salp1 = sqrt(1 - sq(calp1));
00585 }
00586 } else {
00587
00588 real k = Astroid(x, y);
00589
00590 real
00591 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k ),
00592 somg12 = sin(omg12a), comg12 = -cos(omg12a);
00593
00594 salp1 = cbet2 * somg12;
00595 calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
00596 }
00597 }
00598 SinCosNorm(salp1, calp1);
00599 return sig12;
00600 }
00601
00602 Math::real Geodesic::Lambda12(real sbet1, real cbet1, real sbet2, real cbet2,
00603 real salp1, real calp1,
00604 real& salp2, real& calp2,
00605 real& sig12,
00606 real& ssig1, real& csig1,
00607 real& ssig2, real& csig2,
00608 real& eps, real& domg12,
00609 bool diffp, real& dlam12,
00610
00611 real C1a[], real C2a[], real C3a[]) const
00612 throw() {
00613
00614 if (sbet1 == 0 && calp1 == 0)
00615
00616
00617 calp1 = -eps2;
00618
00619 real
00620
00621 salp0 = salp1 * cbet1,
00622 calp0 = Math::hypot(calp1, salp1 * sbet1);
00623
00624 real somg1, comg1, somg2, comg2, omg12, lam12;
00625
00626
00627 ssig1 = sbet1; somg1 = salp0 * sbet1;
00628 csig1 = comg1 = calp1 * cbet1;
00629 SinCosNorm(ssig1, csig1);
00630 SinCosNorm(somg1, comg1);
00631
00632
00633
00634
00635
00636 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
00637
00638
00639
00640
00641 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ?
00642 sqrt(sq(calp1 * cbet1) + (cbet1 < -sbet1 ?
00643 (cbet2 - cbet1) * (cbet1 + cbet2) :
00644 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
00645 abs(calp1);
00646
00647
00648 ssig2 = sbet2; somg2 = salp0 * sbet2;
00649 csig2 = comg2 = calp2 * cbet2;
00650 SinCosNorm(ssig2, csig2);
00651 SinCosNorm(somg2, comg2);
00652
00653
00654 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
00655 csig1 * csig2 + ssig1 * ssig2);
00656
00657
00658 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, real(0)),
00659 comg1 * comg2 + somg1 * somg2);
00660 real B312, h0;
00661 real k2 = sq(calp0) * _ep2;
00662 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
00663 C3f(eps, C3a);
00664 B312 = (SinCosSeries(true, ssig2, csig2, C3a, nC3-1) -
00665 SinCosSeries(true, ssig1, csig1, C3a, nC3-1));
00666 h0 = -_f * A3f(eps);
00667 domg12 = salp0 * h0 * (sig12 + B312);
00668 lam12 = omg12 + domg12;
00669
00670 if (diffp) {
00671 if (calp2 == 0)
00672 dlam12 = - 2 * sqrt(1 - _e2 * sq(cbet1)) / sbet1;
00673 else {
00674 real dummy;
00675 Lengths(eps, sig12, ssig1, csig1, ssig2, csig2,
00676 cbet1, cbet2, dummy, dlam12, dummy,
00677 false, dummy, dummy, C1a, C2a);
00678 dlam12 /= calp2 * cbet2;
00679 }
00680 }
00681
00682 return lam12;
00683 }
00684
00685 Math::real Geodesic::A3f(real eps) const throw() {
00686
00687 real v = 0;
00688 for (int i = nA3x; i; )
00689 v = eps * v + _A3x[--i];
00690 return v;
00691 }
00692
00693 void Geodesic::C3f(real eps, real c[]) const throw() {
00694
00695
00696 for (int j = nC3x, k = nC3 - 1; k; ) {
00697 real t = 0;
00698 for (int i = nC3 - k; i; --i)
00699 t = eps * t + _C3x[--j];
00700 c[k--] = t;
00701 }
00702
00703 real mult = 1;
00704 for (int k = 1; k < nC3; ) {
00705 mult *= eps;
00706 c[k++] *= mult;
00707 }
00708 }
00709
00710 void Geodesic::C4f(real k2, real c[]) const throw() {
00711
00712
00713 for (int j = nC4x, k = nC4; k; ) {
00714 real t = 0;
00715 for (int i = nC4 - k + 1; i; --i)
00716 t = k2 * t + _C4x[--j];
00717 c[--k] = t;
00718 }
00719
00720 real mult = 1;
00721 for (int k = 1; k < nC4; ) {
00722 mult *= k2;
00723 c[k++] *= mult;
00724 }
00725 }
00726
00727
00728
00729
00730 Math::real Geodesic::A1m1f(real eps) throw() {
00731 real
00732 eps2 = sq(eps),
00733 t;
00734 switch (nA1/2) {
00735 case 0:
00736 t = 0;
00737 break;
00738 case 1:
00739 t = eps2/4;
00740 break;
00741 case 2:
00742 t = eps2*(eps2+16)/64;
00743 break;
00744 case 3:
00745 t = eps2*(eps2*(eps2+4)+64)/256;
00746 break;
00747 case 4:
00748 t = eps2*(eps2*(eps2*(25*eps2+64)+256)+4096)/16384;
00749 break;
00750 default:
00751 STATIC_ASSERT(nA1 >= 0 && nA1 <= 8, "Bad value of nA1");
00752 t = 0;
00753 }
00754 return (t + eps) / (1 - eps);
00755 }
00756
00757
00758 void Geodesic::C1f(real eps, real c[]) throw() {
00759 real
00760 eps2 = sq(eps),
00761 d = eps;
00762 switch (nC1) {
00763 case 0:
00764 break;
00765 case 1:
00766 c[1] = -d/2;
00767 break;
00768 case 2:
00769 c[1] = -d/2;
00770 d *= eps;
00771 c[2] = -d/16;
00772 break;
00773 case 3:
00774 c[1] = d*(3*eps2-8)/16;
00775 d *= eps;
00776 c[2] = -d/16;
00777 d *= eps;
00778 c[3] = -d/48;
00779 break;
00780 case 4:
00781 c[1] = d*(3*eps2-8)/16;
00782 d *= eps;
00783 c[2] = d*(eps2-2)/32;
00784 d *= eps;
00785 c[3] = -d/48;
00786 d *= eps;
00787 c[4] = -5*d/512;
00788 break;
00789 case 5:
00790 c[1] = d*((6-eps2)*eps2-16)/32;
00791 d *= eps;
00792 c[2] = d*(eps2-2)/32;
00793 d *= eps;
00794 c[3] = d*(9*eps2-16)/768;
00795 d *= eps;
00796 c[4] = -5*d/512;
00797 d *= eps;
00798 c[5] = -7*d/1280;
00799 break;
00800 case 6:
00801 c[1] = d*((6-eps2)*eps2-16)/32;
00802 d *= eps;
00803 c[2] = d*((64-9*eps2)*eps2-128)/2048;
00804 d *= eps;
00805 c[3] = d*(9*eps2-16)/768;
00806 d *= eps;
00807 c[4] = d*(3*eps2-5)/512;
00808 d *= eps;
00809 c[5] = -7*d/1280;
00810 d *= eps;
00811 c[6] = -7*d/2048;
00812 break;
00813 case 7:
00814 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048;
00815 d *= eps;
00816 c[2] = d*((64-9*eps2)*eps2-128)/2048;
00817 d *= eps;
00818 c[3] = d*((72-9*eps2)*eps2-128)/6144;
00819 d *= eps;
00820 c[4] = d*(3*eps2-5)/512;
00821 d *= eps;
00822 c[5] = d*(35*eps2-56)/10240;
00823 d *= eps;
00824 c[6] = -7*d/2048;
00825 d *= eps;
00826 c[7] = -33*d/14336;
00827 break;
00828 case 8:
00829 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048;
00830 d *= eps;
00831 c[2] = d*(eps2*(eps2*(7*eps2-18)+128)-256)/4096;
00832 d *= eps;
00833 c[3] = d*((72-9*eps2)*eps2-128)/6144;
00834 d *= eps;
00835 c[4] = d*((96-11*eps2)*eps2-160)/16384;
00836 d *= eps;
00837 c[5] = d*(35*eps2-56)/10240;
00838 d *= eps;
00839 c[6] = d*(9*eps2-14)/4096;
00840 d *= eps;
00841 c[7] = -33*d/14336;
00842 d *= eps;
00843 c[8] = -429*d/262144;
00844 break;
00845 default:
00846 STATIC_ASSERT(nC1 >= 0 && nC1 <= 8, "Bad value of nC1");
00847 }
00848 }
00849
00850
00851 void Geodesic::C1pf(real eps, real c[]) throw() {
00852 real
00853 eps2 = sq(eps),
00854 d = eps;
00855 switch (nC1p) {
00856 case 0:
00857 break;
00858 case 1:
00859 c[1] = d/2;
00860 break;
00861 case 2:
00862 c[1] = d/2;
00863 d *= eps;
00864 c[2] = 5*d/16;
00865 break;
00866 case 3:
00867 c[1] = d*(16-9*eps2)/32;
00868 d *= eps;
00869 c[2] = 5*d/16;
00870 d *= eps;
00871 c[3] = 29*d/96;
00872 break;
00873 case 4:
00874 c[1] = d*(16-9*eps2)/32;
00875 d *= eps;
00876 c[2] = d*(30-37*eps2)/96;
00877 d *= eps;
00878 c[3] = 29*d/96;
00879 d *= eps;
00880 c[4] = 539*d/1536;
00881 break;
00882 case 5:
00883 c[1] = d*(eps2*(205*eps2-432)+768)/1536;
00884 d *= eps;
00885 c[2] = d*(30-37*eps2)/96;
00886 d *= eps;
00887 c[3] = d*(116-225*eps2)/384;
00888 d *= eps;
00889 c[4] = 539*d/1536;
00890 d *= eps;
00891 c[5] = 3467*d/7680;
00892 break;
00893 case 6:
00894 c[1] = d*(eps2*(205*eps2-432)+768)/1536;
00895 d *= eps;
00896 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
00897 d *= eps;
00898 c[3] = d*(116-225*eps2)/384;
00899 d *= eps;
00900 c[4] = d*(2695-7173*eps2)/7680;
00901 d *= eps;
00902 c[5] = 3467*d/7680;
00903 d *= eps;
00904 c[6] = 38081*d/61440;
00905 break;
00906 case 7:
00907 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728;
00908 d *= eps;
00909 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
00910 d *= eps;
00911 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288;
00912 d *= eps;
00913 c[4] = d*(2695-7173*eps2)/7680;
00914 d *= eps;
00915 c[5] = d*(41604-141115*eps2)/92160;
00916 d *= eps;
00917 c[6] = 38081*d/61440;
00918 d *= eps;
00919 c[7] = 459485*d/516096;
00920 break;
00921 case 8:
00922 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728;
00923 d *= eps;
00924 c[2] = d*(eps2*((120150-86171*eps2)*eps2-142080)+115200)/368640;
00925 d *= eps;
00926 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288;
00927 d *= eps;
00928 c[4] = d*(eps2*(1082857*eps2-688608)+258720)/737280;
00929 d *= eps;
00930 c[5] = d*(41604-141115*eps2)/92160;
00931 d *= eps;
00932 c[6] = d*(533134-2200311*eps2)/860160;
00933 d *= eps;
00934 c[7] = 459485*d/516096;
00935 d *= eps;
00936 c[8] = 109167851*d/82575360;
00937 break;
00938 default:
00939 STATIC_ASSERT(nC1p >= 0 && nC1p <= 8, "Bad value of nC1p");
00940 }
00941 }
00942
00943
00944 Math::real Geodesic::A2m1f(real eps) throw() {
00945 real
00946 eps2 = sq(eps),
00947 t;
00948 switch (nA2/2) {
00949 case 0:
00950 t = 0;
00951 break;
00952 case 1:
00953 t = eps2/4;
00954 break;
00955 case 2:
00956 t = eps2*(9*eps2+16)/64;
00957 break;
00958 case 3:
00959 t = eps2*(eps2*(25*eps2+36)+64)/256;
00960 break;
00961 case 4:
00962 t = eps2*(eps2*(eps2*(1225*eps2+1600)+2304)+4096)/16384;
00963 break;
00964 default:
00965 STATIC_ASSERT(nA2 >= 0 && nA2 <= 8, "Bad value of nA2");
00966 t = 0;
00967 }
00968 return t * (1 - eps) - eps;
00969 }
00970
00971
00972 void Geodesic::C2f(real eps, real c[]) throw() {
00973 real
00974 eps2 = sq(eps),
00975 d = eps;
00976 switch (nC2) {
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] = 3*d/16;
00986 break;
00987 case 3:
00988 c[1] = d*(eps2+8)/16;
00989 d *= eps;
00990 c[2] = 3*d/16;
00991 d *= eps;
00992 c[3] = 5*d/48;
00993 break;
00994 case 4:
00995 c[1] = d*(eps2+8)/16;
00996 d *= eps;
00997 c[2] = d*(eps2+6)/32;
00998 d *= eps;
00999 c[3] = 5*d/48;
01000 d *= eps;
01001 c[4] = 35*d/512;
01002 break;
01003 case 5:
01004 c[1] = d*(eps2*(eps2+2)+16)/32;
01005 d *= eps;
01006 c[2] = d*(eps2+6)/32;
01007 d *= eps;
01008 c[3] = d*(15*eps2+80)/768;
01009 d *= eps;
01010 c[4] = 35*d/512;
01011 d *= eps;
01012 c[5] = 63*d/1280;
01013 break;
01014 case 6:
01015 c[1] = d*(eps2*(eps2+2)+16)/32;
01016 d *= eps;
01017 c[2] = d*(eps2*(35*eps2+64)+384)/2048;
01018 d *= eps;
01019 c[3] = d*(15*eps2+80)/768;
01020 d *= eps;
01021 c[4] = d*(7*eps2+35)/512;
01022 d *= eps;
01023 c[5] = 63*d/1280;
01024 d *= eps;
01025 c[6] = 77*d/2048;
01026 break;
01027 case 7:
01028 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048;
01029 d *= eps;
01030 c[2] = d*(eps2*(35*eps2+64)+384)/2048;
01031 d *= eps;
01032 c[3] = d*(eps2*(69*eps2+120)+640)/6144;
01033 d *= eps;
01034 c[4] = d*(7*eps2+35)/512;
01035 d *= eps;
01036 c[5] = d*(105*eps2+504)/10240;
01037 d *= eps;
01038 c[6] = 77*d/2048;
01039 d *= eps;
01040 c[7] = 429*d/14336;
01041 break;
01042 case 8:
01043 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048;
01044 d *= eps;
01045 c[2] = d*(eps2*(eps2*(47*eps2+70)+128)+768)/4096;
01046 d *= eps;
01047 c[3] = d*(eps2*(69*eps2+120)+640)/6144;
01048 d *= eps;
01049 c[4] = d*(eps2*(133*eps2+224)+1120)/16384;
01050 d *= eps;
01051 c[5] = d*(105*eps2+504)/10240;
01052 d *= eps;
01053 c[6] = d*(33*eps2+154)/4096;
01054 d *= eps;
01055 c[7] = 429*d/14336;
01056 d *= eps;
01057 c[8] = 6435*d/262144;
01058 break;
01059 default:
01060 STATIC_ASSERT(nC2 >= 0 && nC2 <= 8, "Bad value of nC2");
01061 }
01062 }
01063
01064
01065 void Geodesic::A3coeff() throw() {
01066 switch (nA3) {
01067 case 0:
01068 break;
01069 case 1:
01070 _A3x[0] = 1;
01071 break;
01072 case 2:
01073 _A3x[0] = 1;
01074 _A3x[1] = -1/real(2);
01075 break;
01076 case 3:
01077 _A3x[0] = 1;
01078 _A3x[1] = (_n-1)/2;
01079 _A3x[2] = -1/real(4);
01080 break;
01081 case 4:
01082 _A3x[0] = 1;
01083 _A3x[1] = (_n-1)/2;
01084 _A3x[2] = (-_n-2)/8;
01085 _A3x[3] = -1/real(16);
01086 break;
01087 case 5:
01088 _A3x[0] = 1;
01089 _A3x[1] = (_n-1)/2;
01090 _A3x[2] = (_n*(3*_n-1)-2)/8;
01091 _A3x[3] = (-3*_n-1)/16;
01092 _A3x[4] = -3/real(64);
01093 break;
01094 case 6:
01095 _A3x[0] = 1;
01096 _A3x[1] = (_n-1)/2;
01097 _A3x[2] = (_n*(3*_n-1)-2)/8;
01098 _A3x[3] = ((-_n-3)*_n-1)/16;
01099 _A3x[4] = (-2*_n-3)/64;
01100 _A3x[5] = -3/real(128);
01101 break;
01102 case 7:
01103 _A3x[0] = 1;
01104 _A3x[1] = (_n-1)/2;
01105 _A3x[2] = (_n*(3*_n-1)-2)/8;
01106 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16;
01107 _A3x[4] = ((-10*_n-2)*_n-3)/64;
01108 _A3x[5] = (-5*_n-3)/128;
01109 _A3x[6] = -5/real(256);
01110 break;
01111 case 8:
01112 _A3x[0] = 1;
01113 _A3x[1] = (_n-1)/2;
01114 _A3x[2] = (_n*(3*_n-1)-2)/8;
01115 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16;
01116 _A3x[4] = (_n*((-5*_n-20)*_n-4)-6)/128;
01117 _A3x[5] = ((-5*_n-10)*_n-6)/256;
01118 _A3x[6] = (-15*_n-20)/1024;
01119 _A3x[7] = -25/real(2048);
01120 break;
01121 default:
01122 STATIC_ASSERT(nA3 >= 0 && nA3 <= 8, "Bad value of nA3");
01123 }
01124 }
01125
01126
01127 void Geodesic::C3coeff() throw() {
01128 switch (nC3) {
01129 case 0:
01130 break;
01131 case 1:
01132 break;
01133 case 2:
01134 _C3x[0] = 1/real(4);
01135 break;
01136 case 3:
01137 _C3x[0] = (1-_n)/4;
01138 _C3x[1] = 1/real(8);
01139 _C3x[2] = 1/real(16);
01140 break;
01141 case 4:
01142 _C3x[0] = (1-_n)/4;
01143 _C3x[1] = 1/real(8);
01144 _C3x[2] = 3/real(64);
01145 _C3x[3] = (2-3*_n)/32;
01146 _C3x[4] = 3/real(64);
01147 _C3x[5] = 5/real(192);
01148 break;
01149 case 5:
01150 _C3x[0] = (1-_n)/4;
01151 _C3x[1] = (1-_n*_n)/8;
01152 _C3x[2] = (3*_n+3)/64;
01153 _C3x[3] = 5/real(128);
01154 _C3x[4] = ((_n-3)*_n+2)/32;
01155 _C3x[5] = (3-2*_n)/64;
01156 _C3x[6] = 3/real(128);
01157 _C3x[7] = (5-9*_n)/192;
01158 _C3x[8] = 3/real(128);
01159 _C3x[9] = 7/real(512);
01160 break;
01161 case 6:
01162 _C3x[0] = (1-_n)/4;
01163 _C3x[1] = (1-_n*_n)/8;
01164 _C3x[2] = ((3-_n)*_n+3)/64;
01165 _C3x[3] = (2*_n+5)/128;
01166 _C3x[4] = 3/real(128);
01167 _C3x[5] = ((_n-3)*_n+2)/32;
01168 _C3x[6] = ((-3*_n-2)*_n+3)/64;
01169 _C3x[7] = (_n+3)/128;
01170 _C3x[8] = 5/real(256);
01171 _C3x[9] = (_n*(5*_n-9)+5)/192;
01172 _C3x[10] = (9-10*_n)/384;
01173 _C3x[11] = 7/real(512);
01174 _C3x[12] = (7-14*_n)/512;
01175 _C3x[13] = 7/real(512);
01176 _C3x[14] = 21/real(2560);
01177 break;
01178 case 7:
01179 _C3x[0] = (1-_n)/4;
01180 _C3x[1] = (1-_n*_n)/8;
01181 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64;
01182 _C3x[3] = (_n*(2*_n+2)+5)/128;
01183 _C3x[4] = (11*_n+12)/512;
01184 _C3x[5] = 21/real(1024);
01185 _C3x[6] = ((_n-3)*_n+2)/32;
01186 _C3x[7] = (_n*(_n*(2*_n-3)-2)+3)/64;
01187 _C3x[8] = ((2-9*_n)*_n+6)/256;
01188 _C3x[9] = (_n+5)/256;
01189 _C3x[10] = 27/real(2048);
01190 _C3x[11] = (_n*((5-_n)*_n-9)+5)/192;
01191 _C3x[12] = ((-6*_n-10)*_n+9)/384;
01192 _C3x[13] = (21-4*_n)/1536;
01193 _C3x[14] = 3/real(256);
01194 _C3x[15] = (_n*(10*_n-14)+7)/512;
01195 _C3x[16] = (7-10*_n)/512;
01196 _C3x[17] = 9/real(1024);
01197 _C3x[18] = (21-45*_n)/2560;
01198 _C3x[19] = 9/real(1024);
01199 _C3x[20] = 11/real(2048);
01200 break;
01201 case 8:
01202 _C3x[0] = (1-_n)/4;
01203 _C3x[1] = (1-_n*_n)/8;
01204 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64;
01205 _C3x[3] = (_n*((2-2*_n)*_n+2)+5)/128;
01206 _C3x[4] = (_n*(3*_n+11)+12)/512;
01207 _C3x[5] = (10*_n+21)/1024;
01208 _C3x[6] = 243/real(16384);
01209 _C3x[7] = ((_n-3)*_n+2)/32;
01210 _C3x[8] = (_n*(_n*(2*_n-3)-2)+3)/64;
01211 _C3x[9] = (_n*((-6*_n-9)*_n+2)+6)/256;
01212 _C3x[10] = ((1-2*_n)*_n+5)/256;
01213 _C3x[11] = (69*_n+108)/8192;
01214 _C3x[12] = 187/real(16384);
01215 _C3x[13] = (_n*((5-_n)*_n-9)+5)/192;
01216 _C3x[14] = (_n*(_n*(10*_n-6)-10)+9)/384;
01217 _C3x[15] = ((-77*_n-8)*_n+42)/3072;
01218 _C3x[16] = (12-_n)/1024;
01219 _C3x[17] = 139/real(16384);
01220 _C3x[18] = (_n*((20-7*_n)*_n-28)+14)/1024;
01221 _C3x[19] = ((-7*_n-40)*_n+28)/2048;
01222 _C3x[20] = (72-43*_n)/8192;
01223 _C3x[21] = 127/real(16384);
01224 _C3x[22] = (_n*(75*_n-90)+42)/5120;
01225 _C3x[23] = (9-15*_n)/1024;
01226 _C3x[24] = 99/real(16384);
01227 _C3x[25] = (44-99*_n)/8192;
01228 _C3x[26] = 99/real(16384);
01229 _C3x[27] = 429/real(114688);
01230 break;
01231 default:
01232 STATIC_ASSERT(nC3 >= 0 && nC3 <= 8, "Bad value of nC3");
01233 }
01234 }
01235
01236
01237 void Geodesic::C4coeff() throw() {
01238 switch (nC4) {
01239 case 0:
01240 break;
01241 case 1:
01242 _C4x[0] = 2/real(3);
01243 break;
01244 case 2:
01245 _C4x[0] = (10-_ep2)/15;
01246 _C4x[1] = -1/real(20);
01247 _C4x[2] = 1/real(180);
01248 break;
01249 case 3:
01250 _C4x[0] = (_ep2*(4*_ep2-7)+70)/105;
01251 _C4x[1] = (4*_ep2-7)/140;
01252 _C4x[2] = 1/real(42);
01253 _C4x[3] = (7-4*_ep2)/1260;
01254 _C4x[4] = -1/real(252);
01255 _C4x[5] = 1/real(2100);
01256 break;
01257 case 4:
01258 _C4x[0] = (_ep2*((12-8*_ep2)*_ep2-21)+210)/315;
01259 _C4x[1] = ((12-8*_ep2)*_ep2-21)/420;
01260 _C4x[2] = (3-2*_ep2)/126;
01261 _C4x[3] = -1/real(72);
01262 _C4x[4] = (_ep2*(8*_ep2-12)+21)/3780;
01263 _C4x[5] = (2*_ep2-3)/756;
01264 _C4x[6] = 1/real(360);
01265 _C4x[7] = (3-2*_ep2)/6300;
01266 _C4x[8] = -1/real(1800);
01267 _C4x[9] = 1/real(17640);
01268 break;
01269 case 5:
01270 _C4x[0] = (_ep2*(_ep2*(_ep2*(64*_ep2-88)+132)-231)+2310)/3465;
01271 _C4x[1] = (_ep2*(_ep2*(64*_ep2-88)+132)-231)/4620;
01272 _C4x[2] = (_ep2*(16*_ep2-22)+33)/1386;
01273 _C4x[3] = (8*_ep2-11)/792;
01274 _C4x[4] = 1/real(110);
01275 _C4x[5] = (_ep2*((88-64*_ep2)*_ep2-132)+231)/41580;
01276 _C4x[6] = ((22-16*_ep2)*_ep2-33)/8316;
01277 _C4x[7] = (11-8*_ep2)/3960;
01278 _C4x[8] = -1/real(495);
01279 _C4x[9] = (_ep2*(16*_ep2-22)+33)/69300;
01280 _C4x[10] = (8*_ep2-11)/19800;
01281 _C4x[11] = 1/real(1925);
01282 _C4x[12] = (11-8*_ep2)/194040;
01283 _C4x[13] = -1/real(10780);
01284 _C4x[14] = 1/real(124740);
01285 break;
01286 case 6:
01287 _C4x[0] = (_ep2*(_ep2*(_ep2*((832-640*_ep2)*_ep2-1144)+1716)-3003)+
01288 30030)/45045;
01289 _C4x[1] = (_ep2*(_ep2*((832-640*_ep2)*_ep2-1144)+1716)-3003)/60060;
01290 _C4x[2] = (_ep2*((208-160*_ep2)*_ep2-286)+429)/18018;
01291 _C4x[3] = ((104-80*_ep2)*_ep2-143)/10296;
01292 _C4x[4] = (13-10*_ep2)/1430;
01293 _C4x[5] = -1/real(156);
01294 _C4x[6] = (_ep2*(_ep2*(_ep2*(640*_ep2-832)+1144)-1716)+3003)/540540;
01295 _C4x[7] = (_ep2*(_ep2*(160*_ep2-208)+286)-429)/108108;
01296 _C4x[8] = (_ep2*(80*_ep2-104)+143)/51480;
01297 _C4x[9] = (10*_ep2-13)/6435;
01298 _C4x[10] = 5/real(3276);
01299 _C4x[11] = (_ep2*((208-160*_ep2)*_ep2-286)+429)/900900;
01300 _C4x[12] = ((104-80*_ep2)*_ep2-143)/257400;
01301 _C4x[13] = (13-10*_ep2)/25025;
01302 _C4x[14] = -1/real(2184);
01303 _C4x[15] = (_ep2*(80*_ep2-104)+143)/2522520;
01304 _C4x[16] = (10*_ep2-13)/140140;
01305 _C4x[17] = 5/real(45864);
01306 _C4x[18] = (13-10*_ep2)/1621620;
01307 _C4x[19] = -1/real(58968);
01308 _C4x[20] = 1/real(792792);
01309 break;
01310 case 7:
01311 _C4x[0] = (_ep2*(_ep2*(_ep2*(_ep2*(_ep2*(512*_ep2-640)+832)-1144)+1716)-
01312 3003)+30030)/45045;
01313 _C4x[1] = (_ep2*(_ep2*(_ep2*(_ep2*(512*_ep2-640)+832)-1144)+1716)-
01314 3003)/60060;
01315 _C4x[2] = (_ep2*(_ep2*(_ep2*(128*_ep2-160)+208)-286)+429)/18018;
01316 _C4x[3] = (_ep2*(_ep2*(64*_ep2-80)+104)-143)/10296;
01317 _C4x[4] = (_ep2*(8*_ep2-10)+13)/1430;
01318 _C4x[5] = (4*_ep2-5)/780;
01319 _C4x[6] = 1/real(210);
01320 _C4x[7] = (_ep2*(_ep2*(_ep2*((640-512*_ep2)*_ep2-832)+1144)-1716)+
01321 3003)/540540;
01322 _C4x[8] = (_ep2*(_ep2*((160-128*_ep2)*_ep2-208)+286)-429)/108108;
01323 _C4x[9] = (_ep2*((80-64*_ep2)*_ep2-104)+143)/51480;
01324 _C4x[10] = ((10-8*_ep2)*_ep2-13)/6435;
01325 _C4x[11] = (5-4*_ep2)/3276;
01326 _C4x[12] = -1/real(840);
01327 _C4x[13] = (_ep2*(_ep2*(_ep2*(128*_ep2-160)+208)-286)+429)/900900;
01328 _C4x[14] = (_ep2*(_ep2*(64*_ep2-80)+104)-143)/257400;
01329 _C4x[15] = (_ep2*(8*_ep2-10)+13)/25025;
01330 _C4x[16] = (4*_ep2-5)/10920;
01331 _C4x[17] = 1/real(2520);
01332 _C4x[18] = (_ep2*((80-64*_ep2)*_ep2-104)+143)/2522520;
01333 _C4x[19] = ((10-8*_ep2)*_ep2-13)/140140;
01334 _C4x[20] = (5-4*_ep2)/45864;
01335 _C4x[21] = -1/real(8820);
01336 _C4x[22] = (_ep2*(8*_ep2-10)+13)/1621620;
01337 _C4x[23] = (4*_ep2-5)/294840;
01338 _C4x[24] = 1/real(41580);
01339 _C4x[25] = (5-4*_ep2)/3963960;
01340 _C4x[26] = -1/real(304920);
01341 _C4x[27] = 1/real(4684680);
01342 break;
01343 case 8:
01344 _C4x[0] = (_ep2*(_ep2*(_ep2*(_ep2*(_ep2*((8704-7168*_ep2)*_ep2-10880)+
01345 14144)-19448)+29172)-51051)+510510)/765765;
01346 _C4x[1] = (_ep2*(_ep2*(_ep2*(_ep2*((8704-7168*_ep2)*_ep2-10880)+14144)-
01347 19448)+29172)-51051)/1021020;
01348 _C4x[2] = (_ep2*(_ep2*(_ep2*((2176-1792*_ep2)*_ep2-2720)+3536)-4862)+
01349 7293)/306306;
01350 _C4x[3] = (_ep2*(_ep2*((1088-896*_ep2)*_ep2-1360)+1768)-2431)/175032;
01351 _C4x[4] = (_ep2*((136-112*_ep2)*_ep2-170)+221)/24310;
01352 _C4x[5] = ((68-56*_ep2)*_ep2-85)/13260;
01353 _C4x[6] = (17-14*_ep2)/3570;
01354 _C4x[7] = -1/real(272);
01355 _C4x[8] = (_ep2*(_ep2*(_ep2*(_ep2*(_ep2*(7168*_ep2-8704)+10880)-14144)+
01356 19448)-29172)+51051)/9189180;
01357 _C4x[9] = (_ep2*(_ep2*(_ep2*(_ep2*(1792*_ep2-2176)+2720)-3536)+4862)-
01358 7293)/1837836;
01359 _C4x[10] = (_ep2*(_ep2*(_ep2*(896*_ep2-1088)+1360)-1768)+2431)/875160;
01360 _C4x[11] = (_ep2*(_ep2*(112*_ep2-136)+170)-221)/109395;
01361 _C4x[12] = (_ep2*(56*_ep2-68)+85)/55692;
01362 _C4x[13] = (14*_ep2-17)/14280;
01363 _C4x[14] = 7/real(7344);
01364 _C4x[15] = (_ep2*(_ep2*(_ep2*((2176-1792*_ep2)*_ep2-2720)+3536)-4862)+
01365 7293)/15315300;
01366 _C4x[16] = (_ep2*(_ep2*((1088-896*_ep2)*_ep2-1360)+1768)-2431)/4375800;
01367 _C4x[17] = (_ep2*((136-112*_ep2)*_ep2-170)+221)/425425;
01368 _C4x[18] = ((68-56*_ep2)*_ep2-85)/185640;
01369 _C4x[19] = (17-14*_ep2)/42840;
01370 _C4x[20] = -7/real(20400);
01371 _C4x[21] = (_ep2*(_ep2*(_ep2*(896*_ep2-1088)+1360)-1768)+2431)/42882840;
01372 _C4x[22] = (_ep2*(_ep2*(112*_ep2-136)+170)-221)/2382380;
01373 _C4x[23] = (_ep2*(56*_ep2-68)+85)/779688;
01374 _C4x[24] = (14*_ep2-17)/149940;
01375 _C4x[25] = 1/real(8976);
01376 _C4x[26] = (_ep2*((136-112*_ep2)*_ep2-170)+221)/27567540;
01377 _C4x[27] = ((68-56*_ep2)*_ep2-85)/5012280;
01378 _C4x[28] = (17-14*_ep2)/706860;
01379 _C4x[29] = -7/real(242352);
01380 _C4x[30] = (_ep2*(56*_ep2-68)+85)/67387320;
01381 _C4x[31] = (14*_ep2-17)/5183640;
01382 _C4x[32] = 7/real(1283568);
01383 _C4x[33] = (17-14*_ep2)/79639560;
01384 _C4x[34] = -1/real(1516944);
01385 _C4x[35] = 1/real(26254800);
01386 break;
01387 default:
01388 STATIC_ASSERT(nC3 >= 0 && nC4 <= 8, "Bad value of nC4");
01389 }
01390 }
01391 }