137 #if defined(_MSC_VER)
140 # pragma warning (disable: 4127 4701)
147 vector<Math::real>& SphericalEngine::sqrttable() {
148 static vector<real> sqrttable(0);
152 template<
bool gradp, SphericalEngine::normalization norm,
int L>
154 real x, real y, real z, real a,
155 real& gradx, real& grady, real& gradz)
157 GEOGRAPHICLIB_STATIC_ASSERT(L > 0,
"L must be positive");
158 GEOGRAPHICLIB_STATIC_ASSERT(norm == FULL || norm == SCHMIDT,
159 "Unknown normalization");
160 int N = c[0].
nmx(), M = c[0].
mmx();
164 cl = p != 0 ? x / p : 1,
165 sl = p != 0 ? y / p : 0,
167 t = r != 0 ? z / r : 0,
168 u = r != 0 ? max(p / r, eps()) : 1,
176 real vc = 0, vc2 = 0, vs = 0, vs2 = 0;
179 real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0;
180 real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0;
181 real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0;
183 const vector<real>& root( sqrttable() );
184 for (
int m = M; m >= 0; --m) {
187 wc = 0, wc2 = 0, ws = 0, ws2 = 0,
188 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0,
189 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
190 for (
int l = 0; l < L; ++l)
191 k[l] = c[l].index(N, m) + 1;
192 for (
int n = N; n >= m; --n) {
196 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);
197 Ax = q * w * root[2 * n + 3];
199 B = - q2 * root[2 * n + 5] /
200 (w * root[n - m + 2] * root[n + m + 2]);
203 w = root[n - m + 1] * root[n + m + 1];
204 Ax = q * (2 * n + 1) / w;
206 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);
211 for (
int l = 1; l < L; ++l)
212 R += c[l].Cv(--k[l], n, m, f[l]);
214 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
216 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
217 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
221 for (
int l = 1; l < L; ++l)
222 R += c[l].Sv(k[l], n, m, f[l]);
224 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
226 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
227 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
237 v = root[2] * root[2 * m + 3] / root[m + 1];
239 B = - v * root[2 * m + 5] / (root[8] * root[m + 2]) * uq2;
242 v = root[2] * root[2 * m + 1] / root[m + 1];
244 B = - v * root[2 * m + 3] / (root[8] * root[m + 2]) * uq2;
248 v = A * vc + B * vc2 + wc ; vc2 = vc ; vc = v;
249 v = A * vs + B * vs2 + ws ; vs2 = vs ; vs = v;
252 wtc += m * tu * wc; wts += m * tu * ws;
253 v = A * vrc + B * vrc2 + wrc; vrc2 = vrc; vrc = v;
254 v = A * vrs + B * vrs2 + wrs; vrs2 = vrs; vrs = v;
255 v = A * vtc + B * vtc2 + wtc; vtc2 = vtc; vtc = v;
256 v = A * vts + B * vts2 + wts; vts2 = vts; vts = v;
257 v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;
258 v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;
265 B = - root[15]/2 * uq2;
269 B = - root[3]/2 * uq2;
274 vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);
281 vrc = - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);
282 vtc = qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);
283 vlc = qs / u * ( A * (cl * vlc + sl * vls) + B * vlc2);
290 gradx = cl * (u * vrc + t * vtc) - sl * vlc;
291 grady = sl * (u * vrc + t * vtc) + cl * vlc;
292 gradz = t * vrc - u * vtc ;
297 template<
bool gradp, SphericalEngine::normalization norm,
int L>
299 real p, real z, real a) {
301 GEOGRAPHICLIB_STATIC_ASSERT(L > 0,
"L must be positive");
302 GEOGRAPHICLIB_STATIC_ASSERT(norm == FULL || norm == SCHMIDT,
303 "Unknown normalization");
304 int N = c[0].
nmx(), M = c[0].
mmx();
308 t = r != 0 ? z / r : 0,
309 u = r != 0 ? max(p / r, eps()) : 1,
316 const vector<real>& root( sqrttable() );
317 for (
int m = M; m >= 0; --m) {
320 wc = 0, wc2 = 0, ws = 0, ws2 = 0,
321 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0,
322 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;
323 for (
int l = 0; l < L; ++l)
324 k[l] = c[l].index(N, m) + 1;
325 for (
int n = N; n >= m; --n) {
329 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);
330 Ax = q * w * root[2 * n + 3];
332 B = - q2 * root[2 * n + 5] /
333 (w * root[n - m + 2] * root[n + m + 2]);
336 w = root[n - m + 1] * root[n + m + 1];
337 Ax = q * (2 * n + 1) / w;
339 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);
344 for (
int l = 1; l < L; ++l)
345 R += c[l].Cv(--k[l], n, m, f[l]);
347 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
349 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
350 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
354 for (
int l = 1; l < L; ++l)
355 R += c[l].Sv(k[l], n, m, f[l]);
357 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
359 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
360 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
365 circ.SetCoeff(m, wc, ws);
368 wtc += m * tu * wc; wts += m * tu * ws;
369 circ.SetCoeff(m, wc, ws, wrc, wrs, wtc, wts);
378 vector<real>& root( sqrttable() );
379 int L = max(2 * N + 5, 15) + 1, oldL = int(root.size());
383 for (
int l = oldL; l < L; ++l)
384 root[l] = sqrt(
real(l));
388 std::vector<real>& C,
389 std::vector<real>& S,
392 if (!((N >= M && M >= 0) || (N == -1 && M == -1)))
398 Utility::readarray<int, int, false>(stream, nm, 2);
399 int N0 = nm[0], M0 = nm[1];
400 if (!((N0 >= M0 && M0 >= 0) || (N0 == -1 && M0 == -1)))
404 N = truncate ? min(N, N0) : N0;
405 M = truncate ? min(M, M0) : M0;
411 Utility::readarray<double, real, false>(stream, C);
412 if (skip) stream.seekg(std::ios::streamoff(skip), ios::cur);
413 Utility::readarray<double, real, false>(stream, S);
414 if (skip) stream.seekg(std::ios::streamoff(skip), ios::cur);
416 for (
int m = 0, k = 0; m <= M; ++m) {
417 Utility::readarray<double, real, false>(stream, &C[k], N + 1 - m);
418 stream.seekg((N0 - N) *
sizeof(
double), ios::cur);
421 if (skip) stream.seekg(std::ios::streamoff(skip), ios::cur);
422 for (
int m = 1, k = 0; m <= M; ++m) {
423 Utility::readarray<double, real, false>(stream, &S[k], N + 1 - m);
424 stream.seekg((N0 - N) *
sizeof(
double), ios::cur);
427 if (skip) stream.seekg(std::ios::streamoff(skip), ios::cur);
434 SphericalEngine::Value<true, SphericalEngine::FULL, 1>
435 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
437 SphericalEngine::Value<false, SphericalEngine::FULL, 1>
438 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
440 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
441 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
443 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
444 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
447 SphericalEngine::Value<true, SphericalEngine::FULL, 2>
448 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
450 SphericalEngine::Value<false, SphericalEngine::FULL, 2>
451 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
453 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
454 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
456 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
457 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
460 SphericalEngine::Value<true, SphericalEngine::FULL, 3>
461 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
463 SphericalEngine::Value<false, SphericalEngine::FULL, 3>
464 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
466 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
467 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
469 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
470 (
const coeff[],
const real[], real, real, real, real, real&, real&, real&);
473 SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
474 (
const coeff[],
const real[], real, real, real);
476 SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
477 (
const coeff[],
const real[], real, real, real);
479 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
480 (
const coeff[],
const real[], real, real, real);
482 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
483 (
const coeff[],
const real[], real, real, real);
486 SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
487 (
const coeff[],
const real[], real, real, real);
489 SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
490 (
const coeff[],
const real[], real, real, real);
492 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
493 (
const coeff[],
const real[], real, real, real);
495 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
496 (
const coeff[],
const real[], real, real, real);
499 SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
500 (
const coeff[],
const real[], real, real, real);
502 SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
503 (
const coeff[],
const real[], real, real, real);
505 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
506 (
const coeff[],
const real[], real, real, real);
508 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
509 (
const coeff[],
const real[], real, real, real);