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