139# pragma warning (disable: 4701)
146 vector<Math::real>& SphericalEngine::sqrttable() {
147 static vector<real> sqrttable(0);
151 template<
bool gradp, SphericalEngine::normalization norm,
int L>
153 real x, real y, real z, real a,
154 real& gradx, real& grady, real& gradz)
156 static_assert(L > 0,
"L must be positive");
157 static_assert(norm ==
FULL || norm ==
SCHMIDT,
"Unknown normalization");
158 int N = c[0].
nmx(), M = c[0].
mmx();
162 cl = p != 0 ? x / p : 1,
163 sl = p != 0 ? y / p : 0,
165 t = r != 0 ? z / r : 0,
166 u = r != 0 ? fmax(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 const vector<real>& root( sqrttable() );
182 for (
int m = M; m >= 0; --m) {
185 wc = 0, wc2 = 0, ws = 0, ws2 = 0,
186 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0,
187 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
188 for (
int l = 0; l < L; ++l)
189 k[l] = c[l].index(N, m) + 1;
190 for (
int n = N; n >= m; --n) {
194 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);
195 Ax = q * w * root[2 * n + 3];
197 B = - q2 * root[2 * n + 5] /
198 (w * root[n - m + 2] * root[n + m + 2]);
201 w = root[n - m + 1] * root[n + m + 1];
202 Ax = q * (2 * n + 1) / w;
204 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);
209 for (
int l = 1; l < L; ++l)
210 R += c[l].Cv(--k[l], n, m, f[l]);
212 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
214 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
215 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
219 for (
int l = 1; l < L; ++l)
220 R += c[l].Sv(k[l], n, m, f[l]);
222 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
224 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
225 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
235 v = root[2] * root[2 * m + 3] / root[m + 1];
237 B = - v * root[2 * m + 5] / (root[8] * root[m + 2]) * uq2;
240 v = root[2] * root[2 * m + 1] / root[m + 1];
242 B = - v * root[2 * m + 3] / (root[8] * root[m + 2]) * uq2;
246 v = A * vc + B * vc2 + wc ; vc2 = vc ; vc = v;
247 v = A * vs + B * vs2 + ws ; vs2 = vs ; vs = v;
250 wtc += m * tu * wc; wts += m * tu * ws;
251 v = A * vrc + B * vrc2 + wrc; vrc2 = vrc; vrc = v;
252 v = A * vrs + B * vrs2 + wrs; vrs2 = vrs; vrs = v;
253 v = A * vtc + B * vtc2 + wtc; vtc2 = vtc; vtc = v;
254 v = A * vts + B * vts2 + wts; vts2 = vts; vts = v;
255 v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;
256 v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;
263 B = - root[15]/2 * uq2;
267 B = - root[3]/2 * uq2;
272 vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);
279 vrc = - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);
280 vtc = qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);
281 vlc = qs / u * ( A * (cl * vlc + sl * vls) + B * vlc2);
288 gradx = cl * (u * vrc + t * vtc) - sl * vlc;
289 grady = sl * (u * vrc + t * vtc) + cl * vlc;
290 gradz = t * vrc - u * vtc ;
295 template<
bool gradp, SphericalEngine::normalization norm,
int L>
297 real p, real z, real a) {
299 static_assert(L > 0,
"L must be positive");
300 static_assert(norm ==
FULL || norm ==
SCHMIDT,
"Unknown normalization");
301 int N = c[0].
nmx(), M = c[0].
mmx();
305 t = r != 0 ? z / r : 0,
306 u = r != 0 ? fmax(p / r, eps()) : 1,
313 const vector<real>& root( sqrttable() );
314 for (
int m = M; m >= 0; --m) {
317 wc = 0, wc2 = 0, ws = 0, ws2 = 0,
318 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0,
319 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
320 for (
int l = 0; l < L; ++l)
321 k[l] = c[l].index(N, m) + 1;
322 for (
int n = N; n >= m; --n) {
326 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);
327 Ax = q * w * root[2 * n + 3];
329 B = - q2 * root[2 * n + 5] /
330 (w * root[n - m + 2] * root[n + m + 2]);
333 w = root[n - m + 1] * root[n + m + 1];
334 Ax = q * (2 * n + 1) / w;
336 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);
341 for (
int l = 1; l < L; ++l)
342 R += c[l].Cv(--k[l], n, m, f[l]);
344 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
346 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
347 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
351 for (
int l = 1; l < L; ++l)
352 R += c[l].Sv(k[l], n, m, f[l]);
354 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
356 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
357 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
362 circ.SetCoeff(m, wc, ws);
365 wtc += m * tu * wc; wts += m * tu * ws;
366 circ.SetCoeff(m, wc, ws, wrc, wrs, wtc, wts);
375 vector<real>& root( sqrttable() );
376 int L = max(2 * N + 5, 15) + 1, oldL = int(root.size());
380 for (
int l = oldL; l < L; ++l)
381 root[l] = sqrt(real(l));
389 if (!((
N >= M && M >= 0) || (
N == -1 && M == -1)))
395 Utility::readarray<int, int, false>(stream, nm, 2);
396 int N0 = nm[0], M0 = nm[1];
397 if (!((N0 >= M0 && M0 >= 0) || (N0 == -1 && M0 == -1)))
401 N = truncate ? min(
N, N0) : N0;
402 M = truncate ? min(M, M0) : M0;
408 Utility::readarray<double, real, false>(stream, C);
409 if (skip) stream.seekg(streamoff(skip), ios::cur);
410 Utility::readarray<double, real, false>(stream, S);
411 if (skip) stream.seekg(streamoff(skip), ios::cur);
413 for (
int m = 0, k = 0; m <= M; ++m) {
414 Utility::readarray<double, real, false>(stream, &C[k],
N + 1 - m);
415 stream.seekg((N0 -
N) *
sizeof(
double), ios::cur);
418 if (skip) stream.seekg(streamoff(skip), ios::cur);
419 for (
int m = 1, k = 0; m <= M; ++m) {
420 Utility::readarray<double, real, false>(stream, &S[k],
N + 1 - m);
421 stream.seekg((N0 -
N) *
sizeof(
double), ios::cur);
424 if (skip) stream.seekg(streamoff(skip), ios::cur);
431 SphericalEngine::Value<true, SphericalEngine::FULL, 1>
432 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
434 SphericalEngine::Value<false, SphericalEngine::FULL, 1>
435 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
437 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
438 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
440 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
441 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
444 SphericalEngine::Value<true, SphericalEngine::FULL, 2>
445 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
447 SphericalEngine::Value<false, SphericalEngine::FULL, 2>
448 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
450 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
451 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
453 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
454 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
457 SphericalEngine::Value<true, SphericalEngine::FULL, 3>
458 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
460 SphericalEngine::Value<false, SphericalEngine::FULL, 3>
461 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
463 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
464 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
466 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
467 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
470 SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
471 (
const coeff[],
const real[], real, real, real);
473 SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
474 (
const coeff[],
const real[], real, real, real);
476 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
477 (
const coeff[],
const real[], real, real, real);
479 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
480 (
const coeff[],
const real[], real, real, real);
483 SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
484 (
const coeff[],
const real[], real, real, real);
486 SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
487 (
const coeff[],
const real[], real, real, real);
489 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
490 (
const coeff[],
const real[], real, real, real);
492 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
493 (
const coeff[],
const real[], real, real, real);
496 SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
497 (
const coeff[],
const real[], real, real, real);
499 SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
500 (
const coeff[],
const real[], real, real, real);
502 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
503 (
const coeff[],
const real[], real, real, real);
505 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
506 (
const coeff[],
const real[], real, real, real);
Header for GeographicLib::CircularEngine class.
#define GEOGRAPHICLIB_EXPORT
Header for GeographicLib::SphericalEngine class.
Header for GeographicLib::Utility class.
Spherical harmonic sums for a circle.
Exception handling for GeographicLib.
Package up coefficients for SphericalEngine.
Math::real Sv(int k) const
static int Csize(int N, int M)
static int Ssize(int N, int M)
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)
Math::real Cv(int k) const
static Math::real Value(const coeff c[], const real f[], real x, real y, real z, real a, real &gradx, real &grady, real &gradz)
static void RootTable(int N)
static CircularEngine Circle(const coeff c[], const real f[], real p, real z, real a)
static std::string str(T x, int p=-1)
Namespace for GeographicLib.