GeographicLib  1.35
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) <charles@karney.com> and licensed under
6  * the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
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 
135 #include <GeographicLib/Utility.hpp>
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 
143 namespace GeographicLib {
144 
145  using namespace std;
146 
147  const Math::real SphericalEngine::scale_ =
148  pow(real(numeric_limits<real>::radix),
149  -3 * numeric_limits<real>::max_exponent / 5);
150  const Math::real SphericalEngine::eps_ =
151  numeric_limits<real>::epsilon() * sqrt(numeric_limits<real>::epsilon());
152 
153  const vector<Math::real> SphericalEngine::Z_(0);
154  vector<Math::real> SphericalEngine::root_(0);
155 
156  template<bool gradp, SphericalEngine::normalization norm, int L>
157  Math::real SphericalEngine::Value(const coeff c[], const real f[],
158  real x, real y, real z, real a,
159  real& gradx, real& grady, real& gradz)
160  throw() {
161  STATIC_ASSERT(L > 0, "L must be positive");
162  STATIC_ASSERT(norm == FULL || norm == SCHMIDT, "Unknown normalization");
163  int N = c[0].nmx(), M = c[0].mmx();
164 
165  real
166  p = Math::hypot(x, y),
167  cl = p ? x / p : 1, // cos(lambda); at pole, pick lambda = 0
168  sl = p ? y / p : 0, // sin(lambda)
169  r = Math::hypot(z, p),
170  t = r ? z / r : 0, // cos(theta); at origin, pick theta = pi/2
171  u = r ? max(p / r, eps_) : 1, // sin(theta); but avoid the pole
172  q = a / r;
173  real
174  q2 = Math::sq(q),
175  uq = u * q,
176  uq2 = Math::sq(uq),
177  tu = t / u;
178  // Initialize outer sum
179  real vc = 0, vc2 = 0, vs = 0, vs2 = 0; // v [N + 1], v [N + 2]
180  // vr, vt, vl and similar w variable accumulate the sums for the
181  // derivatives wrt r, theta, and lambda, respectively.
182  real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0; // vr[N + 1], vr[N + 2]
183  real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0; // vt[N + 1], vt[N + 2]
184  real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0; // vl[N + 1], vl[N + 2]
185  int k[L];
186  for (int m = M; m >= 0; --m) { // m = M .. 0
187  // Initialize inner sum
188  real wc = 0, wc2 = 0, ws = 0, ws2 = 0; // w [N - m + 1], w [N - m + 2]
189  real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0; // wr[N - m + 1], wr[N - m + 2]
190  real wtc = 0, wtc2 = 0, wts = 0, wts2 = 0; // wt[N - m + 1], wt[N - m + 2]
191  for (int l = 0; l < L; ++l)
192  k[l] = c[l].index(N, m) + 1;
193  for (int n = N; n >= m; --n) { // n = N .. m; l = N - m .. 0
194  real w, A, Ax, B, R; // alpha[l], beta[l + 1]
195  switch (norm) {
196  case FULL:
197  w = root_[2 * n + 1] / (root_[n - m + 1] * root_[n + m + 1]);
198  Ax = q * w * root_[2 * n + 3];
199  A = t * Ax;
200  B = - q2 * root_[2 * n + 5] /
201  (w * root_[n - m + 2] * root_[n + m + 2]);
202  break;
203  case SCHMIDT:
204  w = root_[n - m + 1] * root_[n + m + 1];
205  Ax = q * (2 * n + 1) / w;
206  A = t * Ax;
207  B = - q2 * w / (root_[n - m + 2] * root_[n + m + 2]);
208  break;
209  default: break; // To suppress warning message from Visual Studio
210  }
211  R = c[0].Cv(--k[0]);
212  for (int l = 1; l < L; ++l)
213  R += c[l].Cv(--k[l], n, m, f[l]);
214  R *= scale_;
215  w = A * wc + B * wc2 + R; wc2 = wc; wc = w;
216  if (gradp) {
217  w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;
218  w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;
219  }
220  if (m) {
221  R = c[0].Sv(k[0]);
222  for (int l = 1; l < L; ++l)
223  R += c[l].Sv(k[l], n, m, f[l]);
224  R *= scale_;
225  w = A * ws + B * ws2 + R; ws2 = ws; ws = w;
226  if (gradp) {
227  w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;
228  w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;
229  }
230  }
231  }
232  // Now Sc[m] = wc, Ss[m] = ws
233  // Sc'[m] = wtc, Ss'[m] = wtc
234  if (m) {
235  real v, A, B; // alpha[m], beta[m + 1]
236  switch (norm) {
237  case FULL:
238  v = root_[2] * root_[2 * m + 3] / root_[m + 1];
239  A = cl * v * uq;
240  B = - v * root_[2 * m + 5] / (root_[8] * root_[m + 2]) * uq2;
241  break;
242  case SCHMIDT:
243  v = root_[2] * root_[2 * m + 1] / root_[m + 1];
244  A = cl * v * uq;
245  B = - v * root_[2 * m + 3] / (root_[8] * root_[m + 2]) * uq2;
246  break;
247  default: break; // To suppress warning message from Visual Studio
248  }
249  v = A * vc + B * vc2 + wc ; vc2 = vc ; vc = v;
250  v = A * vs + B * vs2 + ws ; vs2 = vs ; vs = v;
251  if (gradp) {
252  // Include the terms Sc[m] * P'[m,m](t) and Ss[m] * P'[m,m](t)
253  wtc += m * tu * wc; wts += m * tu * ws;
254  v = A * vrc + B * vrc2 + wrc; vrc2 = vrc; vrc = v;
255  v = A * vrs + B * vrs2 + wrs; vrs2 = vrs; vrs = v;
256  v = A * vtc + B * vtc2 + wtc; vtc2 = vtc; vtc = v;
257  v = A * vts + B * vts2 + wts; vts2 = vts; vts = v;
258  v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;
259  v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;
260  }
261  } else {
262  real A, B, qs;
263  switch (norm) {
264  case FULL:
265  A = root_[3] * uq; // F[1]/(q*cl) or F[1]/(q*sl)
266  B = - root_[15]/2 * uq2; // beta[1]/q
267  break;
268  case SCHMIDT:
269  A = uq;
270  B = - root_[3]/2 * uq2;
271  break;
272  default: break; // To suppress warning message from Visual Studio
273  }
274  qs = q / scale_;
275  vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);
276  if (gradp) {
277  qs /= r;
278  // The components of the gradient in spherical coordinates are
279  // r: dV/dr
280  // theta: 1/r * dV/dtheta
281  // lambda: 1/(r*u) * dV/dlambda
282  vrc = - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);
283  vtc = qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);
284  vlc = qs / u * ( A * (cl * vlc + sl * vls) + B * vlc2);
285  }
286  }
287  }
288 
289  if (gradp) {
290  // Rotate into cartesian (geocentric) coordinates
291  gradx = cl * (u * vrc + t * vtc) - sl * vlc;
292  grady = sl * (u * vrc + t * vtc) + cl * vlc;
293  gradz = t * vrc - u * vtc ;
294  }
295  return vc;
296  }
297 
298  template<bool gradp, SphericalEngine::normalization norm, int L>
299  CircularEngine SphericalEngine::Circle(const coeff c[], const real f[],
300  real p, real z, real a) {
301 
302  STATIC_ASSERT(L > 0, "L must be positive");
303  STATIC_ASSERT(norm == FULL || norm == SCHMIDT, "Unknown normalization");
304  int N = c[0].nmx(), M = c[0].mmx();
305 
306  real
307  r = Math::hypot(z, p),
308  t = r ? z / r : 0, // cos(theta); at origin, pick theta = pi/2
309  u = r ? max(p / r, eps_) : 1, // sin(theta); but avoid the pole
310  q = a / r;
311  real
312  q2 = Math::sq(q),
313  tu = t / u;
314  CircularEngine circ(M, gradp, norm, a, r, u, t);
315  int k[L];
316  for (int m = M; m >= 0; --m) { // m = M .. 0
317  // Initialize inner sum
318  real wc = 0, wc2 = 0, ws = 0, ws2 = 0; // w [N - m + 1], w [N - m + 2]
319  real wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0; // wr[N - m + 1], wr[N - m + 2]
320  real 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  int L = max(2 * N + 5, 15) + 1, oldL = int(root_.size());
377  if (oldL >= L)
378  return;
379  root_.resize(L);
380  for (int l = oldL; l < L; ++l)
381  root_[l] = sqrt(real(l));
382  }
383 
384  void SphericalEngine::coeff::readcoeffs(std::istream& stream, int& N, int& M,
385  std::vector<real>& C,
386  std::vector<real>& S) {
387  int nm[2];
388  Utility::readarray<int, int, false>(stream, nm, 2);
389  N = nm[0]; M = nm[1];
390  if (!(N >= M && M >= -1 && N * M >= 0))
391  // The last condition is that M = -1 implies N = -1 and vice versa.
392  throw GeographicErr("Bad degree and order " +
393  Utility::str(N) + " " + Utility::str(M));
394  C.resize(SphericalEngine::coeff::Csize(N, M));
395  Utility::readarray<double, real, false>(stream, C);
396  S.resize(SphericalEngine::coeff::Ssize(N, M));
397  Utility::readarray<double, real, false>(stream, S);
398  return;
399  }
400 
401  /// \cond SKIP
403  SphericalEngine::Value<true, SphericalEngine::FULL, 1>
404  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
406  SphericalEngine::Value<false, SphericalEngine::FULL, 1>
407  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
409  SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
410  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
412  SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
413  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
414 
416  SphericalEngine::Value<true, SphericalEngine::FULL, 2>
417  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
419  SphericalEngine::Value<false, SphericalEngine::FULL, 2>
420  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
422  SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
423  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
425  SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
426  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
427 
429  SphericalEngine::Value<true, SphericalEngine::FULL, 3>
430  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
432  SphericalEngine::Value<false, SphericalEngine::FULL, 3>
433  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
435  SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
436  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
438  SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
439  (const coeff[], const real[], real, real, real, real, real&, real&, real&);
440 
442  SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
443  (const coeff[], const real[], real, real, real);
445  SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
446  (const coeff[], const real[], real, real, real);
448  SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
449  (const coeff[], const real[], real, real, real);
451  SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
452  (const coeff[], const real[], real, real, real);
453 
455  SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
456  (const coeff[], const real[], real, real, real);
458  SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
459  (const coeff[], const real[], real, real, real);
461  SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
462  (const coeff[], const real[], real, real, real);
464  SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
465  (const coeff[], const real[], real, real, real);
466 
468  SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
469  (const coeff[], const real[], real, real, real);
471  SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
472  (const coeff[], const real[], real, real, real);
474  SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
475  (const coeff[], const real[], real, real, real);
477  SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
478  (const coeff[], const real[], real, real, real);
479  /// \endcond
480 
481 } // namespace GeographicLib
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
Header for GeographicLib::Utility class.
static CircularEngine Circle(const coeff c[], const real f[], real p, real z, real a)
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S)
static T hypot(T x, T y)
Definition: Math.hpp:165
Package up coefficients for SphericalEngine.
static T sq(T x)
Definition: Math.hpp:153
static std::string str(T x, int p=-1)
Definition: Utility.hpp:266
#define STATIC_ASSERT(cond, reason)
Definition: Constants.hpp:37
Header for GeographicLib::CircularEngine class.
Spherical harmonic sums for a circle.
Exception handling for GeographicLib.
Definition: Constants.hpp:320
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.