00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 #include <GeographicLib/SphericalEngine.hpp>
00134 #include <GeographicLib/CircularEngine.hpp>
00135 #include <GeographicLib/Utility.hpp>
00136
00137 #if defined(_MSC_VER)
00138
00139
00140 # pragma warning (disable: 4127 4701)
00141 #endif
00142
00143 namespace GeographicLib {
00144
00145 using namespace std;
00146
00147 const vector<Math::real> SphericalEngine::Z_(0);
00148 vector<Math::real> SphericalEngine::root_(0);
00149
00150 template<bool gradp, SphericalEngine::normalization norm, int L>
00151 Math::real SphericalEngine::Value(const coeff c[], const real f[],
00152 real x, real y, real z, real a,
00153 real& gradx, real& grady, real& gradz)
00154 {
00155 GEOGRAPHICLIB_STATIC_ASSERT(L > 0, "L must be positive");
00156 GEOGRAPHICLIB_STATIC_ASSERT(norm == FULL || norm == SCHMIDT,
00157 "Unknown normalization");
00158 int N = c[0].nmx(), M = c[0].mmx();
00159
00160 real
00161 p = Math::hypot(x, y),
00162 cl = p ? x / p : 1,
00163 sl = p ? y / p : 0,
00164 r = Math::hypot(z, p),
00165 t = r ? z / r : 0,
00166 u = r ? max(p / r, eps()) : 1,
00167 q = a / r;
00168 real
00169 q2 = Math::sq(q),
00170 uq = u * q,
00171 uq2 = Math::sq(uq),
00172 tu = t / u;
00173
00174 real vc = 0, vc2 = 0, vs = 0, vs2 = 0;
00175
00176
00177 real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0;
00178 real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0;
00179 real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0;
00180 int k[L];
00181 for (int m = M; m >= 0; --m) {
00182
00183 real wc = 0, wc2 = 0, ws = 0, ws2 = 0;
00184 real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0;
00185 real wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
00186 for (int l = 0; l < L; ++l)
00187 k[l] = c[l].index(N, m) + 1;
00188 for (int n = N; n >= m; --n) {
00189 real w, A, Ax, B, R;
00190 switch (norm) {
00191 case FULL:
00192 w = root_[2 * n + 1] / (root_[n - m + 1] * root_[n + m + 1]);
00193 Ax = q * w * root_[2 * n + 3];
00194 A = t * Ax;
00195 B = - q2 * root_[2 * n + 5] /
00196 (w * root_[n - m + 2] * root_[n + m + 2]);
00197 break;
00198 case SCHMIDT:
00199 w = root_[n - m + 1] * root_[n + m + 1];
00200 Ax = q * (2 * n + 1) / w;
00201 A = t * Ax;
00202 B = - q2 * w / (root_[n - m + 2] * root_[n + m + 2]);
00203 break;
00204 default: break;
00205 }
00206 R = c[0].Cv(--k[0]);
00207 for (int l = 1; l < L; ++l)
00208 R += c[l].Cv(--k[l], n, m, f[l]);
00209 R *= scale();
00210 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
00211 if (gradp) {
00212 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
00213 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
00214 }
00215 if (m) {
00216 R = c[0].Sv(k[0]);
00217 for (int l = 1; l < L; ++l)
00218 R += c[l].Sv(k[l], n, m, f[l]);
00219 R *= scale();
00220 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
00221 if (gradp) {
00222 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
00223 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
00224 }
00225 }
00226 }
00227
00228
00229 if (m) {
00230 real v, A, B;
00231 switch (norm) {
00232 case FULL:
00233 v = root_[2] * root_[2 * m + 3] / root_[m + 1];
00234 A = cl * v * uq;
00235 B = - v * root_[2 * m + 5] / (root_[8] * root_[m + 2]) * uq2;
00236 break;
00237 case SCHMIDT:
00238 v = root_[2] * root_[2 * m + 1] / root_[m + 1];
00239 A = cl * v * uq;
00240 B = - v * root_[2 * m + 3] / (root_[8] * root_[m + 2]) * uq2;
00241 break;
00242 default: break;
00243 }
00244 v = A * vc + B * vc2 + wc ; vc2 = vc ; vc = v;
00245 v = A * vs + B * vs2 + ws ; vs2 = vs ; vs = v;
00246 if (gradp) {
00247
00248 wtc += m * tu * wc; wts += m * tu * ws;
00249 v = A * vrc + B * vrc2 + wrc; vrc2 = vrc; vrc = v;
00250 v = A * vrs + B * vrs2 + wrs; vrs2 = vrs; vrs = v;
00251 v = A * vtc + B * vtc2 + wtc; vtc2 = vtc; vtc = v;
00252 v = A * vts + B * vts2 + wts; vts2 = vts; vts = v;
00253 v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;
00254 v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;
00255 }
00256 } else {
00257 real A, B, qs;
00258 switch (norm) {
00259 case FULL:
00260 A = root_[3] * uq;
00261 B = - root_[15]/2 * uq2;
00262 break;
00263 case SCHMIDT:
00264 A = uq;
00265 B = - root_[3]/2 * uq2;
00266 break;
00267 default: break;
00268 }
00269 qs = q / scale();
00270 vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);
00271 if (gradp) {
00272 qs /= r;
00273
00274
00275
00276
00277 vrc = - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);
00278 vtc = qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);
00279 vlc = qs / u * ( A * (cl * vlc + sl * vls) + B * vlc2);
00280 }
00281 }
00282 }
00283
00284 if (gradp) {
00285
00286 gradx = cl * (u * vrc + t * vtc) - sl * vlc;
00287 grady = sl * (u * vrc + t * vtc) + cl * vlc;
00288 gradz = t * vrc - u * vtc ;
00289 }
00290 return vc;
00291 }
00292
00293 template<bool gradp, SphericalEngine::normalization norm, int L>
00294 CircularEngine SphericalEngine::Circle(const coeff c[], const real f[],
00295 real p, real z, real a) {
00296
00297 GEOGRAPHICLIB_STATIC_ASSERT(L > 0, "L must be positive");
00298 GEOGRAPHICLIB_STATIC_ASSERT(norm == FULL || norm == SCHMIDT,
00299 "Unknown normalization");
00300 int N = c[0].nmx(), M = c[0].mmx();
00301
00302 real
00303 r = Math::hypot(z, p),
00304 t = r ? z / r : 0,
00305 u = r ? max(p / r, eps()) : 1,
00306 q = a / r;
00307 real
00308 q2 = Math::sq(q),
00309 tu = t / u;
00310 CircularEngine circ(M, gradp, norm, a, r, u, t);
00311 int k[L];
00312 for (int m = M; m >= 0; --m) {
00313
00314 real wc = 0, wc2 = 0, ws = 0, ws2 = 0;
00315 real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0;
00316 real wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
00317 for (int l = 0; l < L; ++l)
00318 k[l] = c[l].index(N, m) + 1;
00319 for (int n = N; n >= m; --n) {
00320 real w, A, Ax, B, R;
00321 switch (norm) {
00322 case FULL:
00323 w = root_[2 * n + 1] / (root_[n - m + 1] * root_[n + m + 1]);
00324 Ax = q * w * root_[2 * n + 3];
00325 A = t * Ax;
00326 B = - q2 * root_[2 * n + 5] /
00327 (w * root_[n - m + 2] * root_[n + m + 2]);
00328 break;
00329 case SCHMIDT:
00330 w = root_[n - m + 1] * root_[n + m + 1];
00331 Ax = q * (2 * n + 1) / w;
00332 A = t * Ax;
00333 B = - q2 * w / (root_[n - m + 2] * root_[n + m + 2]);
00334 break;
00335 default: break;
00336 }
00337 R = c[0].Cv(--k[0]);
00338 for (int l = 1; l < L; ++l)
00339 R += c[l].Cv(--k[l], n, m, f[l]);
00340 R *= scale();
00341 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
00342 if (gradp) {
00343 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
00344 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
00345 }
00346 if (m) {
00347 R = c[0].Sv(k[0]);
00348 for (int l = 1; l < L; ++l)
00349 R += c[l].Sv(k[l], n, m, f[l]);
00350 R *= scale();
00351 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
00352 if (gradp) {
00353 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
00354 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
00355 }
00356 }
00357 }
00358 if (!gradp)
00359 circ.SetCoeff(m, wc, ws);
00360 else {
00361
00362 wtc += m * tu * wc; wts += m * tu * ws;
00363 circ.SetCoeff(m, wc, ws, wrc, wrs, wtc, wts);
00364 }
00365 }
00366
00367 return circ;
00368 }
00369
00370 void SphericalEngine::RootTable(int N) {
00371
00372 int L = max(2 * N + 5, 15) + 1, oldL = int(root_.size());
00373 if (oldL >= L)
00374 return;
00375 root_.resize(L);
00376 for (int l = oldL; l < L; ++l)
00377 root_[l] = sqrt(real(l));
00378 }
00379
00380 void SphericalEngine::coeff::readcoeffs(std::istream& stream, int& N, int& M,
00381 std::vector<real>& C,
00382 std::vector<real>& S) {
00383 int nm[2];
00384 Utility::readarray<int, int, false>(stream, nm, 2);
00385 N = nm[0]; M = nm[1];
00386 if (!(N >= M && M >= -1 && N * M >= 0))
00387
00388 throw GeographicErr("Bad degree and order " +
00389 Utility::str(N) + " " + Utility::str(M));
00390 C.resize(SphericalEngine::coeff::Csize(N, M));
00391 Utility::readarray<double, real, false>(stream, C);
00392 S.resize(SphericalEngine::coeff::Ssize(N, M));
00393 Utility::readarray<double, real, false>(stream, S);
00394 return;
00395 }
00396
00397
00398 template Math::real GEOGRAPHICLIB_EXPORT
00399 SphericalEngine::Value<true, SphericalEngine::FULL, 1>
00400 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00401 template Math::real GEOGRAPHICLIB_EXPORT
00402 SphericalEngine::Value<false, SphericalEngine::FULL, 1>
00403 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00404 template Math::real GEOGRAPHICLIB_EXPORT
00405 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
00406 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00407 template Math::real GEOGRAPHICLIB_EXPORT
00408 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
00409 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00410
00411 template Math::real GEOGRAPHICLIB_EXPORT
00412 SphericalEngine::Value<true, SphericalEngine::FULL, 2>
00413 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00414 template Math::real GEOGRAPHICLIB_EXPORT
00415 SphericalEngine::Value<false, SphericalEngine::FULL, 2>
00416 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00417 template Math::real GEOGRAPHICLIB_EXPORT
00418 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
00419 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00420 template Math::real GEOGRAPHICLIB_EXPORT
00421 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
00422 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00423
00424 template Math::real GEOGRAPHICLIB_EXPORT
00425 SphericalEngine::Value<true, SphericalEngine::FULL, 3>
00426 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00427 template Math::real GEOGRAPHICLIB_EXPORT
00428 SphericalEngine::Value<false, SphericalEngine::FULL, 3>
00429 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00430 template Math::real GEOGRAPHICLIB_EXPORT
00431 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
00432 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00433 template Math::real GEOGRAPHICLIB_EXPORT
00434 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
00435 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
00436
00437 template CircularEngine GEOGRAPHICLIB_EXPORT
00438 SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
00439 (const coeff[], const real[], real, real, real);
00440 template CircularEngine GEOGRAPHICLIB_EXPORT
00441 SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
00442 (const coeff[], const real[], real, real, real);
00443 template CircularEngine GEOGRAPHICLIB_EXPORT
00444 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
00445 (const coeff[], const real[], real, real, real);
00446 template CircularEngine GEOGRAPHICLIB_EXPORT
00447 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
00448 (const coeff[], const real[], real, real, real);
00449
00450 template CircularEngine GEOGRAPHICLIB_EXPORT
00451 SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
00452 (const coeff[], const real[], real, real, real);
00453 template CircularEngine GEOGRAPHICLIB_EXPORT
00454 SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
00455 (const coeff[], const real[], real, real, real);
00456 template CircularEngine GEOGRAPHICLIB_EXPORT
00457 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
00458 (const coeff[], const real[], real, real, real);
00459 template CircularEngine GEOGRAPHICLIB_EXPORT
00460 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
00461 (const coeff[], const real[], real, real, real);
00462
00463 template CircularEngine GEOGRAPHICLIB_EXPORT
00464 SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
00465 (const coeff[], const real[], real, real, real);
00466 template CircularEngine GEOGRAPHICLIB_EXPORT
00467 SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
00468 (const coeff[], const real[], real, real, real);
00469 template CircularEngine GEOGRAPHICLIB_EXPORT
00470 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
00471 (const coeff[], const real[], real, real, real);
00472 template CircularEngine GEOGRAPHICLIB_EXPORT
00473 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
00474 (const coeff[], const real[], real, real, real);
00475
00476
00477 }