GeographicLib  1.21
Geodesic.cpp
Go to the documentation of this file.
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