GeographicLib 2.1.2
SphericalEngine.cpp
Go to the documentation of this file.
1/**
2 * \file SphericalEngine.cpp
3 * \brief Implementation for GeographicLib::SphericalEngine class
4 *
5 * Copyright (c) Charles Karney (2011-2022) <charles@karney.com> and licensed
6 * under the MIT/X11 License. For more information, see
7 * https://geographiclib.sourceforge.io/
8 *
9 * The general sum is\verbatim
10 V(r, theta, lambda) = sum(n = 0..N) sum(m = 0..n)
11 q^(n+1) * (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) * P[n,m](t)
12\endverbatim
13 * where <tt>t = cos(theta)</tt>, <tt>q = a/r</tt>. In addition write <tt>u =
14 * sin(theta)</tt>.
15 *
16 * <tt>P[n,m]</tt> is a normalized associated Legendre function of degree
17 * <tt>n</tt> and order <tt>m</tt>. Here the formulas are given for full
18 * normalized functions (usually denoted <tt>Pbar</tt>).
19 *
20 * Rewrite outer sum\verbatim
21 V(r, theta, lambda) = sum(m = 0..N) * P[m,m](t) * q^(m+1) *
22 [Sc[m] * cos(m*lambda) + Ss[m] * sin(m*lambda)]
23\endverbatim
24 * where the inner sums are\verbatim
25 Sc[m] = sum(n = m..N) q^(n-m) * C[n,m] * P[n,m](t)/P[m,m](t)
26 Ss[m] = sum(n = m..N) q^(n-m) * S[n,m] * P[n,m](t)/P[m,m](t)
27\endverbatim
28 * Evaluate sums via Clenshaw method. The overall framework is similar to
29 * Deakin with the following changes:
30 * - Clenshaw summation is used to roll the computation of
31 * <tt>cos(m*lambda)</tt> and <tt>sin(m*lambda)</tt> into the evaluation of
32 * the outer sum (rather than independently computing an array of these
33 * trigonometric terms).
34 * - Scale the coefficients to guard against overflow when <tt>N</tt> is large.
35 * .
36 * For the general framework of Clenshaw, see
37 * http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
38 *
39 * Let\verbatim
40 S = sum(k = 0..N) c[k] * F[k](x)
41 F[n+1](x) = alpha[n](x) * F[n](x) + beta[n](x) * F[n-1](x)
42\endverbatim
43 * Evaluate <tt>S</tt> with\verbatim
44 y[N+2] = y[N+1] = 0
45 y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k]
46 S = c[0] * F[0] + y[1] * F[1] + beta[1] * F[0] * y[2]
47\endverbatim
48 * \e IF <tt>F[0](x) = 1</tt> and <tt>beta(0,x) = 0</tt>, then <tt>F[1](x) =
49 * alpha(0,x)</tt> and we can continue the recursion for <tt>y[k]</tt> until
50 * <tt>y[0]</tt>, giving\verbatim
51 S = y[0]
52\endverbatim
53 *
54 * Evaluating the inner sum\verbatim
55 l = n-m; n = l+m
56 Sc[m] = sum(l = 0..N-m) C[l+m,m] * q^l * P[l+m,m](t)/P[m,m](t)
57 F[l] = q^l * P[l+m,m](t)/P[m,m](t)
58\endverbatim
59 * Holmes + Featherstone, Eq. (11), give\verbatim
60 P[n,m] = sqrt((2*n-1)*(2*n+1)/((n-m)*(n+m))) * t * P[n-1,m] -
61 sqrt((2*n+1)*(n+m-1)*(n-m-1)/((n-m)*(n+m)*(2*n-3))) * P[n-2,m]
62\endverbatim
63 * thus\verbatim
64 alpha[l] = t * q * sqrt(((2*n+1)*(2*n+3))/
65 ((n-m+1)*(n+m+1)))
66 beta[l+1] = - q^2 * sqrt(((n-m+1)*(n+m+1)*(2*n+5))/
67 ((n-m+2)*(n+m+2)*(2*n+1)))
68\endverbatim
69 * In this case, <tt>F[0] = 1</tt> and <tt>beta[0] = 0</tt>, so the <tt>Sc[m]
70 * = y[0]</tt>.
71 *
72 * Evaluating the outer sum\verbatim
73 V = sum(m = 0..N) Sc[m] * q^(m+1) * cos(m*lambda) * P[m,m](t)
74 + sum(m = 0..N) Ss[m] * q^(m+1) * cos(m*lambda) * P[m,m](t)
75 F[m] = q^(m+1) * cos(m*lambda) * P[m,m](t) [or sin(m*lambda)]
76\endverbatim
77 * Holmes + Featherstone, Eq. (13), give\verbatim
78 P[m,m] = u * sqrt((2*m+1)/((m>1?2:1)*m)) * P[m-1,m-1]
79\endverbatim
80 * also, we have\verbatim
81 cos((m+1)*lambda) = 2*cos(lambda)*cos(m*lambda) - cos((m-1)*lambda)
82\endverbatim
83 * thus\verbatim
84 alpha[m] = 2*cos(lambda) * sqrt((2*m+3)/(2*(m+1))) * u * q
85 = cos(lambda) * sqrt( 2*(2*m+3)/(m+1) ) * u * q
86 beta[m+1] = -sqrt((2*m+3)*(2*m+5)/(4*(m+1)*(m+2))) * u^2 * q^2
87 * (m == 0 ? sqrt(2) : 1)
88\endverbatim
89 * Thus\verbatim
90 F[0] = q [or 0]
91 F[1] = cos(lambda) * sqrt(3) * u * q^2 [or sin(lambda)]
92 beta[1] = - sqrt(15/4) * u^2 * q^2
93\endverbatim
94 *
95 * Here is how the various components of the gradient are computed
96 *
97 * Differentiate wrt <tt>r</tt>\verbatim
98 d q^(n+1) / dr = (-1/r) * (n+1) * q^(n+1)
99\endverbatim
100 * so multiply <tt>C[n,m]</tt> by <tt>n+1</tt> in inner sum and multiply the
101 * sum by <tt>-1/r</tt>.
102 *
103 * Differentiate wrt <tt>lambda</tt>\verbatim
104 d cos(m*lambda) = -m * sin(m*lambda)
105 d sin(m*lambda) = m * cos(m*lambda)
106\endverbatim
107 * so multiply terms by <tt>m</tt> in outer sum and swap sine and cosine
108 * variables.
109 *
110 * Differentiate wrt <tt>theta</tt>\verbatim
111 dV/dtheta = V' = -u * dV/dt = -u * V'
112\endverbatim
113 * here <tt>'</tt> denotes differentiation wrt <tt>theta</tt>.\verbatim
114 d/dtheta (Sc[m] * P[m,m](t)) = Sc'[m] * P[m,m](t) + Sc[m] * P'[m,m](t)
115\endverbatim
116 * Now <tt>P[m,m](t) = const * u^m</tt>, so <tt>P'[m,m](t) = m * t/u *
117 * P[m,m](t)</tt>, thus\verbatim
118 d/dtheta (Sc[m] * P[m,m](t)) = (Sc'[m] + m * t/u * Sc[m]) * P[m,m](t)
119\endverbatim
120 * Clenshaw recursion for <tt>Sc[m]</tt> reads\verbatim
121 y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k]
122\endverbatim
123 * Substituting <tt>alpha[k] = const * t</tt>, <tt>alpha'[k] = -u/t *
124 * alpha[k]</tt>, <tt>beta'[k] = c'[k] = 0</tt> gives\verbatim
125 y'[k] = alpha[k] * y'[k+1] + beta[k+1] * y'[k+2] - u/t * alpha[k] * y[k+1]
126\endverbatim
127 *
128 * Finally, given the derivatives of <tt>V</tt>, we can compute the components
129 * of the gradient in spherical coordinates and transform the result into
130 * cartesian coordinates.
131 **********************************************************************/
132
136
137#if defined(_MSC_VER)
138// Squelch warnings about constant conditional expressions and potentially
139// uninitialized local variables
140# pragma warning (disable: 4127 4701)
141#endif
142
143namespace GeographicLib {
144
145 using namespace std;
146
147 vector<Math::real>& SphericalEngine::sqrttable() {
148 static vector<real> sqrttable(0);
149 return sqrttable;
150 }
151
152 template<bool gradp, SphericalEngine::normalization norm, int L>
153 Math::real SphericalEngine::Value(const coeff c[], const real f[],
154 real x, real y, real z, real a,
155 real& gradx, real& grady, real& gradz)
156 {
157 static_assert(L > 0, "L must be positive");
158 static_assert(norm == FULL || norm == SCHMIDT, "Unknown normalization");
159 int N = c[0].nmx(), M = c[0].mmx();
160
161 real
162 p = hypot(x, y),
163 cl = p != 0 ? x / p : 1, // cos(lambda); at pole, pick lambda = 0
164 sl = p != 0 ? y / p : 0, // sin(lambda)
165 r = hypot(z, p),
166 t = r != 0 ? z / r : 0, // cos(theta); at origin, pick theta = pi/2
167 u = r != 0 ? fmax(p / r, eps()) : 1, // sin(theta); but avoid the pole
168 q = a / r;
169 real
170 q2 = Math::sq(q),
171 uq = u * q,
172 uq2 = Math::sq(uq),
173 tu = t / u;
174 // Initialize outer sum
175 real vc = 0, vc2 = 0, vs = 0, vs2 = 0; // v [N + 1], v [N + 2]
176 // vr, vt, vl and similar w variable accumulate the sums for the
177 // derivatives wrt r, theta, and lambda, respectively.
178 real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0; // vr[N + 1], vr[N + 2]
179 real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0; // vt[N + 1], vt[N + 2]
180 real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0; // vl[N + 1], vl[N + 2]
181 int k[L];
182 const vector<real>& root( sqrttable() );
183 for (int m = M; m >= 0; --m) { // m = M .. 0
184 // Initialize inner sum
185 real
186 wc = 0, wc2 = 0, ws = 0, ws2 = 0, // w [N - m + 1], w [N - m + 2]
187 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0, // wr[N - m + 1], wr[N - m + 2]
188 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0; // wt[N - m + 1], wt[N - m + 2]
189 for (int l = 0; l < L; ++l)
190 k[l] = c[l].index(N, m) + 1;
191 for (int n = N; n >= m; --n) { // n = N .. m; l = N - m .. 0
192 real w, A, Ax, B, R; // alpha[l], beta[l + 1]
193 switch (norm) {
194 case FULL:
195 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);
196 Ax = q * w * root[2 * n + 3];
197 A = t * Ax;
198 B = - q2 * root[2 * n + 5] /
199 (w * root[n - m + 2] * root[n + m + 2]);
200 break;
201 case SCHMIDT:
202 w = root[n - m + 1] * root[n + m + 1];
203 Ax = q * (2 * n + 1) / w;
204 A = t * Ax;
205 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);
206 break;
207 default: break; // To suppress warning message from Visual Studio
208 }
209 R = c[0].Cv(--k[0]);
210 for (int l = 1; l < L; ++l)
211 R += c[l].Cv(--k[l], n, m, f[l]);
212 R *= scale();
213 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
214 if (gradp) {
215 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
216 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
217 }
218 if (m) {
219 R = c[0].Sv(k[0]);
220 for (int l = 1; l < L; ++l)
221 R += c[l].Sv(k[l], n, m, f[l]);
222 R *= scale();
223 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
224 if (gradp) {
225 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
226 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
227 }
228 }
229 }
230 // Now Sc[m] = wc, Ss[m] = ws
231 // Sc'[m] = wtc, Ss'[m] = wtc
232 if (m) {
233 real v, A, B; // alpha[m], beta[m + 1]
234 switch (norm) {
235 case FULL:
236 v = root[2] * root[2 * m + 3] / root[m + 1];
237 A = cl * v * uq;
238 B = - v * root[2 * m + 5] / (root[8] * root[m + 2]) * uq2;
239 break;
240 case SCHMIDT:
241 v = root[2] * root[2 * m + 1] / root[m + 1];
242 A = cl * v * uq;
243 B = - v * root[2 * m + 3] / (root[8] * root[m + 2]) * uq2;
244 break;
245 default: break; // To suppress warning message from Visual Studio
246 }
247 v = A * vc + B * vc2 + wc ; vc2 = vc ; vc = v;
248 v = A * vs + B * vs2 + ws ; vs2 = vs ; vs = v;
249 if (gradp) {
250 // Include the terms Sc[m] * P'[m,m](t) and Ss[m] * P'[m,m](t)
251 wtc += m * tu * wc; wts += m * tu * ws;
252 v = A * vrc + B * vrc2 + wrc; vrc2 = vrc; vrc = v;
253 v = A * vrs + B * vrs2 + wrs; vrs2 = vrs; vrs = v;
254 v = A * vtc + B * vtc2 + wtc; vtc2 = vtc; vtc = v;
255 v = A * vts + B * vts2 + wts; vts2 = vts; vts = v;
256 v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;
257 v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;
258 }
259 } else {
260 real A, B, qs;
261 switch (norm) {
262 case FULL:
263 A = root[3] * uq; // F[1]/(q*cl) or F[1]/(q*sl)
264 B = - root[15]/2 * uq2; // beta[1]/q
265 break;
266 case SCHMIDT:
267 A = uq;
268 B = - root[3]/2 * uq2;
269 break;
270 default: break; // To suppress warning message from Visual Studio
271 }
272 qs = q / scale();
273 vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);
274 if (gradp) {
275 qs /= r;
276 // The components of the gradient in spherical coordinates are
277 // r: dV/dr
278 // theta: 1/r * dV/dtheta
279 // lambda: 1/(r*u) * dV/dlambda
280 vrc = - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);
281 vtc = qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);
282 vlc = qs / u * ( A * (cl * vlc + sl * vls) + B * vlc2);
283 }
284 }
285 }
286
287 if (gradp) {
288 // Rotate into cartesian (geocentric) coordinates
289 gradx = cl * (u * vrc + t * vtc) - sl * vlc;
290 grady = sl * (u * vrc + t * vtc) + cl * vlc;
291 gradz = t * vrc - u * vtc ;
292 }
293 return vc;
294 }
295
296 template<bool gradp, SphericalEngine::normalization norm, int L>
298 real p, real z, real a) {
299
300 static_assert(L > 0, "L must be positive");
301 static_assert(norm == FULL || norm == SCHMIDT, "Unknown normalization");
302 int N = c[0].nmx(), M = c[0].mmx();
303
304 real
305 r = hypot(z, p),
306 t = r != 0 ? z / r : 0, // cos(theta); at origin, pick theta = pi/2
307 u = r != 0 ? fmax(p / r, eps()) : 1, // sin(theta); but avoid the pole
308 q = a / r;
309 real
310 q2 = Math::sq(q),
311 tu = t / u;
312 CircularEngine circ(M, gradp, norm, a, r, u, t);
313 int k[L];
314 const vector<real>& root( sqrttable() );
315 for (int m = M; m >= 0; --m) { // m = M .. 0
316 // Initialize inner sum
317 real
318 wc = 0, wc2 = 0, ws = 0, ws2 = 0, // w [N - m + 1], w [N - m + 2]
319 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0, // wr[N - m + 1], wr[N - m + 2]
320 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0; // wt[N - m + 1], wt[N - m + 2]
321 for (int l = 0; l < L; ++l)
322 k[l] = c[l].index(N, m) + 1;
323 for (int n = N; n >= m; --n) { // n = N .. m; l = N - m .. 0
324 real w, A, Ax, B, R; // alpha[l], beta[l + 1]
325 switch (norm) {
326 case FULL:
327 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);
328 Ax = q * w * root[2 * n + 3];
329 A = t * Ax;
330 B = - q2 * root[2 * n + 5] /
331 (w * root[n - m + 2] * root[n + m + 2]);
332 break;
333 case SCHMIDT:
334 w = root[n - m + 1] * root[n + m + 1];
335 Ax = q * (2 * n + 1) / w;
336 A = t * Ax;
337 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);
338 break;
339 default: break; // To suppress warning message from Visual Studio
340 }
341 R = c[0].Cv(--k[0]);
342 for (int l = 1; l < L; ++l)
343 R += c[l].Cv(--k[l], n, m, f[l]);
344 R *= scale();
345 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
346 if (gradp) {
347 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
348 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
349 }
350 if (m) {
351 R = c[0].Sv(k[0]);
352 for (int l = 1; l < L; ++l)
353 R += c[l].Sv(k[l], n, m, f[l]);
354 R *= scale();
355 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
356 if (gradp) {
357 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
358 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
359 }
360 }
361 }
362 if (!gradp)
363 circ.SetCoeff(m, wc, ws);
364 else {
365 // Include the terms Sc[m] * P'[m,m](t) and Ss[m] * P'[m,m](t)
366 wtc += m * tu * wc; wts += m * tu * ws;
367 circ.SetCoeff(m, wc, ws, wrc, wrs, wtc, wts);
368 }
369 }
370
371 return circ;
372 }
373
375 // Need square roots up to max(2 * N + 5, 15).
376 vector<real>& root( sqrttable() );
377 int L = max(2 * N + 5, 15) + 1, oldL = int(root.size());
378 if (oldL >= L)
379 return;
380 root.resize(L);
381 for (int l = oldL; l < L; ++l)
382 root[l] = sqrt(real(l));
383 }
384
385 void SphericalEngine::coeff::readcoeffs(istream& stream, int& N, int& M,
386 vector<real>& C,
387 vector<real>& S,
388 bool truncate) {
389 if (truncate) {
390 if (!((N >= M && M >= 0) || (N == -1 && M == -1)))
391 // The last condition is that M = -1 implies N = -1.
392 throw GeographicErr("Bad requested degree and order " +
393 Utility::str(N) + " " + Utility::str(M));
394 }
395 int nm[2];
396 Utility::readarray<int, int, false>(stream, nm, 2);
397 int N0 = nm[0], M0 = nm[1];
398 if (!((N0 >= M0 && M0 >= 0) || (N0 == -1 && M0 == -1)))
399 // The last condition is that M0 = -1 implies N0 = -1.
400 throw GeographicErr("Bad degree and order " +
401 Utility::str(N0) + " " + Utility::str(M0));
402 N = truncate ? min(N, N0) : N0;
403 M = truncate ? min(M, M0) : M0;
404 C.resize(SphericalEngine::coeff::Csize(N, M));
405 S.resize(SphericalEngine::coeff::Ssize(N, M));
406 int skip = (SphericalEngine::coeff::Csize(N0, M0) -
407 SphericalEngine::coeff::Csize(N0, M )) * sizeof(double);
408 if (N == N0) {
409 Utility::readarray<double, real, false>(stream, C);
410 if (skip) stream.seekg(streamoff(skip), ios::cur);
411 Utility::readarray<double, real, false>(stream, S);
412 if (skip) stream.seekg(streamoff(skip), ios::cur);
413 } else {
414 for (int m = 0, k = 0; m <= M; ++m) {
415 Utility::readarray<double, real, false>(stream, &C[k], N + 1 - m);
416 stream.seekg((N0 - N) * sizeof(double), ios::cur);
417 k += N + 1 - m;
418 }
419 if (skip) stream.seekg(streamoff(skip), ios::cur);
420 for (int m = 1, k = 0; m <= M; ++m) {
421 Utility::readarray<double, real, false>(stream, &S[k], N + 1 - m);
422 stream.seekg((N0 - N) * sizeof(double), ios::cur);
423 k += N + 1 - m;
424 }
425 if (skip) stream.seekg(streamoff(skip), ios::cur);
426 }
427 return;
428 }
429
430 /// \cond SKIP
432 SphericalEngine::Value<true, SphericalEngine::FULL, 1>
433 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
435 SphericalEngine::Value<false, SphericalEngine::FULL, 1>
436 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
438 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
439 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
441 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
442 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
443
445 SphericalEngine::Value<true, SphericalEngine::FULL, 2>
446 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
448 SphericalEngine::Value<false, SphericalEngine::FULL, 2>
449 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
451 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
452 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
454 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
455 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
456
458 SphericalEngine::Value<true, SphericalEngine::FULL, 3>
459 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
461 SphericalEngine::Value<false, SphericalEngine::FULL, 3>
462 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
464 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
465 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
467 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
468 (const coeff[], const real[], real, real, real, real, real&, real&, real&);
469
471 SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
472 (const coeff[], const real[], real, real, real);
474 SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
475 (const coeff[], const real[], real, real, real);
477 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
478 (const coeff[], const real[], real, real, real);
480 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
481 (const coeff[], const real[], real, real, real);
482
484 SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
485 (const coeff[], const real[], real, real, real);
487 SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
488 (const coeff[], const real[], real, real, real);
490 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
491 (const coeff[], const real[], real, real, real);
493 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
494 (const coeff[], const real[], real, real, real);
495
497 SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
498 (const coeff[], const real[], real, real, real);
500 SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
501 (const coeff[], const real[], real, real, real);
503 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
504 (const coeff[], const real[], real, real, real);
506 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
507 (const coeff[], const real[], real, real, real);
508 /// \endcond
509
510} // namespace GeographicLib
Header for GeographicLib::CircularEngine class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
Header for GeographicLib::SphericalEngine class.
Header for GeographicLib::Utility class.
Spherical harmonic sums for a circle.
Exception handling for GeographicLib.
Definition: Constants.hpp:316
static T sq(T x)
Definition: Math.hpp:212
Package up coefficients for SphericalEngine.
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)
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 CircularEngine Circle(const coeff c[], const real f[], real p, real z, real a)
static std::string str(T x, int p=-1)
Definition: Utility.hpp:161
Namespace for GeographicLib.
Definition: Accumulator.cpp:12