GeographicLib  1.21
SphericalEngine.cpp
Go to the documentation of this file.
00001 /**
00002  * \file SphericalEngine.cpp
00003  * \brief Implementation for GeographicLib::SphericalEngine class
00004  *
00005  * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under
00006  * the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  *
00009  * The general sum is\verbatim
00010  V(r, theta, lambda) = sum(n = 0..N) sum(m = 0..n)
00011    q^(n+1) * (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) * P[n,m](t)
00012 \endverbatim
00013  * where <tt>t = cos(theta)</tt>, <tt>q = a/r</tt>.  In addition write <tt>u =
00014  * sin(theta)</tt>.
00015  *
00016  * <tt>P[n,m]</tt> is a normalized associated Legendre function of degree
00017  * <tt>n</tt> and order <tt>m</tt>.  Here the formulas are given for full
00018  * normalized functions (usually denoted <tt>Pbar</tt>).
00019  *
00020  * Rewrite outer sum\verbatim
00021  V(r, theta, lambda) = sum(m = 0..N) * P[m,m](t) * q^(m+1) *
00022     [Sc[m] * cos(m*lambda) + Ss[m] * sin(m*lambda)]
00023 \endverbatim
00024  * where the inner sums are\verbatim
00025    Sc[m] = sum(n = m..N) q^(n-m) * C[n,m] * P[n,m](t)/P[m,m](t)
00026    Ss[m] = sum(n = m..N) q^(n-m) * S[n,m] * P[n,m](t)/P[m,m](t)
00027 \endverbatim
00028  * Evaluate sums via Clenshaw method.  The overall framework is similar to
00029  * Deakin with the following changes:
00030  * - Clenshaw summation is used to roll the computation of
00031  *   <tt>cos(m*lambda)</tt> and <tt>sin(m*lambda)</tt> into the evaluation of
00032  *   the outer sum (rather than independently computing an array of these
00033  *   trigonometric terms).
00034  * - Scale the coefficients to guard against overflow when <tt>N</tt> is large.
00035  * .
00036  * For the general framework of Clenshaw, see
00037  * http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
00038  *
00039  * Let\verbatim
00040     S = sum(k = 0..N) c[k] * F[k](x)
00041     F[n+1](x) = alpha[n](x) * F[n](x) + beta[n](x) * F[n-1](x)
00042 \endverbatim
00043  * Evaluate <tt>S</tt> with\verbatim
00044     y[N+2] = y[N+1] = 0
00045     y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k]
00046     S = c[0] * F[0] + y[1] * F[1] + beta[1] * F[0] * y[2]
00047 \endverbatim
00048  * \e IF <tt>F[0](x) = 1</tt> and <tt>beta(0,x) = 0</tt>, then <tt>F[1](x) =
00049  * alpha(0,x)</tt> and we can continue the recursion for <tt>y[k]</tt> until
00050  * <tt>y[0]</tt>, giving\verbatim
00051     S = y[0]
00052 \endverbatim
00053  *
00054  * Evaluating the inner sum\verbatim
00055  l = n-m; n = l+m
00056  Sc[m] = sum(l = 0..N-m) C[l+m,m] * q^l * P[l+m,m](t)/P[m,m](t)
00057  F[l] = q^l * P[l+m,m](t)/P[m,m](t)
00058 \endverbatim
00059  * Holmes + Featherstone, Eq. (11), give\verbatim
00060    P[n,m] = sqrt((2*n-1)*(2*n+1)/((n-m)*(n+m))) * t * P[n-1,m] -
00061             sqrt((2*n+1)*(n+m-1)*(n-m-1)/((n-m)*(n+m)*(2*n-3))) * P[n-2,m]
00062 \endverbatim
00063  * thus\verbatim
00064    alpha[l] = t * q * sqrt(((2*n+1)*(2*n+3))/
00065                            ((n-m+1)*(n+m+1)))
00066    beta[l+1] = - q^2 * sqrt(((n-m+1)*(n+m+1)*(2*n+5))/
00067                             ((n-m+2)*(n+m+2)*(2*n+1)))
00068 \endverbatim
00069  * In this case, <tt>F[0] = 1</tt> and <tt>beta[0] = 0</tt>, so the <tt>Sc[m]
00070  * = y[0]</tt>.
00071  *
00072  * Evaluating the outer sum\verbatim
00073  V = sum(m = 0..N) Sc[m] * q^(m+1) * cos(m*lambda) * P[m,m](t)
00074    + sum(m = 0..N) Ss[m] * q^(m+1) * cos(m*lambda) * P[m,m](t)
00075  F[m] = q^(m+1) * cos(m*lambda) * P[m,m](t) [or sin(m*lambda)]
00076 \endverbatim
00077  * Holmes + Featherstone, Eq. (13), give\verbatim
00078    P[m,m] = u * sqrt((2*m+1)/((m>1?2:1)*m)) * P[m-1,m-1]
00079 \endverbatim
00080  * also, we have\verbatim
00081    cos((m+1)*lambda) = 2*cos(lambda)*cos(m*lambda) - cos((m-1)*lambda)
00082 \endverbatim
00083  * thus\verbatim
00084    alpha[m] = 2*cos(lambda) * sqrt((2*m+3)/(2*(m+1))) * u * q
00085             =   cos(lambda) * sqrt( 2*(2*m+3)/(m+1) ) * u * q
00086    beta[m+1] = -sqrt((2*m+3)*(2*m+5)/(4*(m+1)*(m+2))) * u^2 * q^2
00087                * (m == 0 ? sqrt(2) : 1)
00088 \endverbatim
00089  * Thus\verbatim
00090  F[0] = q                                [or 0]
00091  F[1] = cos(lambda) * sqrt(3) * u * q^2  [or sin(lambda)]
00092  beta[1] = - sqrt(15/4) * u^2 * q^2
00093 \endverbatim
00094  *
00095  * Here is how the various components of the gradient are computed
00096  *
00097  * Differentiate wrt <tt>r</tt>\verbatim
00098    d q^(n+1) / dr = (-1/r) * (n+1) * q^(n+1)
00099 \endverbatim
00100  * so multiply <tt>C[n,m]</tt> by <tt>n+1</tt> in inner sum and multiply the
00101  * sum by <tt>-1/r</tt>.
00102  *
00103  * Differentiate wrt <tt>lambda</tt>\verbatim
00104    d cos(m*lambda) = -m * sin(m*lambda)
00105    d sin(m*lambda) =  m * cos(m*lambda)
00106 \endverbatim
00107  * so multiply terms by <tt>m</tt> in outer sum and swap sine and cosine
00108  * variables.
00109  *
00110  * Differentiate wrt <tt>theta</tt>\verbatim
00111   dV/dtheta = V' = -u * dV/dt = -u * V'
00112 \endverbatim
00113  * here <tt>'</tt> denotes differentiation wrt <tt>theta</tt>.\verbatim
00114    d/dtheta (Sc[m] * P[m,m](t)) = Sc'[m] * P[m,m](t) + Sc[m] * P'[m,m](t)
00115 \endverbatim
00116  * Now <tt>P[m,m](t) = const * u^m</tt>, so <tt>P'[m,m](t) = m * t/u *
00117  * P[m,m](t)</tt>, thus\verbatim
00118    d/dtheta (Sc[m] * P[m,m](t)) = (Sc'[m] + m * t/u * Sc[m]) * P[m,m](t)
00119 \endverbatim
00120  * Clenshaw recursion for <tt>Sc[m]</tt> reads\verbatim
00121     y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k]
00122 \endverbatim
00123  * Substituting <tt>alpha[k] = const * t</tt>, <tt>alpha'[k] = -u/t *
00124  * alpha[k]</tt>, <tt>beta'[k] = c'[k] = 0</tt> gives\verbatim
00125     y'[k] = alpha[k] * y'[k+1] + beta[k+1] * y'[k+2] - u/t * alpha[k] * y[k+1]
00126 \endverbatim
00127  *
00128  * Finally, given the derivatives of <tt>V</tt>, we can compute the components
00129  * of the gradient in spherical coordinates and transform the result into
00130  * cartesian coordinates.
00131  **********************************************************************/
00132 
00133 #include <GeographicLib/SphericalEngine.hpp>
00134 #include <limits>
00135 #include <algorithm>
00136 #include <GeographicLib/CircularEngine.hpp>
00137 #include <GeographicLib/Utility.hpp>
00138 
00139 #define GEOGRAPHICLIB_SPHERICALENGINE_CPP \
00140   "$Id: 861a718d09c83cdd9bf58939d938a25797b9f306 $"
00141 
00142 RCSID_DECL(GEOGRAPHICLIB_SPHERICALENGINE_CPP)
00143 RCSID_DECL(GEOGRAPHICLIB_SPHERICALENGINE_HPP)
00144 
00145 namespace GeographicLib {
00146 
00147   using namespace std;
00148 
00149   const Math::real SphericalEngine::scale_ =
00150     pow(real(numeric_limits<real>::radix),
00151         -3 * numeric_limits<real>::max_exponent / 5);
00152   const Math::real SphericalEngine::eps_ =
00153     numeric_limits<real>::epsilon() * sqrt(numeric_limits<real>::epsilon());
00154 
00155   const vector<Math::real> SphericalEngine::Z_(0);
00156   vector<Math::real> SphericalEngine::root_(0);
00157 
00158   template<bool gradp, SphericalEngine::normalization norm, int L>
00159   Math::real SphericalEngine::Value(const coeff c[], const real f[],
00160                                     real x, real y, real z, real a,
00161                                     real& gradx, real& grady, real& gradz)
00162     throw() {
00163     STATIC_ASSERT(L > 0, "L must be positive");
00164     STATIC_ASSERT(norm == FULL || norm == SCHMIDT, "Unknown normalization");
00165     int N = c[0].nmx(), M = c[0].mmx();
00166 
00167     real
00168       p = Math::hypot(x, y),
00169       cl = p ? x / p : 1,       // cos(lambda); at pole, pick lambda = 0
00170       sl = p ? y / p : 0,       // sin(lambda)
00171       r = Math::hypot(z, p),
00172       t = r ? z / r : 0,            // cos(theta); at origin, pick theta = pi/2
00173       u = r ? max(p / r, eps_) : 1, // sin(theta); but avoid the pole
00174       q = a / r;
00175     real
00176       q2 = Math::sq(q),
00177       uq = u * q,
00178       uq2 = Math::sq(uq),
00179       tu = t / u;
00180     // Initialize outer sum
00181     real vc  = 0, vc2  = 0, vs  = 0, vs2  = 0;   // v [N + 1], v [N + 2]
00182     // vr, vt, vl and similar w variable accumulate the sums for the
00183     // derivatives wrt r, theta, and lambda, respectively.
00184     real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0;   // vr[N + 1], vr[N + 2]
00185     real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0;   // vt[N + 1], vt[N + 2]
00186     real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0;   // vl[N + 1], vl[N + 2]
00187     int k[L];
00188     for (int m = M; m >= 0; --m) {   // m = M .. 0
00189       // Initialize inner sum
00190       real wc  = 0, wc2  = 0, ws  = 0, ws2  = 0; // w [N - m + 1], w [N - m + 2]
00191       real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0; // wr[N - m + 1], wr[N - m + 2]
00192       real wtc = 0, wtc2 = 0, wts = 0, wts2 = 0; // wt[N - m + 1], wt[N - m + 2]
00193       for (int l = 0; l < L; ++l)
00194         k[l] = c[l].index(N, m) + 1;
00195       for (int n = N; n >= m; --n) {             // n = N .. m; l = N - m .. 0
00196         real w, A, Ax, B, R;    // alpha[l], beta[l + 1]
00197         switch (norm) {
00198         case FULL:
00199           w = root_[2 * n + 1] / (root_[n - m + 1] * root_[n + m + 1]);
00200           Ax = q * w * root_[2 * n + 3];
00201           A = t * Ax;
00202           B = - q2 * root_[2 * n + 5] /
00203             (w * root_[n - m + 2] * root_[n + m + 2]);
00204           break;
00205         case SCHMIDT:
00206           w = root_[n - m + 1] * root_[n + m + 1];
00207           Ax = q * (2 * n + 1) / w;
00208           A = t * Ax;
00209           B = - q2 * w / (root_[n - m + 2] * root_[n + m + 2]);
00210           break;
00211         }
00212         R = c[0].Cv(--k[0]);
00213         for (int l = 1; l < L; ++l)
00214           R += c[l].Cv(--k[l], n, m, f[l]);
00215         R *= scale_;
00216         w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
00217         if (gradp) {
00218           w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
00219           w = A * wtc + B * wtc2 -  u*Ax * wc2; wtc2 = wtc; wtc = w;
00220         }
00221         if (m) {
00222           R = c[0].Sv(k[0]);
00223           for (int l = 1; l < L; ++l)
00224             R += c[l].Sv(k[l], n, m, f[l]);
00225           R *= scale_;
00226           w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
00227           if (gradp) {
00228             w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
00229             w = A * wts + B * wts2 -  u*Ax * ws2; wts2 = wts; wts = w;
00230           }
00231         }
00232       }
00233       // Now Sc[m] = wc, Ss[m] = ws
00234       // Sc'[m] = wtc, Ss'[m] = wtc
00235       if (m) {
00236         real v, A, B;           // alpha[m], beta[m + 1]
00237         switch (norm) {
00238         case FULL:
00239           v = root_[2] * root_[2 * m + 3] / root_[m + 1];
00240           A = cl * v * uq;
00241           B = - v * root_[2 * m + 5] / (root_[8] * root_[m + 2]) * uq2;
00242           break;
00243         case SCHMIDT:
00244           v = root_[2] * root_[2 * m + 1] / root_[m + 1];
00245           A = cl * v * uq;
00246           B = - v * root_[2 * m + 3] / (root_[8] * root_[m + 2]) * uq2;
00247           break;
00248         }
00249         v = A * vc  + B * vc2  +  wc ; vc2  = vc ; vc  = v;
00250         v = A * vs  + B * vs2  +  ws ; vs2  = vs ; vs  = v;
00251         if (gradp) {
00252           // Include the terms Sc[m] * P'[m,m](t) and Ss[m] * P'[m,m](t)
00253           wtc += m * tu * wc; wts += m * tu * ws;
00254           v = A * vrc + B * vrc2 +  wrc; vrc2 = vrc; vrc = v;
00255           v = A * vrs + B * vrs2 +  wrs; vrs2 = vrs; vrs = v;
00256           v = A * vtc + B * vtc2 +  wtc; vtc2 = vtc; vtc = v;
00257           v = A * vts + B * vts2 +  wts; vts2 = vts; vts = v;
00258           v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;
00259           v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;
00260         }
00261       } else {
00262         real A, B, qs;
00263         switch (norm) {
00264         case FULL:
00265           A = root_[3] * uq;       // F[1]/(q*cl) or F[1]/(q*sl)
00266           B = - root_[15]/2 * uq2; // beta[1]/q
00267           break;
00268         case SCHMIDT:
00269           A = uq;
00270           B = - root_[3]/2 * uq2;
00271           break;
00272         }
00273         qs = q / scale_;
00274         vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);
00275         if (gradp) {
00276           qs /= r;
00277           // The components of the gradient in spherical coordinates are
00278           // r: dV/dr
00279           // theta: 1/r * dV/dtheta
00280           // lambda: 1/(r*u) * dV/dlambda
00281           vrc =   - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);
00282           vtc =     qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);
00283           vlc = qs / u * (      A * (cl * vlc + sl * vls) + B * vlc2);
00284         }
00285       }
00286     }
00287 
00288     if (gradp) {
00289       // Rotate into cartesian (geocentric) coordinates
00290       gradx = cl * (u * vrc + t * vtc) - sl * vlc;
00291       grady = sl * (u * vrc + t * vtc) + cl * vlc;
00292       gradz =       t * vrc - u * vtc            ;
00293     }
00294     return vc;
00295   }
00296 
00297   template<bool gradp, SphericalEngine::normalization norm, int L>
00298   CircularEngine SphericalEngine::Circle(const coeff c[], const real f[],
00299                                          real p, real z, real a) {
00300 
00301     STATIC_ASSERT(L > 0, "L must be positive");
00302     STATIC_ASSERT(norm == FULL || norm == SCHMIDT, "Unknown normalization");
00303     int N = c[0].nmx(), M = c[0].mmx();
00304 
00305     real
00306       r = Math::hypot(z, p),
00307       t = r ? z / r : 0,            // cos(theta); at origin, pick theta = pi/2
00308       u = r ? max(p / r, eps_) : 1, // sin(theta); but avoid the pole
00309       q = a / r;
00310     real
00311       q2 = Math::sq(q),
00312       tu = t / u;
00313     CircularEngine circ(M, gradp, norm, a, r, u, t);
00314     int k[L];
00315     for (int m = M; m >= 0; --m) {   // m = M .. 0
00316       // Initialize inner sum
00317       real wc  = 0, wc2  = 0, ws  = 0, ws2  = 0; // w [N - m + 1], w [N - m + 2]
00318       real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0; // wr[N - m + 1], wr[N - m + 2]
00319       real wtc = 0, wtc2 = 0, wts = 0, wts2 = 0; // wt[N - m + 1], wt[N - m + 2]
00320       for (int l = 0; l < L; ++l)
00321         k[l] = c[l].index(N, m) + 1;
00322       for (int n = N; n >= m; --n) {             // n = N .. m; l = N - m .. 0
00323         real w, A, Ax, B, R;    // alpha[l], beta[l + 1]
00324         switch (norm) {
00325         case FULL:
00326           w = root_[2 * n + 1] / (root_[n - m + 1] * root_[n + m + 1]);
00327           Ax = q * w * root_[2 * n + 3];
00328           A = t * Ax;
00329           B = - q2 * root_[2 * n + 5] /
00330             (w * root_[n - m + 2] * root_[n + m + 2]);
00331           break;
00332         case SCHMIDT:
00333           w = root_[n - m + 1] * root_[n + m + 1];
00334           Ax = q * (2 * n + 1) / w;
00335           A = t * Ax;
00336           B = - q2 * w / (root_[n - m + 2] * root_[n + m + 2]);
00337           break;
00338         }
00339         R = c[0].Cv(--k[0]);
00340         for (int l = 1; l < L; ++l)
00341           R += c[l].Cv(--k[l], n, m, f[l]);
00342         R *= scale_;
00343         w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
00344         if (gradp) {
00345           w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
00346           w = A * wtc + B * wtc2 -  u*Ax * wc2; wtc2 = wtc; wtc = w;
00347         }
00348         if (m) {
00349           R = c[0].Sv(k[0]);
00350           for (int l = 1; l < L; ++l)
00351             R += c[l].Sv(k[l], n, m, f[l]);
00352           R *= scale_;
00353           w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
00354           if (gradp) {
00355             w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
00356             w = A * wts + B * wts2 -  u*Ax * ws2; wts2 = wts; wts = w;
00357           }
00358         }
00359       }
00360       if (!gradp)
00361         circ.SetCoeff(m, wc, ws);
00362       else {
00363         // Include the terms Sc[m] * P'[m,m](t) and  Ss[m] * P'[m,m](t)
00364         wtc += m * tu * wc; wts += m * tu * ws;
00365         circ.SetCoeff(m, wc, ws, wrc, wrs, wtc, wts);
00366       }
00367     }
00368 
00369     return circ;
00370   }
00371 
00372   void SphericalEngine::RootTable(int N) {
00373     // Need square roots up to max(2 * N + 5, 15).
00374     int L = max(2 * N + 5, 15) + 1, oldL = int(root_.size());
00375     if (oldL >= L)
00376       return;
00377     root_.resize(L);
00378     for (int l = oldL; l < L; ++l)
00379       root_[l] = sqrt(real(l));
00380   }
00381 
00382   void SphericalEngine::coeff::readcoeffs(std::istream& stream, int& N, int& M,
00383                                           std::vector<real>& C,
00384                                           std::vector<real>& S) {
00385     int nm[2];
00386     Utility::readarray<int, int, false>(stream, nm, 2);
00387     N = nm[0]; M = nm[1];
00388     if (!(N >= M && M >= -1 && N * M >= 0))
00389       // The last condition is that M = -1 implies N = -1 and vice versa.
00390       throw GeographicErr("Bad degree and order " +
00391                           Utility::str(N) + " " + Utility::str(M));
00392     C.resize(SphericalEngine::coeff::Csize(N, M));
00393     Utility::readarray<double, real, false>(stream, C);
00394     S.resize(SphericalEngine::coeff::Ssize(N, M));
00395     Utility::readarray<double, real, false>(stream, S);
00396     return;
00397   }
00398 
00399   /// \cond SKIP
00400   template
00401   Math::real SphericalEngine::Value<true, SphericalEngine::FULL, 1>
00402   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00403   template
00404   Math::real SphericalEngine::Value<false, SphericalEngine::FULL, 1>
00405   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00406   template
00407   Math::real SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
00408   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00409   template
00410   Math::real SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
00411   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00412 
00413   template
00414   Math::real SphericalEngine::Value<true, SphericalEngine::FULL, 2>
00415   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00416   template
00417   Math::real SphericalEngine::Value<false, SphericalEngine::FULL, 2>
00418   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00419   template
00420   Math::real SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
00421   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00422   template
00423   Math::real SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
00424   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00425 
00426   template
00427   Math::real SphericalEngine::Value<true, SphericalEngine::FULL, 3>
00428   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00429   template
00430   Math::real SphericalEngine::Value<false, SphericalEngine::FULL, 3>
00431   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00432   template
00433   Math::real SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
00434   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00435   template
00436   Math::real SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
00437   (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00438 
00439   template
00440   CircularEngine SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
00441   (const coeff[], const real[], real, real, real);
00442   template
00443   CircularEngine SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
00444   (const coeff[], const real[], real, real, real);
00445   template
00446   CircularEngine SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
00447   (const coeff[], const real[], real, real, real);
00448   template
00449   CircularEngine SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
00450   (const coeff[], const real[], real, real, real);
00451 
00452   template
00453   CircularEngine SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
00454   (const coeff[], const real[], real, real, real);
00455   template
00456   CircularEngine SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
00457   (const coeff[], const real[], real, real, real);
00458   template
00459   CircularEngine SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
00460   (const coeff[], const real[], real, real, real);
00461   template
00462   CircularEngine SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
00463   (const coeff[], const real[], real, real, real);
00464 
00465   template
00466   CircularEngine SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
00467   (const coeff[], const real[], real, real, real);
00468   template
00469   CircularEngine SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
00470   (const coeff[], const real[], real, real, real);
00471   template
00472   CircularEngine SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
00473   (const coeff[], const real[], real, real, real);
00474   template
00475   CircularEngine SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
00476   (const coeff[], const real[], real, real, real);
00477   /// \endcond
00478 
00479 } // namespace GeographicLib