GeographicLib
1.21
|
00001 /** 00002 * \file Geodesic.cpp 00003 * \brief Implementation for GeographicLib::Geodesic class 00004 * 00005 * Copyright (c) Charles Karney (2009-2012) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 * 00009 * This is a reformulation of the geodesic problem. The notation is as 00010 * follows: 00011 * - at a general point (no suffix or 1 or 2 as suffix) 00012 * - phi = latitude 00013 * - beta = latitude on auxiliary sphere 00014 * - omega = longitude on auxiliary sphere 00015 * - lambda = longitude 00016 * - alpha = azimuth of great circle 00017 * - sigma = arc length along great circle 00018 * - s = distance 00019 * - tau = scaled distance (= sigma at multiples of pi/2) 00020 * - at northwards equator crossing 00021 * - beta = phi = 0 00022 * - omega = lambda = 0 00023 * - alpha = alpha0 00024 * - sigma = s = 0 00025 * - a 12 suffix means a difference, e.g., s12 = s2 - s1. 00026 * - s and c prefixes mean sin and cos 00027 **********************************************************************/ 00028 00029 #include <GeographicLib/Geodesic.hpp> 00030 #include <GeographicLib/GeodesicLine.hpp> 00031 00032 #define GEOGRAPHICLIB_GEODESIC_CPP \ 00033 "$Id: dd137806b8a5ba58211a37eb87e163b8a9bd7aa7 $" 00034 00035 RCSID_DECL(GEOGRAPHICLIB_GEODESIC_CPP) 00036 RCSID_DECL(GEOGRAPHICLIB_GEODESIC_HPP) 00037 00038 namespace GeographicLib { 00039 00040 using namespace std; 00041 00042 // Underflow guard. We require 00043 // tiny_ * epsilon() > 0 00044 // tiny_ + epsilon() == epsilon() 00045 const Math::real Geodesic::tiny_ = sqrt(numeric_limits<real>::min()); 00046 const Math::real Geodesic::tol0_ = numeric_limits<real>::epsilon(); 00047 // Increase multiplier in defn of tol1_ from 100 to 200 to fix inverse case 00048 // 52.784459512564 0 -52.784459512563990912 179.634407464943777557 00049 // which otherwise failed for Visual Studio 10 (Release and Debug) 00050 const Math::real Geodesic::tol1_ = 200 * tol0_; 00051 const Math::real Geodesic::tol2_ = sqrt(numeric_limits<real>::epsilon()); 00052 const Math::real Geodesic::xthresh_ = 1000 * tol2_; 00053 00054 Geodesic::Geodesic(real a, real f) 00055 : _a(a) 00056 , _f(f <= 1 ? f : 1/f) 00057 , _f1(1 - _f) 00058 , _e2(_f * (2 - _f)) 00059 , _ep2(_e2 / Math::sq(_f1)) // e2 / (1 - e2) 00060 , _n(_f / ( 2 - _f)) 00061 , _b(_a * _f1) 00062 , _c2((Math::sq(_a) + Math::sq(_b) * 00063 (_e2 == 0 ? 1 : 00064 (_e2 > 0 ? Math::atanh(sqrt(_e2)) : atan(sqrt(-_e2))) / 00065 sqrt(abs(_e2))))/2) // authalic radius squared 00066 // The sig12 threshold for "really short" 00067 , _etol2(10 * tol2_ / max(real(0.1), sqrt(abs(_e2)))) 00068 { 00069 if (!(Math::isfinite(_a) && _a > 0)) 00070 throw GeographicErr("Major radius is not positive"); 00071 if (!(Math::isfinite(_b) && _b > 0)) 00072 throw GeographicErr("Minor radius is not positive"); 00073 A3coeff(); 00074 C3coeff(); 00075 C4coeff(); 00076 } 00077 00078 const Geodesic Geodesic::WGS84(Constants::WGS84_a<real>(), 00079 Constants::WGS84_f<real>()); 00080 00081 Math::real Geodesic::SinCosSeries(bool sinp, 00082 real sinx, real cosx, 00083 const real c[], int n) throw() { 00084 // Evaluate 00085 // y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) : 00086 // sum(c[i] * cos((2*i+1) * x), i, 0, n-1) : 00087 // using Clenshaw summation. N.B. c[0] is unused for sin series 00088 // Approx operation count = (n + 5) mult and (2 * n + 2) add 00089 c += (n + sinp); // Point to one beyond last element 00090 real 00091 ar = 2 * (cosx - sinx) * (cosx + sinx), // 2 * cos(2 * x) 00092 y0 = n & 1 ? *--c : 0, y1 = 0; // accumulators for sum 00093 // Now n is even 00094 n /= 2; 00095 while (n--) { 00096 // Unroll loop x 2, so accumulators return to their original role 00097 y1 = ar * y0 - y1 + *--c; 00098 y0 = ar * y1 - y0 + *--c; 00099 } 00100 return sinp 00101 ? 2 * sinx * cosx * y0 // sin(2 * x) * y0 00102 : cosx * (y0 - y1); // cos(x) * (y0 - y1) 00103 } 00104 00105 GeodesicLine Geodesic::Line(real lat1, real lon1, real azi1, unsigned caps) 00106 const throw() { 00107 return GeodesicLine(*this, lat1, lon1, azi1, caps); 00108 } 00109 00110 Math::real Geodesic::GenDirect(real lat1, real lon1, real azi1, 00111 bool arcmode, real s12_a12, unsigned outmask, 00112 real& lat2, real& lon2, real& azi2, 00113 real& s12, real& m12, real& M12, real& M21, 00114 real& S12) const throw() { 00115 return GeodesicLine(*this, lat1, lon1, azi1, 00116 // Automatically supply DISTANCE_IN if necessary 00117 outmask | (arcmode ? NONE : DISTANCE_IN)) 00118 . // Note the dot! 00119 GenPosition(arcmode, s12_a12, outmask, 00120 lat2, lon2, azi2, s12, m12, M12, M21, S12); 00121 } 00122 00123 Math::real Geodesic::GenInverse(real lat1, real lon1, real lat2, real lon2, 00124 unsigned outmask, 00125 real& s12, real& azi1, real& azi2, 00126 real& m12, real& M12, real& M21, real& S12) 00127 const throw() { 00128 outmask &= OUT_ALL; 00129 lon1 = AngNormalize(lon1); 00130 real lon12 = AngNormalize(AngNormalize(lon2) - lon1); 00131 // If very close to being on the same meridian, then make it so. 00132 // Not sure this is necessary... 00133 lon12 = AngRound(lon12); 00134 // Make longitude difference positive. 00135 int lonsign = lon12 >= 0 ? 1 : -1; 00136 lon12 *= lonsign; 00137 if (lon12 == 180) 00138 lonsign = 1; 00139 // If really close to the equator, treat as on equator. 00140 lat1 = AngRound(lat1); 00141 lat2 = AngRound(lat2); 00142 // Swap points so that point with higher (abs) latitude is point 1 00143 int swapp = abs(lat1) >= abs(lat2) ? 1 : -1; 00144 if (swapp < 0) { 00145 lonsign *= -1; 00146 swap(lat1, lat2); 00147 } 00148 // Make lat1 <= 0 00149 int latsign = lat1 < 0 ? 1 : -1; 00150 lat1 *= latsign; 00151 lat2 *= latsign; 00152 // Now we have 00153 // 00154 // 0 <= lon12 <= 180 00155 // -90 <= lat1 <= 0 00156 // lat1 <= lat2 <= -lat1 00157 // 00158 // longsign, swapp, latsign register the transformation to bring the 00159 // coordinates to this canonical form. In all cases, 1 means no change was 00160 // made. We make these transformations so that there are few cases to 00161 // check, e.g., on verifying quadrants in atan2. In addition, this 00162 // enforces some symmetries in the results returned. 00163 00164 real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x; 00165 00166 phi = lat1 * Math::degree<real>(); 00167 // Ensure cbet1 = +epsilon at poles 00168 sbet1 = _f1 * sin(phi); 00169 cbet1 = lat1 == -90 ? tiny_ : cos(phi); 00170 SinCosNorm(sbet1, cbet1); 00171 00172 phi = lat2 * Math::degree<real>(); 00173 // Ensure cbet2 = +epsilon at poles 00174 sbet2 = _f1 * sin(phi); 00175 cbet2 = abs(lat2) == 90 ? tiny_ : cos(phi); 00176 SinCosNorm(sbet2, cbet2); 00177 00178 // If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the 00179 // |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is 00180 // a better measure. This logic is used in assigning calp2 in Lambda12. 00181 // Sometimes these quantities vanish and in that case we force bet2 = +/- 00182 // bet1 exactly. An example where is is necessary is the inverse problem 00183 // 48.522876735459 0 -48.52287673545898293 179.599720456223079643 00184 // which failed with Visual Studio 10 (Release and Debug) 00185 00186 if (cbet1 < -sbet1) { 00187 if (cbet2 == cbet1) 00188 sbet2 = sbet2 < 0 ? sbet1 : -sbet1; 00189 } else { 00190 if (abs(sbet2) == -sbet1) 00191 cbet2 = cbet1; 00192 } 00193 00194 real 00195 lam12 = lon12 * Math::degree<real>(), 00196 slam12 = lon12 == 180 ? 0 : sin(lam12), 00197 clam12 = cos(lam12); // lon12 == 90 isn't interesting 00198 00199 real a12, sig12, calp1, salp1, calp2, salp2; 00200 // index zero elements of these arrays are unused 00201 real C1a[nC1_ + 1], C2a[nC2_ + 1], C3a[nC3_]; 00202 00203 bool meridian = lat1 == -90 || slam12 == 0; 00204 00205 if (meridian) { 00206 00207 // Endpoints are on a single full meridian, so the geodesic might lie on 00208 // a meridian. 00209 00210 calp1 = clam12; salp1 = slam12; // Head to the target longitude 00211 calp2 = 1; salp2 = 0; // At the target we're heading north 00212 00213 real 00214 // tan(bet) = tan(sig) * cos(alp) 00215 ssig1 = sbet1, csig1 = calp1 * cbet1, 00216 ssig2 = sbet2, csig2 = calp2 * cbet2; 00217 00218 // sig12 = sig2 - sig1 00219 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)), 00220 csig1 * csig2 + ssig1 * ssig2); 00221 { 00222 real dummy; 00223 Lengths(_n, sig12, ssig1, csig1, ssig2, csig2, 00224 cbet1, cbet2, s12x, m12x, dummy, 00225 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a); 00226 } 00227 // Add the check for sig12 since zero length geodesics might yield m12 < 00228 // 0. Test case was 00229 // 00230 // echo 20.001 0 20.001 0 | Geod -i 00231 // 00232 // In fact, we will have sig12 > pi/2 for meridional geodesic which is 00233 // not a shortest path. 00234 if (sig12 < 1 || m12x >= 0) { 00235 m12x *= _a; 00236 s12x *= _b; 00237 a12 = sig12 / Math::degree<real>(); 00238 } else 00239 // m12 < 0, i.e., prolate and too close to anti-podal 00240 meridian = false; 00241 } 00242 00243 real omg12; 00244 if (!meridian && 00245 sbet1 == 0 && // and sbet2 == 0 00246 // Mimic the way Lambda12 works with calp1 = 0 00247 (_f <= 0 || lam12 <= Math::pi<real>() - _f * Math::pi<real>())) { 00248 00249 // Geodesic runs along equator 00250 calp1 = calp2 = 0; salp1 = salp2 = 1; 00251 s12x = _a * lam12; 00252 m12x = _b * sin(lam12 / _f1); 00253 if (outmask & GEODESICSCALE) 00254 M12 = M21 = cos(lam12 / _f1); 00255 a12 = lon12 / _f1; 00256 sig12 = omg12 = lam12 / _f1; 00257 00258 } else if (!meridian) { 00259 00260 // Now point1 and point2 belong within a hemisphere bounded by a 00261 // meridian and geodesic is neither meridional or equatorial. 00262 00263 // Figure a starting point for Newton's method 00264 sig12 = InverseStart(sbet1, cbet1, sbet2, cbet2, 00265 lam12, 00266 salp1, calp1, salp2, calp2, 00267 C1a, C2a); 00268 00269 if (sig12 >= 0) { 00270 // Short lines (InverseStart sets salp2, calp2) 00271 real wm = sqrt(1 - _e2 * Math::sq((cbet1 + cbet2) / 2)); 00272 s12x = sig12 * _a * wm; 00273 m12x = Math::sq(wm) * _a / _f1 * sin(sig12 * _f1 / wm); 00274 if (outmask & GEODESICSCALE) 00275 M12 = M21 = cos(sig12 * _f1 / wm); 00276 a12 = sig12 / Math::degree<real>(); 00277 omg12 = lam12 / wm; 00278 } else { 00279 00280 // Newton's method 00281 real ssig1, csig1, ssig2, csig2, eps; 00282 real ov = 0; 00283 unsigned numit = 0; 00284 for (unsigned trip = 0; numit < maxit_; ++numit) { 00285 real dv; 00286 real v = Lambda12(sbet1, cbet1, sbet2, cbet2, salp1, calp1, 00287 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, 00288 eps, omg12, trip < 1, dv, C1a, C2a, C3a) - lam12; 00289 if (!(abs(v) > tiny_) || !(trip < 1)) { 00290 if (!(abs(v) <= max(tol1_, ov))) 00291 numit = maxit_; 00292 break; 00293 } 00294 real 00295 dalp1 = -v/dv; 00296 real 00297 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1), 00298 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1; 00299 calp1 = calp1 * cdalp1 - salp1 * sdalp1; 00300 salp1 = max(real(0), nsalp1); 00301 SinCosNorm(salp1, calp1); 00302 // In some regimes we don't get quadratic convergence because slope 00303 // -> 0. So use convergence conditions based on epsilon instead of 00304 // sqrt(epsilon). The first criterion is a test on abs(v) against 00305 // 100 * epsilon. The second takes credit for an anticipated 00306 // reduction in abs(v) by v/ov (due to the latest update in alp1) and 00307 // checks this against epsilon. 00308 if (!(abs(v) >= tol1_ && Math::sq(v) >= ov * tol0_)) ++trip; 00309 ov = abs(v); 00310 } 00311 00312 if (numit >= maxit_) { 00313 // Signal failure. 00314 if (outmask & DISTANCE) 00315 s12 = Math::NaN<real>(); 00316 if (outmask & AZIMUTH) 00317 azi1 = azi2 = Math::NaN<real>(); 00318 if (outmask & REDUCEDLENGTH) 00319 m12 = Math::NaN<real>(); 00320 if (outmask & GEODESICSCALE) 00321 M12 = M21 = Math::NaN<real>(); 00322 if (outmask & AREA) 00323 S12 = Math::NaN<real>(); 00324 return Math::NaN<real>(); 00325 } 00326 00327 { 00328 real dummy; 00329 Lengths(eps, sig12, ssig1, csig1, ssig2, csig2, 00330 cbet1, cbet2, s12x, m12x, dummy, 00331 (outmask & GEODESICSCALE) != 0U, M12, M21, C1a, C2a); 00332 } 00333 m12x *= _a; 00334 s12x *= _b; 00335 a12 = sig12 / Math::degree<real>(); 00336 omg12 = lam12 - omg12; 00337 } 00338 } 00339 00340 if (outmask & DISTANCE) 00341 s12 = 0 + s12x; // Convert -0 to 0 00342 00343 if (outmask & REDUCEDLENGTH) 00344 m12 = 0 + m12x; // Convert -0 to 0 00345 00346 if (outmask & AREA) { 00347 real 00348 // From Lambda12: sin(alp1) * cos(bet1) = sin(alp0) 00349 salp0 = salp1 * cbet1, 00350 calp0 = Math::hypot(calp1, salp1 * sbet1); // calp0 > 0 00351 real alp12; 00352 if (calp0 != 0 && salp0 != 0) { 00353 real 00354 // From Lambda12: tan(bet) = tan(sig) * cos(alp) 00355 ssig1 = sbet1, csig1 = calp1 * cbet1, 00356 ssig2 = sbet2, csig2 = calp2 * cbet2, 00357 k2 = Math::sq(calp0) * _ep2, 00358 // Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0). 00359 A4 = Math::sq(_a) * calp0 * salp0 * _e2; 00360 SinCosNorm(ssig1, csig1); 00361 SinCosNorm(ssig2, csig2); 00362 real C4a[nC4_]; 00363 C4f(k2, C4a); 00364 real 00365 B41 = SinCosSeries(false, ssig1, csig1, C4a, nC4_), 00366 B42 = SinCosSeries(false, ssig2, csig2, C4a, nC4_); 00367 S12 = A4 * (B42 - B41); 00368 } else 00369 // Avoid problems with indeterminate sig1, sig2 on equator 00370 S12 = 0; 00371 00372 if (!meridian && 00373 omg12 < real(0.75) * Math::pi<real>() && // Long difference too big 00374 sbet2 - sbet1 < real(1.75)) { // Lat difference too big 00375 // Use tan(Gamma/2) = tan(omg12/2) 00376 // * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2)) 00377 // with tan(x/2) = sin(x)/(1+cos(x)) 00378 real 00379 somg12 = sin(omg12), domg12 = 1 + cos(omg12), 00380 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2; 00381 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ), 00382 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) ); 00383 } else { 00384 // alp12 = alp2 - alp1, used in atan2 so no need to normalize 00385 real 00386 salp12 = salp2 * calp1 - calp2 * salp1, 00387 calp12 = calp2 * calp1 + salp2 * salp1; 00388 // The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz 00389 // salp12 = -0 and alp12 = -180. However this depends on the sign 00390 // being attached to 0 correctly. The following ensures the correct 00391 // behavior. 00392 if (salp12 == 0 && calp12 < 0) { 00393 salp12 = tiny_ * calp1; 00394 calp12 = -1; 00395 } 00396 alp12 = atan2(salp12, calp12); 00397 } 00398 S12 += _c2 * alp12; 00399 S12 *= swapp * lonsign * latsign; 00400 // Convert -0 to 0 00401 S12 += 0; 00402 } 00403 00404 // Convert calp, salp to azimuth accounting for lonsign, swapp, latsign. 00405 if (swapp < 0) { 00406 swap(salp1, salp2); 00407 swap(calp1, calp2); 00408 if (outmask & GEODESICSCALE) 00409 swap(M12, M21); 00410 } 00411 00412 salp1 *= swapp * lonsign; calp1 *= swapp * latsign; 00413 salp2 *= swapp * lonsign; calp2 *= swapp * latsign; 00414 00415 if (outmask & AZIMUTH) { 00416 // minus signs give range [-180, 180). 0- converts -0 to +0. 00417 azi1 = 0 - atan2(-salp1, calp1) / Math::degree<real>(); 00418 azi2 = 0 - atan2(-salp2, calp2) / Math::degree<real>(); 00419 } 00420 00421 // Returned value in [0, 180] 00422 return a12; 00423 } 00424 00425 void Geodesic::Lengths(real eps, real sig12, 00426 real ssig1, real csig1, real ssig2, real csig2, 00427 real cbet1, real cbet2, 00428 real& s12b, real& m12a, real& m0, 00429 bool scalep, real& M12, real& M21, 00430 // Scratch areas of the right size 00431 real C1a[], real C2a[]) const throw() { 00432 // Return m12a = (reduced length)/_a; also calculate s12b = distance/_b, 00433 // and m0 = coefficient of secular term in expression for reduced length. 00434 C1f(eps, C1a); 00435 C2f(eps, C2a); 00436 real 00437 A1m1 = A1m1f(eps), 00438 AB1 = (1 + A1m1) * (SinCosSeries(true, ssig2, csig2, C1a, nC1_) - 00439 SinCosSeries(true, ssig1, csig1, C1a, nC1_)), 00440 A2m1 = A2m1f(eps), 00441 AB2 = (1 + A2m1) * (SinCosSeries(true, ssig2, csig2, C2a, nC2_) - 00442 SinCosSeries(true, ssig1, csig1, C2a, nC2_)), 00443 cbet1sq = Math::sq(cbet1), 00444 cbet2sq = Math::sq(cbet2), 00445 w1 = sqrt(1 - _e2 * cbet1sq), 00446 w2 = sqrt(1 - _e2 * cbet2sq), 00447 // Make sure it's OK to have repeated dummy arguments 00448 m0x = A1m1 - A2m1, 00449 J12 = m0x * sig12 + (AB1 - AB2); 00450 m0 = m0x; 00451 // Missing a factor of _a. 00452 // Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure accurate 00453 // cancellation in the case of coincident points. 00454 m12a = (w2 * (csig1 * ssig2) - w1 * (ssig1 * csig2)) 00455 - _f1 * csig1 * csig2 * J12; 00456 // Missing a factor of _b 00457 s12b = (1 + A1m1) * sig12 + AB1; 00458 if (scalep) { 00459 real csig12 = csig1 * csig2 + ssig1 * ssig2; 00460 J12 *= _f1; 00461 M12 = csig12 + (_e2 * (cbet1sq - cbet2sq) * ssig2 / (w1 + w2) 00462 - csig2 * J12) * ssig1 / w1; 00463 M21 = csig12 - (_e2 * (cbet1sq - cbet2sq) * ssig1 / (w1 + w2) 00464 - csig1 * J12) * ssig2 / w2; 00465 } 00466 } 00467 00468 Math::real Geodesic::Astroid(real x, real y) throw() { 00469 // Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k. 00470 // This solution is adapted from Geocentric::Reverse. 00471 real k; 00472 real 00473 p = Math::sq(x), 00474 q = Math::sq(y), 00475 r = (p + q - 1) / 6; 00476 if ( !(q == 0 && r <= 0) ) { 00477 real 00478 // Avoid possible division by zero when r = 0 by multiplying equations 00479 // for s and t by r^3 and r, resp. 00480 S = p * q / 4, // S = r^3 * s 00481 r2 = Math::sq(r), 00482 r3 = r * r2, 00483 // The discrimant of the quadratic equation for T3. This is zero on 00484 // the evolute curve p^(1/3)+q^(1/3) = 1 00485 disc = S * (S + 2 * r3); 00486 real u = r; 00487 if (disc >= 0) { 00488 real T3 = S + r3; 00489 // Pick the sign on the sqrt to maximize abs(T3). This minimizes loss 00490 // of precision due to cancellation. The result is unchanged because 00491 // of the way the T is used in definition of u. 00492 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc); // T3 = (r * t)^3 00493 // N.B. cbrt always returns the real root. cbrt(-8) = -2. 00494 real T = Math::cbrt(T3); // T = r * t 00495 // T can be zero; but then r2 / T -> 0. 00496 u += T + (T != 0 ? r2 / T : 0); 00497 } else { 00498 // T is complex, but the way u is defined the result is real. 00499 real ang = atan2(sqrt(-disc), -(S + r3)); 00500 // There are three possible cube roots. We choose the root which 00501 // avoids cancellation. Note that disc < 0 implies that r < 0. 00502 u += 2 * r * cos(ang / 3); 00503 } 00504 real 00505 v = sqrt(Math::sq(u) + q), // guaranteed positive 00506 // Avoid loss of accuracy when u < 0. 00507 uv = u < 0 ? q / (v - u) : u + v, // u+v, guaranteed positive 00508 w = (uv - q) / (2 * v); // positive? 00509 // Rearrange expression for k to avoid loss of accuracy due to 00510 // subtraction. Division by 0 not possible because uv > 0, w >= 0. 00511 k = uv / (sqrt(uv + Math::sq(w)) + w); // guaranteed positive 00512 } else { // q == 0 && r <= 0 00513 // y = 0 with |x| <= 1. Handle this case directly. 00514 // for y small, positive root is k = abs(y)/sqrt(1-x^2) 00515 k = 0; 00516 } 00517 return k; 00518 } 00519 00520 Math::real Geodesic::InverseStart(real sbet1, real cbet1, 00521 real sbet2, real cbet2, 00522 real lam12, 00523 real& salp1, real& calp1, 00524 // Only updated if return val >= 0 00525 real& salp2, real& calp2, 00526 // Scratch areas of the right size 00527 real C1a[], real C2a[]) const throw() { 00528 // Return a starting point for Newton's method in salp1 and calp1 (function 00529 // value is -1). If Newton's method doesn't need to be used, return also 00530 // salp2 and calp2 and function value is sig12. 00531 real 00532 sig12 = -1, // Return value 00533 // bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0] 00534 sbet12 = sbet2 * cbet1 - cbet2 * sbet1, 00535 cbet12 = cbet2 * cbet1 + sbet2 * sbet1; 00536 #if defined(__GNUC__) && __GNUC__ == 4 && \ 00537 (__GNUC_MINOR__ < 6 || defined(__MINGW32__)) 00538 // Volatile declaration needed to fix inverse cases 00539 // 88.202499451857 0 -88.202499451857 179.981022032992859592 00540 // 89.262080389218 0 -89.262080389218 179.992207982775375662 00541 // 89.333123580033 0 -89.333123580032997687 179.99295812360148422 00542 // which otherwise fail with g++ 4.4.4 x86 -O3 (Linux) 00543 // and g++ 4.4.0 (mingw) and g++ 4.6.1 (tdm mingw). 00544 real sbet12a; 00545 { 00546 volatile real xx1 = sbet2 * cbet1; 00547 volatile real xx2 = cbet2 * sbet1; 00548 sbet12a = xx1 + xx2; 00549 } 00550 #else 00551 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1; 00552 #endif 00553 bool shortline = cbet12 >= 0 && sbet12 < real(0.5) && 00554 lam12 <= Math::pi<real>() / 6; 00555 real 00556 omg12 = (!shortline ? lam12 : 00557 lam12 / sqrt(1 - _e2 * Math::sq((cbet1 + cbet2) / 2))), 00558 somg12 = sin(omg12), comg12 = cos(omg12); 00559 00560 salp1 = cbet2 * somg12; 00561 calp1 = comg12 >= 0 ? 00562 sbet12 + cbet2 * sbet1 * Math::sq(somg12) / (1 + comg12) : 00563 sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12); 00564 00565 real 00566 ssig12 = Math::hypot(salp1, calp1), 00567 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12; 00568 00569 if (shortline && ssig12 < _etol2) { 00570 // really short lines 00571 salp2 = cbet1 * somg12; 00572 calp2 = sbet12 - cbet1 * sbet2 * Math::sq(somg12) / (1 + comg12); 00573 SinCosNorm(salp2, calp2); 00574 // Set return value 00575 sig12 = atan2(ssig12, csig12); 00576 } else if (csig12 >= 0 || 00577 ssig12 >= 3 * abs(_f) * Math::pi<real>() * Math::sq(cbet1)) { 00578 // Nothing to do, zeroth order spherical approximation is OK 00579 } else { 00580 // Scale lam12 and bet2 to x, y coordinate system where antipodal point 00581 // is at origin and singular point is at y = 0, x = -1. 00582 real y, lamscale, betscale; 00583 // Volatile declaration needed to fix inverse case 00584 // 56.320923501171 0 -56.320923501171 179.664747671772880215 00585 // which otherwise fails with g++ 4.4.4 x86 -O3 00586 volatile real x; 00587 if (_f >= 0) { // In fact f == 0 does not get here 00588 // x = dlong, y = dlat 00589 { 00590 real 00591 k2 = Math::sq(sbet1) * _ep2, 00592 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2); 00593 lamscale = _f * cbet1 * A3f(eps) * Math::pi<real>(); 00594 } 00595 betscale = lamscale * cbet1; 00596 00597 x = (lam12 - Math::pi<real>()) / lamscale; 00598 y = sbet12a / betscale; 00599 } else { // _f < 0 00600 // x = dlat, y = dlong 00601 real 00602 cbet12a = cbet2 * cbet1 - sbet2 * sbet1, 00603 bet12a = atan2(sbet12a, cbet12a); 00604 real m12a, m0, dummy; 00605 // In the case of lon12 = 180, this repeats a calculation made in 00606 // Inverse. 00607 Lengths(_n, Math::pi<real>() + bet12a, sbet1, -cbet1, sbet2, cbet2, 00608 cbet1, cbet2, dummy, m12a, m0, false, 00609 dummy, dummy, C1a, C2a); 00610 x = -1 + m12a/(_f1 * cbet1 * cbet2 * m0 * Math::pi<real>()); 00611 betscale = x < -real(0.01) ? sbet12a / x : 00612 -_f * Math::sq(cbet1) * Math::pi<real>(); 00613 lamscale = betscale / cbet1; 00614 y = (lam12 - Math::pi<real>()) / lamscale; 00615 } 00616 00617 if (y > -tol1_ && x > -1 - xthresh_) { 00618 // strip near cut 00619 if (_f >= 0) { 00620 salp1 = min(real(1), -real(x)); calp1 = - sqrt(1 - Math::sq(salp1)); 00621 } else { 00622 calp1 = max(real(x > -tol1_ ? 0 : -1), real(x)); 00623 salp1 = sqrt(1 - Math::sq(calp1)); 00624 } 00625 } else { 00626 // Estimate alp1, by solving the astroid problem. 00627 // 00628 // Could estimate alpha1 = theta + pi/2, directly, i.e., 00629 // calp1 = y/k; salp1 = -x/(1+k); for _f >= 0 00630 // calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check) 00631 // 00632 // However, it's better to estimate omg12 from astroid and use 00633 // spherical formula to compute alp1. This reduces the mean number of 00634 // Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12 00635 // (min 0 max 5). The changes in the number of iterations are as 00636 // follows: 00637 // 00638 // change percent 00639 // 1 5 00640 // 0 78 00641 // -1 16 00642 // -2 0.6 00643 // -3 0.04 00644 // -4 0.002 00645 // 00646 // The histogram of iterations is (m = number of iterations estimating 00647 // alp1 directly, n = number of iterations estimating via omg12, total 00648 // number of trials = 148605): 00649 // 00650 // iter m n 00651 // 0 148 186 00652 // 1 13046 13845 00653 // 2 93315 102225 00654 // 3 36189 32341 00655 // 4 5396 7 00656 // 5 455 1 00657 // 6 56 0 00658 // 00659 // Because omg12 is near pi, estimate work with omg12a = pi - omg12 00660 real k = Astroid(x, y); 00661 real 00662 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k ), 00663 somg12 = sin(omg12a), comg12 = -cos(omg12a); 00664 // Update spherical estimate of alp1 using omg12 instead of lam12 00665 salp1 = cbet2 * somg12; 00666 calp1 = sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12); 00667 } 00668 } 00669 SinCosNorm(salp1, calp1); 00670 return sig12; 00671 } 00672 00673 Math::real Geodesic::Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00674 real salp1, real calp1, 00675 real& salp2, real& calp2, 00676 real& sig12, 00677 real& ssig1, real& csig1, 00678 real& ssig2, real& csig2, 00679 real& eps, real& domg12, 00680 bool diffp, real& dlam12, 00681 // Scratch areas of the right size 00682 real C1a[], real C2a[], real C3a[]) const 00683 throw() { 00684 00685 if (sbet1 == 0 && calp1 == 0) 00686 // Break degeneracy of equatorial line. This case has already been 00687 // handled. 00688 calp1 = -tiny_; 00689 00690 real 00691 // sin(alp1) * cos(bet1) = sin(alp0) 00692 salp0 = salp1 * cbet1, 00693 calp0 = Math::hypot(calp1, salp1 * sbet1); // calp0 > 0 00694 00695 real somg1, comg1, somg2, comg2, omg12, lam12; 00696 // tan(bet1) = tan(sig1) * cos(alp1) 00697 // tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1) 00698 ssig1 = sbet1; somg1 = salp0 * sbet1; 00699 csig1 = comg1 = calp1 * cbet1; 00700 SinCosNorm(ssig1, csig1); 00701 // SinCosNorm(somg1, comg1); -- don't need to normalize! 00702 00703 // Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful 00704 // about this case, since this can yield singularities in the Newton 00705 // iteration. 00706 // sin(alp2) * cos(bet2) = sin(alp0) 00707 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1; 00708 // calp2 = sqrt(1 - sq(salp2)) 00709 // = sqrt(sq(calp0) - sq(sbet2)) / cbet2 00710 // and subst for calp0 and rearrange to give (choose positive sqrt 00711 // to give alp2 in [0, pi/2]). 00712 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ? 00713 sqrt(Math::sq(calp1 * cbet1) + 00714 (cbet1 < -sbet1 ? 00715 (cbet2 - cbet1) * (cbet1 + cbet2) : 00716 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 : 00717 abs(calp1); 00718 // tan(bet2) = tan(sig2) * cos(alp2) 00719 // tan(omg2) = sin(alp0) * tan(sig2). 00720 ssig2 = sbet2; somg2 = salp0 * sbet2; 00721 csig2 = comg2 = calp2 * cbet2; 00722 SinCosNorm(ssig2, csig2); 00723 // SinCosNorm(somg2, comg2); -- don't need to normalize! 00724 00725 // sig12 = sig2 - sig1, limit to [0, pi] 00726 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)), 00727 csig1 * csig2 + ssig1 * ssig2); 00728 00729 // omg12 = omg2 - omg1, limit to [0, pi] 00730 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, real(0)), 00731 comg1 * comg2 + somg1 * somg2); 00732 real B312, h0; 00733 real k2 = Math::sq(calp0) * _ep2; 00734 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2); 00735 C3f(eps, C3a); 00736 B312 = (SinCosSeries(true, ssig2, csig2, C3a, nC3_-1) - 00737 SinCosSeries(true, ssig1, csig1, C3a, nC3_-1)); 00738 h0 = -_f * A3f(eps); 00739 domg12 = salp0 * h0 * (sig12 + B312); 00740 lam12 = omg12 + domg12; 00741 00742 if (diffp) { 00743 if (calp2 == 0) 00744 dlam12 = - 2 * sqrt(1 - _e2 * Math::sq(cbet1)) / sbet1; 00745 else { 00746 real dummy; 00747 Lengths(eps, sig12, ssig1, csig1, ssig2, csig2, 00748 cbet1, cbet2, dummy, dlam12, dummy, 00749 false, dummy, dummy, C1a, C2a); 00750 dlam12 /= calp2 * cbet2; 00751 } 00752 } 00753 00754 return lam12; 00755 } 00756 00757 Math::real Geodesic::A3f(real eps) const throw() { 00758 // Evaluation sum(_A3c[k] * eps^k, k, 0, nA3x_-1) by Horner's method 00759 real v = 0; 00760 for (int i = nA3x_; i; ) 00761 v = eps * v + _A3x[--i]; 00762 return v; 00763 } 00764 00765 void Geodesic::C3f(real eps, real c[]) const throw() { 00766 // Evaluation C3 coeffs by Horner's method 00767 // Elements c[1] thru c[nC3_ - 1] are set 00768 for (int j = nC3x_, k = nC3_ - 1; k; ) { 00769 real t = 0; 00770 for (int i = nC3_ - k; i; --i) 00771 t = eps * t + _C3x[--j]; 00772 c[k--] = t; 00773 } 00774 00775 real mult = 1; 00776 for (int k = 1; k < nC3_; ) { 00777 mult *= eps; 00778 c[k++] *= mult; 00779 } 00780 } 00781 00782 void Geodesic::C4f(real k2, real c[]) const throw() { 00783 // Evaluation C4 coeffs by Horner's method 00784 // Elements c[0] thru c[nC4_ - 1] are set 00785 for (int j = nC4x_, k = nC4_; k; ) { 00786 real t = 0; 00787 for (int i = nC4_ - k + 1; i; --i) 00788 t = k2 * t + _C4x[--j]; 00789 c[--k] = t; 00790 } 00791 00792 real mult = 1; 00793 for (int k = 1; k < nC4_; ) { 00794 mult *= k2; 00795 c[k++] *= mult; 00796 } 00797 } 00798 00799 // Generated by Maxima on 2010-09-04 10:26:17-04:00 00800 00801 // The scale factor A1-1 = mean value of I1-1 00802 Math::real Geodesic::A1m1f(real eps) throw() { 00803 real 00804 eps2 = Math::sq(eps), 00805 t; 00806 switch (nA1_/2) { 00807 case 0: 00808 t = 0; 00809 break; 00810 case 1: 00811 t = eps2/4; 00812 break; 00813 case 2: 00814 t = eps2*(eps2+16)/64; 00815 break; 00816 case 3: 00817 t = eps2*(eps2*(eps2+4)+64)/256; 00818 break; 00819 case 4: 00820 t = eps2*(eps2*(eps2*(25*eps2+64)+256)+4096)/16384; 00821 break; 00822 default: 00823 STATIC_ASSERT(nA1_ >= 0 && nA1_ <= 8, "Bad value of nA1_"); 00824 t = 0; 00825 } 00826 return (t + eps) / (1 - eps); 00827 } 00828 00829 // The coefficients C1[l] in the Fourier expansion of B1 00830 void Geodesic::C1f(real eps, real c[]) throw() { 00831 real 00832 eps2 = Math::sq(eps), 00833 d = eps; 00834 switch (nC1_) { 00835 case 0: 00836 break; 00837 case 1: 00838 c[1] = -d/2; 00839 break; 00840 case 2: 00841 c[1] = -d/2; 00842 d *= eps; 00843 c[2] = -d/16; 00844 break; 00845 case 3: 00846 c[1] = d*(3*eps2-8)/16; 00847 d *= eps; 00848 c[2] = -d/16; 00849 d *= eps; 00850 c[3] = -d/48; 00851 break; 00852 case 4: 00853 c[1] = d*(3*eps2-8)/16; 00854 d *= eps; 00855 c[2] = d*(eps2-2)/32; 00856 d *= eps; 00857 c[3] = -d/48; 00858 d *= eps; 00859 c[4] = -5*d/512; 00860 break; 00861 case 5: 00862 c[1] = d*((6-eps2)*eps2-16)/32; 00863 d *= eps; 00864 c[2] = d*(eps2-2)/32; 00865 d *= eps; 00866 c[3] = d*(9*eps2-16)/768; 00867 d *= eps; 00868 c[4] = -5*d/512; 00869 d *= eps; 00870 c[5] = -7*d/1280; 00871 break; 00872 case 6: 00873 c[1] = d*((6-eps2)*eps2-16)/32; 00874 d *= eps; 00875 c[2] = d*((64-9*eps2)*eps2-128)/2048; 00876 d *= eps; 00877 c[3] = d*(9*eps2-16)/768; 00878 d *= eps; 00879 c[4] = d*(3*eps2-5)/512; 00880 d *= eps; 00881 c[5] = -7*d/1280; 00882 d *= eps; 00883 c[6] = -7*d/2048; 00884 break; 00885 case 7: 00886 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048; 00887 d *= eps; 00888 c[2] = d*((64-9*eps2)*eps2-128)/2048; 00889 d *= eps; 00890 c[3] = d*((72-9*eps2)*eps2-128)/6144; 00891 d *= eps; 00892 c[4] = d*(3*eps2-5)/512; 00893 d *= eps; 00894 c[5] = d*(35*eps2-56)/10240; 00895 d *= eps; 00896 c[6] = -7*d/2048; 00897 d *= eps; 00898 c[7] = -33*d/14336; 00899 break; 00900 case 8: 00901 c[1] = d*(eps2*(eps2*(19*eps2-64)+384)-1024)/2048; 00902 d *= eps; 00903 c[2] = d*(eps2*(eps2*(7*eps2-18)+128)-256)/4096; 00904 d *= eps; 00905 c[3] = d*((72-9*eps2)*eps2-128)/6144; 00906 d *= eps; 00907 c[4] = d*((96-11*eps2)*eps2-160)/16384; 00908 d *= eps; 00909 c[5] = d*(35*eps2-56)/10240; 00910 d *= eps; 00911 c[6] = d*(9*eps2-14)/4096; 00912 d *= eps; 00913 c[7] = -33*d/14336; 00914 d *= eps; 00915 c[8] = -429*d/262144; 00916 break; 00917 default: 00918 STATIC_ASSERT(nC1_ >= 0 && nC1_ <= 8, "Bad value of nC1_"); 00919 } 00920 } 00921 00922 // The coefficients C1p[l] in the Fourier expansion of B1p 00923 void Geodesic::C1pf(real eps, real c[]) throw() { 00924 real 00925 eps2 = Math::sq(eps), 00926 d = eps; 00927 switch (nC1p_) { 00928 case 0: 00929 break; 00930 case 1: 00931 c[1] = d/2; 00932 break; 00933 case 2: 00934 c[1] = d/2; 00935 d *= eps; 00936 c[2] = 5*d/16; 00937 break; 00938 case 3: 00939 c[1] = d*(16-9*eps2)/32; 00940 d *= eps; 00941 c[2] = 5*d/16; 00942 d *= eps; 00943 c[3] = 29*d/96; 00944 break; 00945 case 4: 00946 c[1] = d*(16-9*eps2)/32; 00947 d *= eps; 00948 c[2] = d*(30-37*eps2)/96; 00949 d *= eps; 00950 c[3] = 29*d/96; 00951 d *= eps; 00952 c[4] = 539*d/1536; 00953 break; 00954 case 5: 00955 c[1] = d*(eps2*(205*eps2-432)+768)/1536; 00956 d *= eps; 00957 c[2] = d*(30-37*eps2)/96; 00958 d *= eps; 00959 c[3] = d*(116-225*eps2)/384; 00960 d *= eps; 00961 c[4] = 539*d/1536; 00962 d *= eps; 00963 c[5] = 3467*d/7680; 00964 break; 00965 case 6: 00966 c[1] = d*(eps2*(205*eps2-432)+768)/1536; 00967 d *= eps; 00968 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288; 00969 d *= eps; 00970 c[3] = d*(116-225*eps2)/384; 00971 d *= eps; 00972 c[4] = d*(2695-7173*eps2)/7680; 00973 d *= eps; 00974 c[5] = 3467*d/7680; 00975 d *= eps; 00976 c[6] = 38081*d/61440; 00977 break; 00978 case 7: 00979 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728; 00980 d *= eps; 00981 c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288; 00982 d *= eps; 00983 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288; 00984 d *= eps; 00985 c[4] = d*(2695-7173*eps2)/7680; 00986 d *= eps; 00987 c[5] = d*(41604-141115*eps2)/92160; 00988 d *= eps; 00989 c[6] = 38081*d/61440; 00990 d *= eps; 00991 c[7] = 459485*d/516096; 00992 break; 00993 case 8: 00994 c[1] = d*(eps2*((9840-4879*eps2)*eps2-20736)+36864)/73728; 00995 d *= eps; 00996 c[2] = d*(eps2*((120150-86171*eps2)*eps2-142080)+115200)/368640; 00997 d *= eps; 00998 c[3] = d*(eps2*(8703*eps2-7200)+3712)/12288; 00999 d *= eps; 01000 c[4] = d*(eps2*(1082857*eps2-688608)+258720)/737280; 01001 d *= eps; 01002 c[5] = d*(41604-141115*eps2)/92160; 01003 d *= eps; 01004 c[6] = d*(533134-2200311*eps2)/860160; 01005 d *= eps; 01006 c[7] = 459485*d/516096; 01007 d *= eps; 01008 c[8] = 109167851*d/82575360; 01009 break; 01010 default: 01011 STATIC_ASSERT(nC1p_ >= 0 && nC1p_ <= 8, "Bad value of nC1p_"); 01012 } 01013 } 01014 01015 // The scale factor A2-1 = mean value of I2-1 01016 Math::real Geodesic::A2m1f(real eps) throw() { 01017 real 01018 eps2 = Math::sq(eps), 01019 t; 01020 switch (nA2_/2) { 01021 case 0: 01022 t = 0; 01023 break; 01024 case 1: 01025 t = eps2/4; 01026 break; 01027 case 2: 01028 t = eps2*(9*eps2+16)/64; 01029 break; 01030 case 3: 01031 t = eps2*(eps2*(25*eps2+36)+64)/256; 01032 break; 01033 case 4: 01034 t = eps2*(eps2*(eps2*(1225*eps2+1600)+2304)+4096)/16384; 01035 break; 01036 default: 01037 STATIC_ASSERT(nA2_ >= 0 && nA2_ <= 8, "Bad value of nA2_"); 01038 t = 0; 01039 } 01040 return t * (1 - eps) - eps; 01041 } 01042 01043 // The coefficients C2[l] in the Fourier expansion of B2 01044 void Geodesic::C2f(real eps, real c[]) throw() { 01045 real 01046 eps2 = Math::sq(eps), 01047 d = eps; 01048 switch (nC2_) { 01049 case 0: 01050 break; 01051 case 1: 01052 c[1] = d/2; 01053 break; 01054 case 2: 01055 c[1] = d/2; 01056 d *= eps; 01057 c[2] = 3*d/16; 01058 break; 01059 case 3: 01060 c[1] = d*(eps2+8)/16; 01061 d *= eps; 01062 c[2] = 3*d/16; 01063 d *= eps; 01064 c[3] = 5*d/48; 01065 break; 01066 case 4: 01067 c[1] = d*(eps2+8)/16; 01068 d *= eps; 01069 c[2] = d*(eps2+6)/32; 01070 d *= eps; 01071 c[3] = 5*d/48; 01072 d *= eps; 01073 c[4] = 35*d/512; 01074 break; 01075 case 5: 01076 c[1] = d*(eps2*(eps2+2)+16)/32; 01077 d *= eps; 01078 c[2] = d*(eps2+6)/32; 01079 d *= eps; 01080 c[3] = d*(15*eps2+80)/768; 01081 d *= eps; 01082 c[4] = 35*d/512; 01083 d *= eps; 01084 c[5] = 63*d/1280; 01085 break; 01086 case 6: 01087 c[1] = d*(eps2*(eps2+2)+16)/32; 01088 d *= eps; 01089 c[2] = d*(eps2*(35*eps2+64)+384)/2048; 01090 d *= eps; 01091 c[3] = d*(15*eps2+80)/768; 01092 d *= eps; 01093 c[4] = d*(7*eps2+35)/512; 01094 d *= eps; 01095 c[5] = 63*d/1280; 01096 d *= eps; 01097 c[6] = 77*d/2048; 01098 break; 01099 case 7: 01100 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048; 01101 d *= eps; 01102 c[2] = d*(eps2*(35*eps2+64)+384)/2048; 01103 d *= eps; 01104 c[3] = d*(eps2*(69*eps2+120)+640)/6144; 01105 d *= eps; 01106 c[4] = d*(7*eps2+35)/512; 01107 d *= eps; 01108 c[5] = d*(105*eps2+504)/10240; 01109 d *= eps; 01110 c[6] = 77*d/2048; 01111 d *= eps; 01112 c[7] = 429*d/14336; 01113 break; 01114 case 8: 01115 c[1] = d*(eps2*(eps2*(41*eps2+64)+128)+1024)/2048; 01116 d *= eps; 01117 c[2] = d*(eps2*(eps2*(47*eps2+70)+128)+768)/4096; 01118 d *= eps; 01119 c[3] = d*(eps2*(69*eps2+120)+640)/6144; 01120 d *= eps; 01121 c[4] = d*(eps2*(133*eps2+224)+1120)/16384; 01122 d *= eps; 01123 c[5] = d*(105*eps2+504)/10240; 01124 d *= eps; 01125 c[6] = d*(33*eps2+154)/4096; 01126 d *= eps; 01127 c[7] = 429*d/14336; 01128 d *= eps; 01129 c[8] = 6435*d/262144; 01130 break; 01131 default: 01132 STATIC_ASSERT(nC2_ >= 0 && nC2_ <= 8, "Bad value of nC2_"); 01133 } 01134 } 01135 01136 // The scale factor A3 = mean value of I3 01137 void Geodesic::A3coeff() throw() { 01138 switch (nA3_) { 01139 case 0: 01140 break; 01141 case 1: 01142 _A3x[0] = 1; 01143 break; 01144 case 2: 01145 _A3x[0] = 1; 01146 _A3x[1] = -1/real(2); 01147 break; 01148 case 3: 01149 _A3x[0] = 1; 01150 _A3x[1] = (_n-1)/2; 01151 _A3x[2] = -1/real(4); 01152 break; 01153 case 4: 01154 _A3x[0] = 1; 01155 _A3x[1] = (_n-1)/2; 01156 _A3x[2] = (-_n-2)/8; 01157 _A3x[3] = -1/real(16); 01158 break; 01159 case 5: 01160 _A3x[0] = 1; 01161 _A3x[1] = (_n-1)/2; 01162 _A3x[2] = (_n*(3*_n-1)-2)/8; 01163 _A3x[3] = (-3*_n-1)/16; 01164 _A3x[4] = -3/real(64); 01165 break; 01166 case 6: 01167 _A3x[0] = 1; 01168 _A3x[1] = (_n-1)/2; 01169 _A3x[2] = (_n*(3*_n-1)-2)/8; 01170 _A3x[3] = ((-_n-3)*_n-1)/16; 01171 _A3x[4] = (-2*_n-3)/64; 01172 _A3x[5] = -3/real(128); 01173 break; 01174 case 7: 01175 _A3x[0] = 1; 01176 _A3x[1] = (_n-1)/2; 01177 _A3x[2] = (_n*(3*_n-1)-2)/8; 01178 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16; 01179 _A3x[4] = ((-10*_n-2)*_n-3)/64; 01180 _A3x[5] = (-5*_n-3)/128; 01181 _A3x[6] = -5/real(256); 01182 break; 01183 case 8: 01184 _A3x[0] = 1; 01185 _A3x[1] = (_n-1)/2; 01186 _A3x[2] = (_n*(3*_n-1)-2)/8; 01187 _A3x[3] = (_n*(_n*(5*_n-1)-3)-1)/16; 01188 _A3x[4] = (_n*((-5*_n-20)*_n-4)-6)/128; 01189 _A3x[5] = ((-5*_n-10)*_n-6)/256; 01190 _A3x[6] = (-15*_n-20)/1024; 01191 _A3x[7] = -25/real(2048); 01192 break; 01193 default: 01194 STATIC_ASSERT(nA3_ >= 0 && nA3_ <= 8, "Bad value of nA3_"); 01195 } 01196 } 01197 01198 // The coefficients C3[l] in the Fourier expansion of B3 01199 void Geodesic::C3coeff() throw() { 01200 switch (nC3_) { 01201 case 0: 01202 break; 01203 case 1: 01204 break; 01205 case 2: 01206 _C3x[0] = 1/real(4); 01207 break; 01208 case 3: 01209 _C3x[0] = (1-_n)/4; 01210 _C3x[1] = 1/real(8); 01211 _C3x[2] = 1/real(16); 01212 break; 01213 case 4: 01214 _C3x[0] = (1-_n)/4; 01215 _C3x[1] = 1/real(8); 01216 _C3x[2] = 3/real(64); 01217 _C3x[3] = (2-3*_n)/32; 01218 _C3x[4] = 3/real(64); 01219 _C3x[5] = 5/real(192); 01220 break; 01221 case 5: 01222 _C3x[0] = (1-_n)/4; 01223 _C3x[1] = (1-_n*_n)/8; 01224 _C3x[2] = (3*_n+3)/64; 01225 _C3x[3] = 5/real(128); 01226 _C3x[4] = ((_n-3)*_n+2)/32; 01227 _C3x[5] = (3-2*_n)/64; 01228 _C3x[6] = 3/real(128); 01229 _C3x[7] = (5-9*_n)/192; 01230 _C3x[8] = 3/real(128); 01231 _C3x[9] = 7/real(512); 01232 break; 01233 case 6: 01234 _C3x[0] = (1-_n)/4; 01235 _C3x[1] = (1-_n*_n)/8; 01236 _C3x[2] = ((3-_n)*_n+3)/64; 01237 _C3x[3] = (2*_n+5)/128; 01238 _C3x[4] = 3/real(128); 01239 _C3x[5] = ((_n-3)*_n+2)/32; 01240 _C3x[6] = ((-3*_n-2)*_n+3)/64; 01241 _C3x[7] = (_n+3)/128; 01242 _C3x[8] = 5/real(256); 01243 _C3x[9] = (_n*(5*_n-9)+5)/192; 01244 _C3x[10] = (9-10*_n)/384; 01245 _C3x[11] = 7/real(512); 01246 _C3x[12] = (7-14*_n)/512; 01247 _C3x[13] = 7/real(512); 01248 _C3x[14] = 21/real(2560); 01249 break; 01250 case 7: 01251 _C3x[0] = (1-_n)/4; 01252 _C3x[1] = (1-_n*_n)/8; 01253 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64; 01254 _C3x[3] = (_n*(2*_n+2)+5)/128; 01255 _C3x[4] = (11*_n+12)/512; 01256 _C3x[5] = 21/real(1024); 01257 _C3x[6] = ((_n-3)*_n+2)/32; 01258 _C3x[7] = (_n*(_n*(2*_n-3)-2)+3)/64; 01259 _C3x[8] = ((2-9*_n)*_n+6)/256; 01260 _C3x[9] = (_n+5)/256; 01261 _C3x[10] = 27/real(2048); 01262 _C3x[11] = (_n*((5-_n)*_n-9)+5)/192; 01263 _C3x[12] = ((-6*_n-10)*_n+9)/384; 01264 _C3x[13] = (21-4*_n)/1536; 01265 _C3x[14] = 3/real(256); 01266 _C3x[15] = (_n*(10*_n-14)+7)/512; 01267 _C3x[16] = (7-10*_n)/512; 01268 _C3x[17] = 9/real(1024); 01269 _C3x[18] = (21-45*_n)/2560; 01270 _C3x[19] = 9/real(1024); 01271 _C3x[20] = 11/real(2048); 01272 break; 01273 case 8: 01274 _C3x[0] = (1-_n)/4; 01275 _C3x[1] = (1-_n*_n)/8; 01276 _C3x[2] = (_n*((-5*_n-1)*_n+3)+3)/64; 01277 _C3x[3] = (_n*((2-2*_n)*_n+2)+5)/128; 01278 _C3x[4] = (_n*(3*_n+11)+12)/512; 01279 _C3x[5] = (10*_n+21)/1024; 01280 _C3x[6] = 243/real(16384); 01281 _C3x[7] = ((_n-3)*_n+2)/32; 01282 _C3x[8] = (_n*(_n*(2*_n-3)-2)+3)/64; 01283 _C3x[9] = (_n*((-6*_n-9)*_n+2)+6)/256; 01284 _C3x[10] = ((1-2*_n)*_n+5)/256; 01285 _C3x[11] = (69*_n+108)/8192; 01286 _C3x[12] = 187/real(16384); 01287 _C3x[13] = (_n*((5-_n)*_n-9)+5)/192; 01288 _C3x[14] = (_n*(_n*(10*_n-6)-10)+9)/384; 01289 _C3x[15] = ((-77*_n-8)*_n+42)/3072; 01290 _C3x[16] = (12-_n)/1024; 01291 _C3x[17] = 139/real(16384); 01292 _C3x[18] = (_n*((20-7*_n)*_n-28)+14)/1024; 01293 _C3x[19] = ((-7*_n-40)*_n+28)/2048; 01294 _C3x[20] = (72-43*_n)/8192; 01295 _C3x[21] = 127/real(16384); 01296 _C3x[22] = (_n*(75*_n-90)+42)/5120; 01297 _C3x[23] = (9-15*_n)/1024; 01298 _C3x[24] = 99/real(16384); 01299 _C3x[25] = (44-99*_n)/8192; 01300 _C3x[26] = 99/real(16384); 01301 _C3x[27] = 429/real(114688); 01302 break; 01303 default: 01304 STATIC_ASSERT(nC3_ >= 0 && nC3_ <= 8, "Bad value of nC3_"); 01305 } 01306 } 01307 01308 // The coefficients C4[l] in the Fourier expansion of I4 01309 void Geodesic::C4coeff() throw() { 01310 switch (nC4_) { 01311 case 0: 01312 break; 01313 case 1: 01314 _C4x[0] = 2/real(3); 01315 break; 01316 case 2: 01317 _C4x[0] = (10-_ep2)/15; 01318 _C4x[1] = -1/real(20); 01319 _C4x[2] = 1/real(180); 01320 break; 01321 case 3: 01322 _C4x[0] = (_ep2*(4*_ep2-7)+70)/105; 01323 _C4x[1] = (4*_ep2-7)/140; 01324 _C4x[2] = 1/real(42); 01325 _C4x[3] = (7-4*_ep2)/1260; 01326 _C4x[4] = -1/real(252); 01327 _C4x[5] = 1/real(2100); 01328 break; 01329 case 4: 01330 _C4x[0] = (_ep2*((12-8*_ep2)*_ep2-21)+210)/315; 01331 _C4x[1] = ((12-8*_ep2)*_ep2-21)/420; 01332 _C4x[2] = (3-2*_ep2)/126; 01333 _C4x[3] = -1/real(72); 01334 _C4x[4] = (_ep2*(8*_ep2-12)+21)/3780; 01335 _C4x[5] = (2*_ep2-3)/756; 01336 _C4x[6] = 1/real(360); 01337 _C4x[7] = (3-2*_ep2)/6300; 01338 _C4x[8] = -1/real(1800); 01339 _C4x[9] = 1/real(17640); 01340 break; 01341 case 5: 01342 _C4x[0] = (_ep2*(_ep2*(_ep2*(64*_ep2-88)+132)-231)+2310)/3465; 01343 _C4x[1] = (_ep2*(_ep2*(64*_ep2-88)+132)-231)/4620; 01344 _C4x[2] = (_ep2*(16*_ep2-22)+33)/1386; 01345 _C4x[3] = (8*_ep2-11)/792; 01346 _C4x[4] = 1/real(110); 01347 _C4x[5] = (_ep2*((88-64*_ep2)*_ep2-132)+231)/41580; 01348 _C4x[6] = ((22-16*_ep2)*_ep2-33)/8316; 01349 _C4x[7] = (11-8*_ep2)/3960; 01350 _C4x[8] = -1/real(495); 01351 _C4x[9] = (_ep2*(16*_ep2-22)+33)/69300; 01352 _C4x[10] = (8*_ep2-11)/19800; 01353 _C4x[11] = 1/real(1925); 01354 _C4x[12] = (11-8*_ep2)/194040; 01355 _C4x[13] = -1/real(10780); 01356 _C4x[14] = 1/real(124740); 01357 break; 01358 case 6: 01359 _C4x[0] = (_ep2*(_ep2*(_ep2*((832-640*_ep2)*_ep2-1144)+1716)-3003)+ 01360 30030)/45045; 01361 _C4x[1] = (_ep2*(_ep2*((832-640*_ep2)*_ep2-1144)+1716)-3003)/60060; 01362 _C4x[2] = (_ep2*((208-160*_ep2)*_ep2-286)+429)/18018; 01363 _C4x[3] = ((104-80*_ep2)*_ep2-143)/10296; 01364 _C4x[4] = (13-10*_ep2)/1430; 01365 _C4x[5] = -1/real(156); 01366 _C4x[6] = (_ep2*(_ep2*(_ep2*(640*_ep2-832)+1144)-1716)+3003)/540540; 01367 _C4x[7] = (_ep2*(_ep2*(160*_ep2-208)+286)-429)/108108; 01368 _C4x[8] = (_ep2*(80*_ep2-104)+143)/51480; 01369 _C4x[9] = (10*_ep2-13)/6435; 01370 _C4x[10] = 5/real(3276); 01371 _C4x[11] = (_ep2*((208-160*_ep2)*_ep2-286)+429)/900900; 01372 _C4x[12] = ((104-80*_ep2)*_ep2-143)/257400; 01373 _C4x[13] = (13-10*_ep2)/25025; 01374 _C4x[14] = -1/real(2184); 01375 _C4x[15] = (_ep2*(80*_ep2-104)+143)/2522520; 01376 _C4x[16] = (10*_ep2-13)/140140; 01377 _C4x[17] = 5/real(45864); 01378 _C4x[18] = (13-10*_ep2)/1621620; 01379 _C4x[19] = -1/real(58968); 01380 _C4x[20] = 1/real(792792); 01381 break; 01382 case 7: 01383 _C4x[0] = (_ep2*(_ep2*(_ep2*(_ep2*(_ep2*(512*_ep2-640)+832)-1144)+1716)- 01384 3003)+30030)/45045; 01385 _C4x[1] = (_ep2*(_ep2*(_ep2*(_ep2*(512*_ep2-640)+832)-1144)+1716)- 01386 3003)/60060; 01387 _C4x[2] = (_ep2*(_ep2*(_ep2*(128*_ep2-160)+208)-286)+429)/18018; 01388 _C4x[3] = (_ep2*(_ep2*(64*_ep2-80)+104)-143)/10296; 01389 _C4x[4] = (_ep2*(8*_ep2-10)+13)/1430; 01390 _C4x[5] = (4*_ep2-5)/780; 01391 _C4x[6] = 1/real(210); 01392 _C4x[7] = (_ep2*(_ep2*(_ep2*((640-512*_ep2)*_ep2-832)+1144)-1716)+ 01393 3003)/540540; 01394 _C4x[8] = (_ep2*(_ep2*((160-128*_ep2)*_ep2-208)+286)-429)/108108; 01395 _C4x[9] = (_ep2*((80-64*_ep2)*_ep2-104)+143)/51480; 01396 _C4x[10] = ((10-8*_ep2)*_ep2-13)/6435; 01397 _C4x[11] = (5-4*_ep2)/3276; 01398 _C4x[12] = -1/real(840); 01399 _C4x[13] = (_ep2*(_ep2*(_ep2*(128*_ep2-160)+208)-286)+429)/900900; 01400 _C4x[14] = (_ep2*(_ep2*(64*_ep2-80)+104)-143)/257400; 01401 _C4x[15] = (_ep2*(8*_ep2-10)+13)/25025; 01402 _C4x[16] = (4*_ep2-5)/10920; 01403 _C4x[17] = 1/real(2520); 01404 _C4x[18] = (_ep2*((80-64*_ep2)*_ep2-104)+143)/2522520; 01405 _C4x[19] = ((10-8*_ep2)*_ep2-13)/140140; 01406 _C4x[20] = (5-4*_ep2)/45864; 01407 _C4x[21] = -1/real(8820); 01408 _C4x[22] = (_ep2*(8*_ep2-10)+13)/1621620; 01409 _C4x[23] = (4*_ep2-5)/294840; 01410 _C4x[24] = 1/real(41580); 01411 _C4x[25] = (5-4*_ep2)/3963960; 01412 _C4x[26] = -1/real(304920); 01413 _C4x[27] = 1/real(4684680); 01414 break; 01415 case 8: 01416 _C4x[0] = (_ep2*(_ep2*(_ep2*(_ep2*(_ep2*((8704-7168*_ep2)*_ep2-10880)+ 01417 14144)-19448)+29172)-51051)+510510)/765765; 01418 _C4x[1] = (_ep2*(_ep2*(_ep2*(_ep2*((8704-7168*_ep2)*_ep2-10880)+14144)- 01419 19448)+29172)-51051)/1021020; 01420 _C4x[2] = (_ep2*(_ep2*(_ep2*((2176-1792*_ep2)*_ep2-2720)+3536)-4862)+ 01421 7293)/306306; 01422 _C4x[3] = (_ep2*(_ep2*((1088-896*_ep2)*_ep2-1360)+1768)-2431)/175032; 01423 _C4x[4] = (_ep2*((136-112*_ep2)*_ep2-170)+221)/24310; 01424 _C4x[5] = ((68-56*_ep2)*_ep2-85)/13260; 01425 _C4x[6] = (17-14*_ep2)/3570; 01426 _C4x[7] = -1/real(272); 01427 _C4x[8] = (_ep2*(_ep2*(_ep2*(_ep2*(_ep2*(7168*_ep2-8704)+10880)-14144)+ 01428 19448)-29172)+51051)/9189180; 01429 _C4x[9] = (_ep2*(_ep2*(_ep2*(_ep2*(1792*_ep2-2176)+2720)-3536)+4862)- 01430 7293)/1837836; 01431 _C4x[10] = (_ep2*(_ep2*(_ep2*(896*_ep2-1088)+1360)-1768)+2431)/875160; 01432 _C4x[11] = (_ep2*(_ep2*(112*_ep2-136)+170)-221)/109395; 01433 _C4x[12] = (_ep2*(56*_ep2-68)+85)/55692; 01434 _C4x[13] = (14*_ep2-17)/14280; 01435 _C4x[14] = 7/real(7344); 01436 _C4x[15] = (_ep2*(_ep2*(_ep2*((2176-1792*_ep2)*_ep2-2720)+3536)-4862)+ 01437 7293)/15315300; 01438 _C4x[16] = (_ep2*(_ep2*((1088-896*_ep2)*_ep2-1360)+1768)-2431)/4375800; 01439 _C4x[17] = (_ep2*((136-112*_ep2)*_ep2-170)+221)/425425; 01440 _C4x[18] = ((68-56*_ep2)*_ep2-85)/185640; 01441 _C4x[19] = (17-14*_ep2)/42840; 01442 _C4x[20] = -7/real(20400); 01443 _C4x[21] = (_ep2*(_ep2*(_ep2*(896*_ep2-1088)+1360)-1768)+2431)/42882840; 01444 _C4x[22] = (_ep2*(_ep2*(112*_ep2-136)+170)-221)/2382380; 01445 _C4x[23] = (_ep2*(56*_ep2-68)+85)/779688; 01446 _C4x[24] = (14*_ep2-17)/149940; 01447 _C4x[25] = 1/real(8976); 01448 _C4x[26] = (_ep2*((136-112*_ep2)*_ep2-170)+221)/27567540; 01449 _C4x[27] = ((68-56*_ep2)*_ep2-85)/5012280; 01450 _C4x[28] = (17-14*_ep2)/706860; 01451 _C4x[29] = -7/real(242352); 01452 _C4x[30] = (_ep2*(56*_ep2-68)+85)/67387320; 01453 _C4x[31] = (14*_ep2-17)/5183640; 01454 _C4x[32] = 7/real(1283568); 01455 _C4x[33] = (17-14*_ep2)/79639560; 01456 _C4x[34] = -1/real(1516944); 01457 _C4x[35] = 1/real(26254800); 01458 break; 01459 default: 01460 STATIC_ASSERT(nC3_ >= 0 && nC4_ <= 8, "Bad value of nC4_"); 01461 } 01462 } 01463 01464 } // namespace GeographicLib