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/GeodesicExact.hpp>
00030 #include <GeographicLib/GeodesicLineExact.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 GeodesicExact::GeodesicExact(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 C4coeff();
00086 }
00087
00088 const GeodesicExact& GeodesicExact::WGS84() {
00089 static const GeodesicExact wgs84(Constants::WGS84_a(),
00090 Constants::WGS84_f());
00091 return wgs84;
00092 }
00093
00094 Math::real GeodesicExact::CosSeries(real sinx, real cosx,
00095 const real c[], int n) {
00096
00097
00098
00099
00100 c += n ;
00101 real
00102 ar = 2 * (cosx - sinx) * (cosx + sinx),
00103 y0 = n & 1 ? *--c : 0, y1 = 0;
00104
00105 n /= 2;
00106 while (n--) {
00107
00108 y1 = ar * y0 - y1 + *--c;
00109 y0 = ar * y1 - y0 + *--c;
00110 }
00111 return cosx * (y0 - y1);
00112 }
00113
00114 GeodesicLineExact GeodesicExact::Line(real lat1, real lon1, real azi1,
00115 unsigned caps) const {
00116 return GeodesicLineExact(*this, lat1, lon1, azi1, caps);
00117 }
00118
00119 Math::real GeodesicExact::GenDirect(real lat1, real lon1, real azi1,
00120 bool arcmode, real s12_a12,
00121 unsigned outmask,
00122 real& lat2, real& lon2, real& azi2,
00123 real& s12, real& m12,
00124 real& M12, real& M21,
00125 real& S12) const {
00126 return GeodesicLineExact(*this, lat1, lon1, azi1,
00127
00128 outmask | (arcmode ? NONE : DISTANCE_IN))
00129 .
00130 GenPosition(arcmode, s12_a12, outmask,
00131 lat2, lon2, azi2, s12, m12, M12, M21, S12);
00132 }
00133
00134 Math::real GeodesicExact::GenInverse(real lat1, real lon1,
00135 real lat2, real lon2,
00136 unsigned outmask,
00137 real& s12, real& azi1, real& azi2,
00138 real& m12, real& M12, real& M21,
00139 real& S12) const {
00140 outmask &= OUT_ALL;
00141
00142
00143
00144 real lon12 = Math::AngDiff(Math::AngNormalize(lon1),
00145 Math::AngNormalize(lon2));
00146
00147 lon12 = AngRound(lon12);
00148
00149 int lonsign = lon12 >= 0 ? 1 : -1;
00150 lon12 *= lonsign;
00151
00152 lat1 = AngRound(lat1);
00153 lat2 = AngRound(lat2);
00154
00155 int swapp = abs(lat1) >= abs(lat2) ? 1 : -1;
00156 if (swapp < 0) {
00157 lonsign *= -1;
00158 swap(lat1, lat2);
00159 }
00160
00161 int latsign = lat1 < 0 ? 1 : -1;
00162 lat1 *= latsign;
00163 lat2 *= latsign;
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176 real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x;
00177
00178
00179 EllipticFunction E(-_ep2);
00180
00181 phi = lat1 * Math::degree();
00182
00183 sbet1 = _f1 * sin(phi);
00184 cbet1 = lat1 == -90 ? tiny_ : cos(phi);
00185 SinCosNorm(sbet1, cbet1);
00186
00187 phi = lat2 * Math::degree();
00188
00189 sbet2 = _f1 * sin(phi);
00190 cbet2 = abs(lat2) == 90 ? tiny_ : cos(phi);
00191 SinCosNorm(sbet2, cbet2);
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 if (cbet1 < -sbet1) {
00202 if (cbet2 == cbet1)
00203 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
00204 } else {
00205 if (abs(sbet2) == -sbet1)
00206 cbet2 = cbet1;
00207 }
00208
00209 real
00210 dn1 = (_f >= 0 ? sqrt(1 + _ep2 * Math::sq(sbet1)) :
00211 sqrt(1 - _e2 * Math::sq(cbet1)) / _f1),
00212 dn2 = (_f >= 0 ? sqrt(1 + _ep2 * Math::sq(sbet2)) :
00213 sqrt(1 - _e2 * Math::sq(cbet2)) / _f1);
00214
00215 real
00216 lam12 = lon12 * Math::degree(),
00217 slam12 = abs(lon12) == 180 ? 0 : sin(lam12),
00218 clam12 = cos(lam12);
00219
00220
00221 real a12, sig12, calp1, salp1, calp2 = 0, salp2 = 0;
00222
00223 bool meridian = lat1 == -90 || slam12 == 0;
00224
00225 if (meridian) {
00226
00227
00228
00229
00230 calp1 = clam12; salp1 = slam12;
00231 calp2 = 1; salp2 = 0;
00232
00233 real
00234
00235 ssig1 = sbet1, csig1 = calp1 * cbet1,
00236 ssig2 = sbet2, csig2 = calp2 * cbet2;
00237
00238
00239 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
00240 csig1 * csig2 + ssig1 * ssig2);
00241 {
00242 real dummy;
00243 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
00244 cbet1, cbet2, s12x, m12x, dummy,
00245 (outmask & GEODESICSCALE) != 0U, M12, M21);
00246 }
00247
00248
00249
00250
00251
00252
00253
00254 if (sig12 < 1 || m12x >= 0) {
00255 m12x *= _b;
00256 s12x *= _b;
00257 a12 = sig12 / Math::degree();
00258 } else
00259
00260 meridian = false;
00261 }
00262
00263 real omg12 = 0;
00264 if (!meridian &&
00265 sbet1 == 0 &&
00266
00267 (_f <= 0 || lam12 <= Math::pi() - _f * Math::pi())) {
00268
00269
00270 calp1 = calp2 = 0; salp1 = salp2 = 1;
00271 s12x = _a * lam12;
00272 sig12 = omg12 = lam12 / _f1;
00273 m12x = _b * sin(sig12);
00274 if (outmask & GEODESICSCALE)
00275 M12 = M21 = cos(sig12);
00276 a12 = lon12 / _f1;
00277
00278 } else if (!meridian) {
00279
00280
00281
00282
00283
00284 real dnm;
00285 sig12 = InverseStart(E, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
00286 lam12,
00287 salp1, calp1, salp2, calp2, dnm);
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;
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
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 real dv;
00341 real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
00342 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
00343 E, omg12, numit < maxit1_, dv) - lam12;
00344
00345
00346 if (tripb || !(abs(v) >= (tripn ? 8 : 2) * tol0_)) break;
00347
00348 if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
00349 { salp1b = salp1; calp1b = calp1; }
00350 else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
00351 { salp1a = salp1; calp1a = calp1; }
00352 if (numit < maxit1_ && dv > 0) {
00353 real
00354 dalp1 = -v/dv;
00355 real
00356 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
00357 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
00358 if (nsalp1 > 0 && abs(dalp1) < Math::pi()) {
00359 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
00360 salp1 = nsalp1;
00361 SinCosNorm(salp1, calp1);
00362
00363
00364
00365 tripn = abs(v) <= 16 * tol0_;
00366 continue;
00367 }
00368 }
00369
00370
00371
00372
00373
00374
00375
00376
00377 salp1 = (salp1a + salp1b)/2;
00378 calp1 = (calp1a + calp1b)/2;
00379 SinCosNorm(salp1, calp1);
00380 tripn = false;
00381 tripb = (abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
00382 abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
00383 }
00384 {
00385 real dummy;
00386 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
00387 cbet1, cbet2, s12x, m12x, dummy,
00388 (outmask & GEODESICSCALE) != 0U, M12, M21);
00389 }
00390 m12x *= _b;
00391 s12x *= _b;
00392 a12 = sig12 / Math::degree();
00393 }
00394 }
00395
00396 if (outmask & DISTANCE)
00397 s12 = 0 + s12x;
00398
00399 if (outmask & REDUCEDLENGTH)
00400 m12 = 0 + m12x;
00401
00402 if (outmask & AREA) {
00403 real
00404
00405 salp0 = salp1 * cbet1,
00406 calp0 = Math::hypot(calp1, salp1 * sbet1);
00407 real alp12;
00408 if (calp0 != 0 && salp0 != 0) {
00409 real
00410
00411 ssig1 = sbet1, csig1 = calp1 * cbet1,
00412 ssig2 = sbet2, csig2 = calp2 * cbet2,
00413 k2 = Math::sq(calp0) * _ep2,
00414 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
00415
00416 A4 = Math::sq(_a) * calp0 * salp0 * _e2;
00417 SinCosNorm(ssig1, csig1);
00418 SinCosNorm(ssig2, csig2);
00419 real C4a[nC4_];
00420 C4f(eps, C4a);
00421 real
00422 B41 = CosSeries(ssig1, csig1, C4a, nC4_),
00423 B42 = CosSeries(ssig2, csig2, C4a, nC4_);
00424 S12 = A4 * (B42 - B41);
00425 } else
00426
00427 S12 = 0;
00428
00429 if (!meridian &&
00430 omg12 < real(0.75) * Math::pi() &&
00431 sbet2 - sbet1 < real(1.75)) {
00432
00433
00434
00435 real
00436 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
00437 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
00438 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
00439 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
00440 } else {
00441
00442 real
00443 salp12 = salp2 * calp1 - calp2 * salp1,
00444 calp12 = calp2 * calp1 + salp2 * salp1;
00445
00446
00447
00448
00449 if (salp12 == 0 && calp12 < 0) {
00450 salp12 = tiny_ * calp1;
00451 calp12 = -1;
00452 }
00453 alp12 = atan2(salp12, calp12);
00454 }
00455 S12 += _c2 * alp12;
00456 S12 *= swapp * lonsign * latsign;
00457
00458 S12 += 0;
00459 }
00460
00461
00462 if (swapp < 0) {
00463 swap(salp1, salp2);
00464 swap(calp1, calp2);
00465 if (outmask & GEODESICSCALE)
00466 swap(M12, M21);
00467 }
00468
00469 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
00470 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
00471
00472 if (outmask & AZIMUTH) {
00473
00474 azi1 = 0 - atan2(-salp1, calp1) / Math::degree();
00475 azi2 = 0 - atan2(-salp2, calp2) / Math::degree();
00476 }
00477
00478
00479 return a12;
00480 }
00481
00482 void GeodesicExact::Lengths(const EllipticFunction& E,
00483 real sig12,
00484 real ssig1, real csig1, real dn1,
00485 real ssig2, real csig2, real dn2,
00486 real cbet1, real cbet2,
00487 real& s12b, real& m12b, real& m0,
00488 bool scalep, real& M12, real& M21) const {
00489
00490
00491
00492
00493
00494 m0 = - E.k2() * E.D() / (Math::pi() / 2);
00495 real J12 = m0 *
00496 (sig12 + E.deltaD(ssig2, csig2, dn2) - E.deltaD(ssig1, csig1, dn1));
00497
00498
00499
00500 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12;
00501
00502 s12b = E.E() / (Math::pi() / 2) *
00503 (sig12 + E.deltaE(ssig2, csig2, dn2) - E.deltaE(ssig1, csig1, dn1));
00504 if (scalep) {
00505 real csig12 = csig1 * csig2 + ssig1 * ssig2;
00506 real t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
00507 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
00508 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
00509 }
00510 }
00511
00512 Math::real GeodesicExact::Astroid(real x, real y) {
00513
00514
00515 real k;
00516 real
00517 p = Math::sq(x),
00518 q = Math::sq(y),
00519 r = (p + q - 1) / 6;
00520 if ( !(q == 0 && r <= 0) ) {
00521 real
00522
00523
00524 S = p * q / 4,
00525 r2 = Math::sq(r),
00526 r3 = r * r2,
00527
00528
00529 disc = S * (S + 2 * r3);
00530 real u = r;
00531 if (disc >= 0) {
00532 real T3 = S + r3;
00533
00534
00535
00536 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
00537
00538 real T = Math::cbrt(T3);
00539
00540 u += T + (T ? r2 / T : 0);
00541 } else {
00542
00543 real ang = atan2(sqrt(-disc), -(S + r3));
00544
00545
00546 u += 2 * r * cos(ang / 3);
00547 }
00548 real
00549 v = sqrt(Math::sq(u) + q),
00550
00551 uv = u < 0 ? q / (v - u) : u + v,
00552 w = (uv - q) / (2 * v);
00553
00554
00555 k = uv / (sqrt(uv + Math::sq(w)) + w);
00556 } else {
00557
00558
00559 k = 0;
00560 }
00561 return k;
00562 }
00563
00564 Math::real GeodesicExact::InverseStart(EllipticFunction& E,
00565 real sbet1, real cbet1, real dn1,
00566 real sbet2, real cbet2, real dn2,
00567 real lam12,
00568 real& salp1, real& calp1,
00569
00570 real& salp2, real& calp2,
00571
00572 real& dnm)
00573 const {
00574
00575
00576
00577 real
00578 sig12 = -1,
00579
00580 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
00581 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
00582 #if defined(__GNUC__) && __GNUC__ == 4 && \
00583 (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
00584
00585
00586
00587
00588
00589
00590 real sbet12a;
00591 {
00592 GEOGRAPHICLIB_VOLATILE real xx1 = sbet2 * cbet1;
00593 GEOGRAPHICLIB_VOLATILE real xx2 = cbet2 * sbet1;
00594 sbet12a = xx1 + xx2;
00595 }
00596 #else
00597 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
00598 #endif
00599 bool shortline = cbet12 >= 0 && sbet12 < real(0.5) &&
00600 cbet2 * lam12 < real(0.5);
00601 real omg12 = lam12;
00602 if (shortline) {
00603 real sbetm2 = Math::sq(sbet1 + sbet2);
00604
00605
00606 sbetm2 /= sbetm2 + Math::sq(cbet1 + cbet2);
00607 dnm = sqrt(1 + _ep2 * sbetm2);
00608 omg12 /= _f1 * dnm;
00609 }
00610 real somg12 = sin(omg12), comg12 = cos(omg12);
00611
00612 salp1 = cbet2 * somg12;
00613 calp1 = comg12 >= 0 ?
00614 sbet12 + cbet2 * sbet1 * Math::sq(somg12) / (1 + comg12) :
00615 sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
00616
00617 real
00618 ssig12 = Math::hypot(salp1, calp1),
00619 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
00620
00621 if (shortline && ssig12 < _etol2) {
00622
00623 salp2 = cbet1 * somg12;
00624 calp2 = sbet12 - cbet1 * sbet2 *
00625 (comg12 >= 0 ? Math::sq(somg12) / (1 + comg12) : 1 - comg12);
00626 SinCosNorm(salp2, calp2);
00627
00628 sig12 = atan2(ssig12, csig12);
00629 } else if (abs(_n) > real(0.1) ||
00630 csig12 >= 0 ||
00631 ssig12 >= 6 * abs(_n) * Math::pi() * Math::sq(cbet1)) {
00632
00633 } else {
00634
00635
00636 real y, lamscale, betscale;
00637
00638
00639
00640 GEOGRAPHICLIB_VOLATILE real x;
00641 if (_f >= 0) {
00642
00643 {
00644 real k2 = Math::sq(sbet1) * _ep2;
00645 E.Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
00646 lamscale = _e2/_f1 * cbet1 * 2 * E.H();
00647 }
00648 betscale = lamscale * cbet1;
00649
00650 x = (lam12 - Math::pi()) / lamscale;
00651 y = sbet12a / betscale;
00652 } else {
00653
00654 real
00655 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
00656 bet12a = atan2(sbet12a, cbet12a);
00657 real m12b, m0, dummy;
00658
00659
00660 Lengths(E, Math::pi() + bet12a,
00661 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
00662 cbet1, cbet2, dummy, m12b, m0, false,
00663 dummy, dummy);
00664 x = -1 + m12b / (cbet1 * cbet2 * m0 * Math::pi());
00665 betscale = x < -real(0.01) ? sbet12a / x :
00666 -_f * Math::sq(cbet1) * Math::pi();
00667 lamscale = betscale / cbet1;
00668 y = (lam12 - Math::pi()) / lamscale;
00669 }
00670
00671 if (y > -tol1_ && x > -1 - xthresh_) {
00672
00673
00674 if (_f >= 0) {
00675 salp1 = min(real(1), -real(x)); calp1 = - sqrt(1 - Math::sq(salp1));
00676 } else {
00677 calp1 = max(real(x > -tol1_ ? 0 : -1), real(x));
00678 salp1 = sqrt(1 - Math::sq(calp1));
00679 }
00680 } else {
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715 real k = Astroid(x, y);
00716 real
00717 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
00718 somg12 = sin(omg12a); comg12 = -cos(omg12a);
00719
00720 salp1 = cbet2 * somg12;
00721 calp1 = sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
00722 }
00723 }
00724 if (salp1 > 0)
00725 SinCosNorm(salp1, calp1);
00726 else {
00727 salp1 = 1; calp1 = 0;
00728 }
00729 return sig12;
00730 }
00731
00732 Math::real GeodesicExact::Lambda12(real sbet1, real cbet1, real dn1,
00733 real sbet2, real cbet2, real dn2,
00734 real salp1, real calp1,
00735 real& salp2, real& calp2,
00736 real& sig12,
00737 real& ssig1, real& csig1,
00738 real& ssig2, real& csig2,
00739 EllipticFunction& E,
00740 real& omg12,
00741 bool diffp, real& dlam12) const
00742 {
00743
00744 if (sbet1 == 0 && calp1 == 0)
00745
00746
00747 calp1 = -tiny_;
00748
00749 real
00750
00751 salp0 = salp1 * cbet1,
00752 calp0 = Math::hypot(calp1, salp1 * sbet1);
00753
00754 real somg1, comg1, somg2, comg2, cchi1, cchi2, lam12;
00755
00756
00757 ssig1 = sbet1; somg1 = salp0 * sbet1;
00758 csig1 = comg1 = calp1 * cbet1;
00759
00760 cchi1 = _f1 * dn1 * comg1;
00761 SinCosNorm(ssig1, csig1);
00762
00763
00764
00765
00766
00767
00768
00769 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
00770
00771
00772
00773
00774 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ?
00775 sqrt(Math::sq(calp1 * cbet1) +
00776 (cbet1 < -sbet1 ?
00777 (cbet2 - cbet1) * (cbet1 + cbet2) :
00778 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
00779 abs(calp1);
00780
00781
00782 ssig2 = sbet2; somg2 = salp0 * sbet2;
00783 csig2 = comg2 = calp2 * cbet2;
00784
00785 cchi2 = _f1 * dn2 * comg2;
00786 SinCosNorm(ssig2, csig2);
00787
00788
00789
00790
00791 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
00792 csig1 * csig2 + ssig1 * ssig2);
00793
00794
00795 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, real(0)),
00796 comg1 * comg2 + somg1 * somg2);
00797 real k2 = Math::sq(calp0) * _ep2;
00798 E.Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
00799 real chi12 = atan2(max(cchi1 * somg2 - somg1 * cchi2, real(0)),
00800 cchi1 * cchi2 + somg1 * somg2);
00801 lam12 = chi12 -
00802 _e2/_f1 * salp0 * E.H() / (Math::pi() / 2) *
00803 (sig12 + E.deltaH(ssig2, csig2, dn2) - E.deltaH(ssig1, csig1, dn1) );
00804
00805 if (diffp) {
00806 if (calp2 == 0)
00807 dlam12 = - 2 * _f1 * dn1 / sbet1;
00808 else {
00809 real dummy;
00810 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
00811 cbet1, cbet2, dummy, dlam12, dummy,
00812 false, dummy, dummy);
00813 dlam12 *= _f1 / (calp2 * cbet2);
00814 }
00815 }
00816
00817 return lam12;
00818 }
00819
00820 void GeodesicExact::C4f(real eps, real c[]) const {
00821
00822
00823 for (int j = nC4x_, k = nC4_; k; ) {
00824 real t = 0;
00825 for (int i = nC4_ - k + 1; i; --i)
00826 t = eps * t + _C4x[--j];
00827 c[--k] = t;
00828 }
00829
00830 real mult = 1;
00831 for (int k = 1; k < nC4_; ) {
00832 mult *= eps;
00833 c[k++] *= mult;
00834 }
00835 }
00836
00837
00838
00839
00840
00841
00842 void GeodesicExact::C4coeff() {
00843 const real* cc = rawC4coeff();
00844
00845 for (int m = 0, k = 0, h = 0; m < nC4_; ++m) {
00846
00847 for (int j = m; j < nC4_; ++j) {
00848 real t = 0;
00849
00850 for (int l = nC4_ - j; l--;)
00851 t = _n * t + cc[h++];
00852 _C4x[k++] = t/cc[h++];
00853 }
00854 }
00855 return;
00856 }
00857
00858 }