GeographicLib  1.45
SphericalHarmonic.hpp
Go to the documentation of this file.
1 /**
2  * \file SphericalHarmonic.hpp
3  * \brief Header for GeographicLib::SphericalHarmonic 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 
10 #if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC_HPP)
11 #define GEOGRAPHICLIB_SPHERICALHARMONIC_HPP 1
12 
13 #include <vector>
17 
18 namespace GeographicLib {
19 
20  /**
21  * \brief Spherical harmonic series
22  *
23  * This class evaluates the spherical harmonic sum \verbatim
24  V(x, y, z) = sum(n = 0..N)[ q^(n+1) * sum(m = 0..n)[
25  (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) *
26  P[n,m](cos(theta)) ] ]
27  \endverbatim
28  * where
29  * - <i>p</i><sup>2</sup> = <i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>,
30  * - <i>r</i><sup>2</sup> = <i>p</i><sup>2</sup> + <i>z</i><sup>2</sup>,
31  * - \e q = <i>a</i>/<i>r</i>,
32  * - &theta; = atan2(\e p, \e z) = the spherical \e colatitude,
33  * - &lambda; = atan2(\e y, \e x) = the longitude.
34  * - P<sub><i>nm</i></sub>(\e t) is the associated Legendre polynomial of
35  * degree \e n and order \e m.
36  *
37  * Two normalizations are supported for P<sub><i>nm</i></sub>
38  * - fully normalized denoted by SphericalHarmonic::FULL.
39  * - Schmidt semi-normalized denoted by SphericalHarmonic::SCHMIDT.
40  *
41  * Clenshaw summation is used for the sums over both \e n and \e m. This
42  * allows the computation to be carried out without the need for any
43  * temporary arrays. See SphericalEngine.cpp for more information on the
44  * implementation.
45  *
46  * References:
47  * - C. W. Clenshaw,
48  * <a href="https://dx.doi.org/10.1090/S0025-5718-1955-0071856-0">
49  * A note on the summation of Chebyshev series</a>,
50  * %Math. Tables Aids Comput. 9(51), 118--120 (1955).
51  * - R. E. Deakin, Derivatives of the earth's potentials, Geomatics
52  * Research Australasia 68, 31--60, (June 1998).
53  * - W. A. Heiskanen and H. Moritz, Physical Geodesy, (Freeman, San
54  * Francisco, 1967). (See Sec. 1-14, for a definition of Pbar.)
55  * - S. A. Holmes and W. E. Featherstone,
56  * <a href="https://dx.doi.org/10.1007/s00190-002-0216-2">
57  * A unified approach to the Clenshaw summation and the recursive
58  * computation of very high degree and order normalised associated Legendre
59  * functions</a>, J. Geodesy 76(5), 279--299 (2002).
60  * - C. C. Tscherning and K. Poder,
61  * <a href="http://cct.gfy.ku.dk/publ_cct/cct80.pdf">
62  * Some geodetic applications of Clenshaw summation</a>,
63  * Boll. Geod. Sci. Aff. 41(4), 349--375 (1982).
64  *
65  * Example of use:
66  * \include example-SphericalHarmonic.cpp
67  **********************************************************************/
68 
70  public:
71  /**
72  * Supported normalizations for the associated Legendre polynomials.
73  **********************************************************************/
75  /**
76  * Fully normalized associated Legendre polynomials.
77  *
78  * These are defined by
79  * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(\e z)
80  * = (&minus;1)<sup><i>m</i></sup>
81  * sqrt(\e k (2\e n + 1) (\e n &minus; \e m)! / (\e n + \e m)!)
82  * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where
83  * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers
84  * function (also known as the Legendre function on the cut or the
85  * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k
86  * = 1 for \e m = 0 and \e k = 2 otherwise.
87  *
88  * The mean squared value of
89  * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos&theta;)
90  * cos(<i>m</i>&lambda;) and
91  * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos&theta;)
92  * sin(<i>m</i>&lambda;) over the sphere is 1.
93  *
94  * @hideinitializer
95  **********************************************************************/
97  /**
98  * Schmidt semi-normalized associated Legendre polynomials.
99  *
100  * These are defined by
101  * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(\e z)
102  * = (&minus;1)<sup><i>m</i></sup>
103  * sqrt(\e k (\e n &minus; \e m)! / (\e n + \e m)!)
104  * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where
105  * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers
106  * function (also known as the Legendre function on the cut or the
107  * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k
108  * = 1 for \e m = 0 and \e k = 2 otherwise.
109  *
110  * The mean squared value of
111  * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos&theta;)
112  * cos(<i>m</i>&lambda;) and
113  * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos&theta;)
114  * sin(<i>m</i>&lambda;) over the sphere is 1/(2\e n + 1).
115  *
116  * @hideinitializer
117  **********************************************************************/
119  /// \cond SKIP
120  // These are deprecated...
121  full = FULL,
122  schmidt = SCHMIDT,
123  /// \endcond
124  };
125 
126  private:
127  typedef Math::real real;
129  real _a;
130  unsigned _norm;
131 
132  public:
133  /**
134  * Constructor with a full set of coefficients specified.
135  *
136  * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>.
137  * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>.
138  * @param[in] N the maximum degree and order of the sum
139  * @param[in] a the reference radius appearing in the definition of the
140  * sum.
141  * @param[in] norm the normalization for the associated Legendre
142  * polynomials, either SphericalHarmonic::FULL (the default) or
143  * SphericalHarmonic::SCHMIDT.
144  * @exception GeographicErr if \e N does not satisfy \e N &ge; &minus;1.
145  * @exception GeographicErr if \e C or \e S is not big enough to hold the
146  * coefficients.
147  *
148  * The coefficients <i>C</i><sub><i>nm</i></sub> and
149  * <i>S</i><sub><i>nm</i></sub> are stored in the one-dimensional vectors
150  * \e C and \e S which must contain (\e N + 1)(\e N + 2)/2 and \e N (\e N +
151  * 1)/2 elements, respectively, stored in "column-major" order. Thus for
152  * \e N = 3, the order would be:
153  * <i>C</i><sub>00</sub>,
154  * <i>C</i><sub>10</sub>,
155  * <i>C</i><sub>20</sub>,
156  * <i>C</i><sub>30</sub>,
157  * <i>C</i><sub>11</sub>,
158  * <i>C</i><sub>21</sub>,
159  * <i>C</i><sub>31</sub>,
160  * <i>C</i><sub>22</sub>,
161  * <i>C</i><sub>32</sub>,
162  * <i>C</i><sub>33</sub>.
163  * In general the (\e n,\e m) element is at index \e m \e N &minus; \e m
164  * (\e m &minus; 1)/2 + \e n. The layout of \e S is the same except that
165  * the first column is omitted (since the \e m = 0 terms never contribute
166  * to the sum) and the 0th element is <i>S</i><sub>11</sub>
167  *
168  * The class stores <i>pointers</i> to the first elements of \e C and \e S.
169  * These arrays should not be altered or destroyed during the lifetime of a
170  * SphericalHarmonic object.
171  **********************************************************************/
172  SphericalHarmonic(const std::vector<real>& C,
173  const std::vector<real>& S,
174  int N, real a, unsigned norm = FULL)
175  : _a(a)
176  , _norm(norm)
177  { _c[0] = SphericalEngine::coeff(C, S, N); }
178 
179  /**
180  * Constructor with a subset of coefficients specified.
181  *
182  * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>.
183  * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>.
184  * @param[in] N the degree used to determine the layout of \e C and \e S.
185  * @param[in] nmx the maximum degree used in the sum. The sum over \e n is
186  * from 0 thru \e nmx.
187  * @param[in] mmx the maximum order used in the sum. The sum over \e m is
188  * from 0 thru min(\e n, \e mmx).
189  * @param[in] a the reference radius appearing in the definition of the
190  * sum.
191  * @param[in] norm the normalization for the associated Legendre
192  * polynomials, either SphericalHarmonic::FULL (the default) or
193  * SphericalHarmonic::SCHMIDT.
194  * @exception GeographicErr if \e N, \e nmx, and \e mmx do not satisfy
195  * \e N &ge; \e nmx &ge; \e mmx &ge; &minus;1.
196  * @exception GeographicErr if \e C or \e S is not big enough to hold the
197  * coefficients.
198  *
199  * The class stores <i>pointers</i> to the first elements of \e C and \e S.
200  * These arrays should not be altered or destroyed during the lifetime of a
201  * SphericalHarmonic object.
202  **********************************************************************/
203  SphericalHarmonic(const std::vector<real>& C,
204  const std::vector<real>& S,
205  int N, int nmx, int mmx,
206  real a, unsigned norm = FULL)
207  : _a(a)
208  , _norm(norm)
209  { _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); }
210 
211  /**
212  * A default constructor so that the object can be created when the
213  * constructor for another object is initialized. This default object can
214  * then be reset with the default copy assignment operator.
215  **********************************************************************/
217 
218  /**
219  * Compute the spherical harmonic sum.
220  *
221  * @param[in] x cartesian coordinate.
222  * @param[in] y cartesian coordinate.
223  * @param[in] z cartesian coordinate.
224  * @return \e V the spherical harmonic sum.
225  *
226  * This routine requires constant memory and thus never throws an
227  * exception.
228  **********************************************************************/
229  Math::real operator()(real x, real y, real z) const {
230  real f[] = {1};
231  real v = 0;
232  real dummy;
233  switch (_norm) {
234  case FULL:
235  v = SphericalEngine::Value<false, SphericalEngine::FULL, 1>
236  (_c, f, x, y, z, _a, dummy, dummy, dummy);
237  break;
238  case SCHMIDT:
239  v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
240  (_c, f, x, y, z, _a, dummy, dummy, dummy);
241  break;
242  }
243  return v;
244  }
245 
246  /**
247  * Compute a spherical harmonic sum and its gradient.
248  *
249  * @param[in] x cartesian coordinate.
250  * @param[in] y cartesian coordinate.
251  * @param[in] z cartesian coordinate.
252  * @param[out] gradx \e x component of the gradient
253  * @param[out] grady \e y component of the gradient
254  * @param[out] gradz \e z component of the gradient
255  * @return \e V the spherical harmonic sum.
256  *
257  * This is the same as the previous function, except that the components of
258  * the gradients of the sum in the \e x, \e y, and \e z directions are
259  * computed. This routine requires constant memory and thus never throws
260  * an exception.
261  **********************************************************************/
262  Math::real operator()(real x, real y, real z,
263  real& gradx, real& grady, real& gradz) const {
264  real f[] = {1};
265  real v = 0;
266  switch (_norm) {
267  case FULL:
268  v = SphericalEngine::Value<true, SphericalEngine::FULL, 1>
269  (_c, f, x, y, z, _a, gradx, grady, gradz);
270  break;
271  case SCHMIDT:
272  v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
273  (_c, f, x, y, z, _a, gradx, grady, gradz);
274  break;
275  }
276  return v;
277  }
278 
279  /**
280  * Create a CircularEngine to allow the efficient evaluation of several
281  * points on a circle of latitude.
282  *
283  * @param[in] p the radius of the circle.
284  * @param[in] z the height of the circle above the equatorial plane.
285  * @param[in] gradp if true the returned object will be able to compute the
286  * gradient of the sum.
287  * @exception std::bad_alloc if the memory for the CircularEngine can't be
288  * allocated.
289  * @return the CircularEngine object.
290  *
291  * SphericalHarmonic::operator()() exchanges the order of the sums in the
292  * definition, i.e., &sum;<sub><i>n</i> = 0..<i>N</i></sub>
293  * &sum;<sub><i>m</i> = 0..<i>n</i></sub> becomes &sum;<sub><i>m</i> =
294  * 0..<i>N</i></sub> &sum;<sub><i>n</i> = <i>m</i>..<i>N</i></sub>.
295  * SphericalHarmonic::Circle performs the inner sum over degree \e n (which
296  * entails about <i>N</i><sup>2</sup> operations). Calling
297  * CircularEngine::operator()() on the returned object performs the outer
298  * sum over the order \e m (about \e N operations).
299  *
300  * Here's an example of computing the spherical sum at a sequence of
301  * longitudes without using a CircularEngine object \code
302  SphericalHarmonic h(...); // Create the SphericalHarmonic object
303  double r = 2, lat = 33, lon0 = 44, dlon = 0.01;
304  double
305  phi = lat * Math::degree<double>(),
306  z = r * sin(phi), p = r * cos(phi);
307  for (int i = 0; i <= 100; ++i) {
308  real
309  lon = lon0 + i * dlon,
310  lam = lon * Math::degree<double>();
311  std::cout << lon << " " << h(p * cos(lam), p * sin(lam), z) << "\n";
312  }
313  \endcode
314  * Here is the same calculation done using a CircularEngine object. This
315  * will be about <i>N</i>/2 times faster. \code
316  SphericalHarmonic h(...); // Create the SphericalHarmonic object
317  double r = 2, lat = 33, lon0 = 44, dlon = 0.01;
318  double
319  phi = lat * Math::degree<double>(),
320  z = r * sin(phi), p = r * cos(phi);
321  CircularEngine c(h(p, z, false)); // Create the CircularEngine object
322  for (int i = 0; i <= 100; ++i) {
323  real
324  lon = lon0 + i * dlon;
325  std::cout << lon << " " << c(lon) << "\n";
326  }
327  \endcode
328  **********************************************************************/
329  CircularEngine Circle(real p, real z, bool gradp) const {
330  real f[] = {1};
331  switch (_norm) {
332  case FULL:
333  return gradp ?
334  SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
335  (_c, f, p, z, _a) :
336  SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
337  (_c, f, p, z, _a);
338  break;
339  case SCHMIDT:
340  default: // To avoid compiler warnings
341  return gradp ?
342  SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
343  (_c, f, p, z, _a) :
344  SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
345  (_c, f, p, z, _a);
346  break;
347  }
348  }
349 
350  /**
351  * @return the zeroth SphericalEngine::coeff object.
352  **********************************************************************/
354  { return _c[0]; }
355  };
356 
357 } // namespace GeographicLib
358 
359 #endif // GEOGRAPHICLIB_SPHERICALHARMONIC_HPP
Math::real operator()(real x, real y, real z) const
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:90
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
CircularEngine Circle(real p, real z, bool gradp) const
Math::real operator()(real x, real y, real z, real &gradx, real &grady, real &gradz) const
Package up coefficients for SphericalEngine.
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
const SphericalEngine::coeff & Coefficients() const
Header for GeographicLib::CircularEngine class.
Spherical harmonic sums for a circle.
Header for GeographicLib::Constants class.
Spherical harmonic series.
Header for GeographicLib::SphericalEngine class.
SphericalHarmonic(const std::vector< real > &C, const std::vector< real > &S, int N, real a, unsigned norm=FULL)
SphericalHarmonic(const std::vector< real > &C, const std::vector< real > &S, int N, int nmx, int mmx, real a, unsigned norm=FULL)