137 #if defined(_MSC_VER)
140 # pragma warning (disable: 4127 4701)
148 static vector<real> sqrttable(0);
152 template<
bool gradp, SphericalEngine::normalization norm,
int L>
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)
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)
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>
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)
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)
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;
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)
388 std::vector<real>&
C,
389 std::vector<real>&
S) {
391 Utility::readarray<int, int, false>(
stream, nm, 2);
392 N = nm[0];
M = nm[1];
393 if (!(
N >=
M &&
M >= -1 &&
N *
M >= 0))
398 Utility::readarray<double, real, false>(
stream,
C);
400 Utility::readarray<double, real, false>(
stream,
S);
406 SphericalEngine::Value<true, SphericalEngine::FULL, 1>
409 SphericalEngine::Value<false, SphericalEngine::FULL, 1>
412 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
415 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
419 SphericalEngine::Value<true, SphericalEngine::FULL, 2>
422 SphericalEngine::Value<false, SphericalEngine::FULL, 2>
425 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
428 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
432 SphericalEngine::Value<true, SphericalEngine::FULL, 3>
435 SphericalEngine::Value<false, SphericalEngine::FULL, 3>
438 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
441 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
445 SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
448 SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
451 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
454 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
458 SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
461 SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
464 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
467 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
471 SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
474 SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
477 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
480 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>