137 #if defined(_MSC_VER)
140 # pragma warning (disable: 4127 4701)
147 const vector<Math::real> SphericalEngine::Z_(0);
148 vector<Math::real> SphericalEngine::root_(0);
150 template<
bool gradp, SphericalEngine::normalization norm,
int L>
152 real x, real y, real z, real a,
153 real& gradx, real& grady, real& gradz)
155 GEOGRAPHICLIB_STATIC_ASSERT(L > 0,
"L must be positive");
156 GEOGRAPHICLIB_STATIC_ASSERT(norm == FULL || norm == SCHMIDT,
157 "Unknown normalization");
158 int N = c[0].
nmx(), M = c[0].
mmx();
166 u = r ? max(p / r, eps()) : 1,
174 real vc = 0, vc2 = 0, vs = 0, vs2 = 0;
177 real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0;
178 real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0;
179 real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0;
181 for (
int m = M; m >= 0; --m) {
183 real wc = 0, wc2 = 0, ws = 0, ws2 = 0;
184 real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0;
185 real wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
186 for (
int l = 0; l < L; ++l)
187 k[l] = c[l].index(N, m) + 1;
188 for (
int n = N; n >= m; --n) {
192 w = root_[2 * n + 1] / (root_[n - m + 1] * root_[n + m + 1]);
193 Ax = q * w * root_[2 * n + 3];
195 B = - q2 * root_[2 * n + 5] /
196 (w * root_[n - m + 2] * root_[n + m + 2]);
199 w = root_[n - m + 1] * root_[n + m + 1];
200 Ax = q * (2 * n + 1) / w;
202 B = - q2 * w / (root_[n - m + 2] * root_[n + m + 2]);
207 for (
int l = 1; l < L; ++l)
208 R += c[l].Cv(--k[l], n, m, f[l]);
210 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
212 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
213 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
217 for (
int l = 1; l < L; ++l)
218 R += c[l].Sv(k[l], n, m, f[l]);
220 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
222 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
223 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
233 v = root_[2] * root_[2 * m + 3] / root_[m + 1];
235 B = - v * root_[2 * m + 5] / (root_[8] * root_[m + 2]) * uq2;
238 v = root_[2] * root_[2 * m + 1] / root_[m + 1];
240 B = - v * root_[2 * m + 3] / (root_[8] * root_[m + 2]) * uq2;
244 v = A * vc + B * vc2 + wc ; vc2 = vc ; vc = v;
245 v = A * vs + B * vs2 + ws ; vs2 = vs ; vs = v;
248 wtc += m * tu * wc; wts += m * tu * ws;
249 v = A * vrc + B * vrc2 + wrc; vrc2 = vrc; vrc = v;
250 v = A * vrs + B * vrs2 + wrs; vrs2 = vrs; vrs = v;
251 v = A * vtc + B * vtc2 + wtc; vtc2 = vtc; vtc = v;
252 v = A * vts + B * vts2 + wts; vts2 = vts; vts = v;
253 v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;
254 v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;
261 B = - root_[15]/2 * uq2;
265 B = - root_[3]/2 * uq2;
270 vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);
277 vrc = - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);
278 vtc = qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);
279 vlc = qs / u * ( A * (cl * vlc + sl * vls) + B * vlc2);
286 gradx = cl * (u * vrc + t * vtc) - sl * vlc;
287 grady = sl * (u * vrc + t * vtc) + cl * vlc;
288 gradz = t * vrc - u * vtc ;
293 template<
bool gradp, SphericalEngine::normalization norm,
int L>
295 real p, real z, real a) {
297 GEOGRAPHICLIB_STATIC_ASSERT(L > 0,
"L must be positive");
298 GEOGRAPHICLIB_STATIC_ASSERT(norm == FULL || norm == SCHMIDT,
299 "Unknown normalization");
300 int N = c[0].
nmx(), M = c[0].
mmx();
305 u = r ? max(p / r, eps()) : 1,
312 for (
int m = M; m >= 0; --m) {
314 real wc = 0, wc2 = 0, ws = 0, ws2 = 0;
315 real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0;
316 real wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
317 for (
int l = 0; l < L; ++l)
318 k[l] = c[l].index(N, m) + 1;
319 for (
int n = N; n >= m; --n) {
323 w = root_[2 * n + 1] / (root_[n - m + 1] * root_[n + m + 1]);
324 Ax = q * w * root_[2 * n + 3];
326 B = - q2 * root_[2 * n + 5] /
327 (w * root_[n - m + 2] * root_[n + m + 2]);
330 w = root_[n - m + 1] * root_[n + m + 1];
331 Ax = q * (2 * n + 1) / w;
333 B = - q2 * w / (root_[n - m + 2] * root_[n + m + 2]);
338 for (
int l = 1; l < L; ++l)
339 R += c[l].Cv(--k[l], n, m, f[l]);
341 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
343 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
344 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
348 for (
int l = 1; l < L; ++l)
349 R += c[l].Sv(k[l], n, m, f[l]);
351 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
353 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
354 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
359 circ.SetCoeff(m, wc, ws);
362 wtc += m * tu * wc; wts += m * tu * ws;
363 circ.SetCoeff(m, wc, ws, wrc, wrs, wtc, wts);
372 int L = max(2 * N + 5, 15) + 1, oldL = int(root_.size());
376 for (
int l = oldL; l < L; ++l)
377 root_[l] = sqrt(
real(l));
381 std::vector<real>& C,
382 std::vector<real>& S) {
384 Utility::readarray<int, int, false>(stream, nm, 2);
385 N = nm[0]; M = nm[1];
386 if (!(N >= M && M >= -1 && N * M >= 0))
391 Utility::readarray<double, real, false>(stream, C);
393 Utility::readarray<double, real, false>(stream, S);
399 SphericalEngine::Value<true, SphericalEngine::FULL, 1>
402 SphericalEngine::Value<false, SphericalEngine::FULL, 1>
405 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
408 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
412 SphericalEngine::Value<true, SphericalEngine::FULL, 2>
415 SphericalEngine::Value<false, SphericalEngine::FULL, 2>
418 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
421 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
425 SphericalEngine::Value<true, SphericalEngine::FULL, 3>
428 SphericalEngine::Value<false, SphericalEngine::FULL, 3>
431 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
434 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
438 SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
441 SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
444 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
447 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
451 SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
454 SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
457 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
460 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
464 SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
467 SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
470 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
473 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
#define GEOGRAPHICLIB_EXPORT
GeographicLib::Math::real real
Math::real Cv(int k) const
static int Ssize(int N, int M)
Header for GeographicLib::Utility class.
static CircularEngine Circle(const coeff c[], const real f[], real p, real z, real a)
static void RootTable(int N)
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S)
Package up coefficients for SphericalEngine.
Math::real Sv(int k) const
Namespace for GeographicLib.
static std::string str(T x, int p=-1)
Header for GeographicLib::CircularEngine class.
Spherical harmonic sums for a circle.
Exception handling for GeographicLib.
static int Csize(int N, int M)
static Math::real Value(const coeff c[], const real f[], real x, real y, real z, real a, real &gradx, real &grady, real &gradz)
Header for GeographicLib::SphericalEngine class.