GeographicLib
1.21
|
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