00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/CircularEngine.hpp>
00011
00012 namespace GeographicLib {
00013
00014 using namespace std;
00015
00016 Math::real CircularEngine::Value(bool gradp, real cl, real sl,
00017 real& gradx, real& grady, real& gradz)
00018 const {
00019 gradp = _gradp && gradp;
00020 const vector<real>& root_( SphericalEngine::root_ );
00021
00022
00023 real vc = 0, vc2 = 0, vs = 0, vs2 = 0;
00024
00025
00026 real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0;
00027 real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0;
00028 real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0;
00029 for (int m = _M; m >= 0; --m) {
00030
00031
00032 if (m) {
00033 real v, A, B;
00034 switch (_norm) {
00035 case FULL:
00036 v = root_[2] * root_[2 * m + 3] / root_[m + 1];
00037 A = cl * v * _uq;
00038 B = - v * root_[2 * m + 5] / (root_[8] * root_[m + 2]) * _uq2;
00039 break;
00040 case SCHMIDT:
00041 v = root_[2] * root_[2 * m + 1] / root_[m + 1];
00042 A = cl * v * _uq;
00043 B = - v * root_[2 * m + 3] / (root_[8] * root_[m + 2]) * _uq2;
00044 break;
00045 default:
00046 A = B = 0;
00047 }
00048 v = A * vc + B * vc2 + _wc[m] ; vc2 = vc ; vc = v;
00049 v = A * vs + B * vs2 + _ws[m] ; vs2 = vs ; vs = v;
00050 if (gradp) {
00051 v = A * vrc + B * vrc2 + _wrc[m]; vrc2 = vrc; vrc = v;
00052 v = A * vrs + B * vrs2 + _wrs[m]; vrs2 = vrs; vrs = v;
00053 v = A * vtc + B * vtc2 + _wtc[m]; vtc2 = vtc; vtc = v;
00054 v = A * vts + B * vts2 + _wts[m]; vts2 = vts; vts = v;
00055 v = A * vlc + B * vlc2 + m*_ws[m]; vlc2 = vlc; vlc = v;
00056 v = A * vls + B * vls2 - m*_wc[m]; vls2 = vls; vls = v;
00057 }
00058 } else {
00059 real A, B, qs;
00060 switch (_norm) {
00061 case FULL:
00062 A = root_[3] * _uq;
00063 B = - root_[15]/2 * _uq2;
00064 break;
00065 case SCHMIDT:
00066 A = _uq;
00067 B = - root_[3]/2 * _uq2;
00068 break;
00069 default:
00070 A = B = 0;
00071 }
00072 qs = _q / SphericalEngine::scale();
00073 vc = qs * (_wc[m] + A * (cl * vc + sl * vs ) + B * vc2);
00074 if (gradp) {
00075 qs /= _r;
00076
00077
00078
00079
00080 vrc = - qs * (_wrc[m] + A * (cl * vrc + sl * vrs) + B * vrc2);
00081 vtc = qs * (_wtc[m] + A * (cl * vtc + sl * vts) + B * vtc2);
00082 vlc = qs / _u * ( A * (cl * vlc + sl * vls) + B * vlc2);
00083 }
00084 }
00085 }
00086
00087 if (gradp) {
00088
00089 gradx = cl * (_u * vrc + _t * vtc) - sl * vlc;
00090 grady = sl * (_u * vrc + _t * vtc) + cl * vlc;
00091 gradz = _t * vrc - _u * vtc ;
00092 }
00093 return vc;
00094 }
00095
00096 }