GeographicLib  1.21
CircularEngine.cpp
Go to the documentation of this file.
00001 /**
00002  * \file CircularEngine.cpp
00003  * \brief Implementation for GeographicLib::CircularEngine 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 
00010 #include <GeographicLib/CircularEngine.hpp>
00011 #include <limits>
00012 
00013 #define GEOGRAPHICLIB_CIRCULARENGINE_CPP \
00014   "$Id: bdd0d21aa34063706e4042410f06bb0f7844fea9 $"
00015 
00016 RCSID_DECL(GEOGRAPHICLIB_CIRCULARENGINE_CPP)
00017 RCSID_DECL(GEOGRAPHICLIB_CIRCULARENGINE_HPP)
00018 
00019 namespace GeographicLib {
00020 
00021   using namespace std;
00022 
00023   Math::real CircularEngine::Value(bool gradp, real cl, real sl,
00024                                    real& gradx, real& grady, real& gradz)
00025     const throw() {
00026     gradp = _gradp && gradp;
00027     const vector<real>& root_( SphericalEngine::root_ );
00028 
00029     // Initialize outer sum
00030     real vc  = 0, vc2  = 0, vs  = 0, vs2  = 0;   // v [N + 1], v [N + 2]
00031     // vr, vt, vl and similar w variable accumulate the sums for the
00032     // derivatives wrt r, theta, and lambda, respectively.
00033     real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0;   // vr[N + 1], vr[N + 2]
00034     real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0;   // vt[N + 1], vt[N + 2]
00035     real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0;   // vl[N + 1], vl[N + 2]
00036     for (int m = _M; m >= 0; --m) {   // m = M .. 0
00037       // Now Sc[m] = wc, Ss[m] = ws
00038       // Sc'[m] = wtc, Ss'[m] = wtc
00039       if (m) {
00040         real v, A, B;           // alpha[m], beta[m + 1]
00041         switch (_norm) {
00042         case FULL:
00043           v = root_[2] * root_[2 * m + 3] / root_[m + 1];
00044           A = cl * v * _uq;
00045           B = - v * root_[2 * m + 5] / (root_[8] * root_[m + 2]) * _uq2;
00046           break;
00047         case SCHMIDT:
00048           v = root_[2] * root_[2 * m + 1] / root_[m + 1];
00049           A = cl * v * _uq;
00050           B = - v * root_[2 * m + 3] / (root_[8] * root_[m + 2]) * _uq2;
00051           break;
00052         default:
00053           A = B = 0;
00054         }
00055         v = A * vc  + B * vc2  +  _wc[m] ; vc2  = vc ; vc  = v;
00056         v = A * vs  + B * vs2  +  _ws[m] ; vs2  = vs ; vs  = v;
00057         if (gradp) {
00058           v = A * vrc + B * vrc2 +  _wrc[m]; vrc2 = vrc; vrc = v;
00059           v = A * vrs + B * vrs2 +  _wrs[m]; vrs2 = vrs; vrs = v;
00060           v = A * vtc + B * vtc2 +  _wtc[m]; vtc2 = vtc; vtc = v;
00061           v = A * vts + B * vts2 +  _wts[m]; vts2 = vts; vts = v;
00062           v = A * vlc + B * vlc2 + m*_ws[m]; vlc2 = vlc; vlc = v;
00063           v = A * vls + B * vls2 - m*_wc[m]; vls2 = vls; vls = v;
00064         }
00065       } else {
00066         real A, B, qs;
00067         switch (_norm) {
00068         case FULL:
00069           A = root_[3] * _uq;       // F[1]/(q*cl) or F[1]/(q*sl)
00070           B = - root_[15]/2 * _uq2; // beta[1]/q
00071           break;
00072         case SCHMIDT:
00073           A = _uq;
00074           B = - root_[3]/2 * _uq2;
00075           break;
00076         default:
00077           A = B = 0;
00078         }
00079         qs = _q / SphericalEngine::scale_;
00080         vc = qs * (_wc[m] + A * (cl * vc + sl * vs ) + B * vc2);
00081         if (gradp) {
00082           qs /= _r;
00083           // The components of the gradient in circular coordinates are
00084           // r: dV/dr
00085           // theta: 1/r * dV/dtheta
00086           // lambda: 1/(r*u) * dV/dlambda
00087           vrc =    - qs * (_wrc[m] + A * (cl * vrc + sl * vrs) + B * vrc2);
00088           vtc =      qs * (_wtc[m] + A * (cl * vtc + sl * vts) + B * vtc2);
00089           vlc = qs / _u * (          A * (cl * vlc + sl * vls) + B * vlc2);
00090         }
00091       }
00092     }
00093 
00094     if (gradp) {
00095       // Rotate into cartesian (geocentric) coordinates
00096       gradx = cl * (_u * vrc + _t * vtc) - sl * vlc;
00097       grady = sl * (_u * vrc + _t * vtc) + cl * vlc;
00098       gradz =           _t * vrc - _u * vtc                ;
00099     }
00100     return vc;
00101   }
00102 
00103 } // namespace GeographicLib