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