GeographicLib  1.38
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SphericalEngine.hpp
Go to the documentation of this file.
1 /**
2  * \file SphericalEngine.hpp
3  * \brief Header for GeographicLib::SphericalEngine class
4  *
5  * Copyright (c) Charles Karney (2011-2012) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_SPHERICALENGINE_HPP)
11 #define GEOGRAPHICLIB_SPHERICALENGINE_HPP 1
12 
13 #include <vector>
14 #include <istream>
16 
17 #if defined(_MSC_VER)
18 // Squelch warnings about dll vs vector
19 # pragma warning (push)
20 # pragma warning (disable: 4251)
21 #endif
22 
23 namespace GeographicLib {
24 
25  class CircularEngine;
26 
27  /**
28  * \brief The evaluation engine for SphericalHarmonic
29  *
30  * This serves as the backend to SphericalHarmonic, SphericalHarmonic1, and
31  * SphericalHarmonic2. Typically end-users will not have to access this
32  * class directly.
33  *
34  * See SphericalEngine.cpp for more information on the implementation.
35  *
36  * Example of use:
37  * \include example-SphericalEngine.cpp
38  **********************************************************************/
39 
41  private:
42  typedef Math::real real;
43  // A table of the square roots of integers
44  static std::vector<real> root_;
45  friend class CircularEngine; // CircularEngine needs access to root_, scale_
46  // An internal scaling of the coefficients to avoid overflow in
47  // intermediate calculations.
48  static real scale() {
49  using std::pow;
50  return pow(real(std::numeric_limits<real>::radix),
51  -3 * (std::numeric_limits<real>::max_exponent < (1<<14) ?
52  std::numeric_limits<real>::max_exponent : (1<<14)) / 5);
53  }
54  // Move latitudes near the pole off the axis by this amount.
55  static real eps() {
56  using std::sqrt;
57  return std::numeric_limits<real>::epsilon() *
58  sqrt(std::numeric_limits<real>::epsilon());
59  }
60  static const std::vector<real> Z_;
61  SphericalEngine(); // Disable constructor
62  public:
63  /**
64  * Supported normalizations for associated Legendre polynomials.
65  **********************************************************************/
67  /**
68  * Fully normalized associated Legendre polynomials. See
69  * SphericalHarmonic::FULL for documentation.
70  *
71  * @hideinitializer
72  **********************************************************************/
73  FULL = 0,
74  /**
75  * Schmidt semi-normalized associated Legendre polynomials. See
76  * SphericalHarmonic::SCHMIDT for documentation.
77  *
78  * @hideinitializer
79  **********************************************************************/
80  SCHMIDT = 1,
81  /// \cond SKIP
82  // These are deprecated...
83  full = FULL,
84  schmidt = SCHMIDT,
85  /// \endcond
86  };
87 
88  /**
89  * \brief Package up coefficients for SphericalEngine
90  *
91  * This packages up the \e C, \e S coefficients and information about how
92  * the coefficients are stored into a single structure. This allows a
93  * vector of type SphericalEngine::coeff to be passed to
94  * SphericalEngine::Value. This class also includes functions to aid
95  * indexing into \e C and \e S.
96  *
97  * The storage layout of the coefficients is documented in
98  * SphericalHarmonic and SphericalHarmonic::SphericalHarmonic.
99  **********************************************************************/
101  private:
102  int _Nx, _nmx, _mmx;
103  std::vector<real>::const_iterator _Cnm;
104  std::vector<real>::const_iterator _Snm;
105  public:
106  /**
107  * A default constructor
108  **********************************************************************/
110  : _Nx(-1)
111  , _nmx(-1)
112  , _mmx(-1)
113  , _Cnm(Z_.begin())
114  , _Snm(Z_.begin()) {}
115  /**
116  * The general constructor.
117  *
118  * @param[in] C a vector of coefficients for the cosine terms.
119  * @param[in] S a vector of coefficients for the sine terms.
120  * @param[in] N the degree giving storage layout for \e C and \e S.
121  * @param[in] nmx the maximum degree to be used.
122  * @param[in] mmx the maximum order to be used.
123  * @exception GeographicErr if \e N, \e nmx, and \e mmx do not satisfy
124  * \e N &ge; \e nmx &ge; \e mmx &ge; &minus;1.
125  * @exception GeographicErr if \e C or \e S is not big enough to hold the
126  * coefficients.
127  * @exception std::bad_alloc if the memory for the square root table
128  * can't be allocated.
129  **********************************************************************/
130  coeff(const std::vector<real>& C,
131  const std::vector<real>& S,
132  int N, int nmx, int mmx)
133  : _Nx(N)
134  , _nmx(nmx)
135  , _mmx(mmx)
136  , _Cnm(C.begin())
137  , _Snm(S.begin())
138  {
139  if (!(_Nx >= _nmx && _nmx >= _mmx && _mmx >= -1))
140  throw GeographicErr("Bad indices for coeff");
141  if (!(index(_nmx, _mmx) < int(C.size()) &&
142  index(_nmx, _mmx) < int(S.size()) + (_Nx + 1)))
143  throw GeographicErr("Arrays too small in coeff");
145  }
146  /**
147  * The constructor for full coefficient vectors.
148  *
149  * @param[in] C a vector of coefficients for the cosine terms.
150  * @param[in] S a vector of coefficients for the sine terms.
151  * @param[in] N the maximum degree and order.
152  * @exception GeographicErr if \e N does not satisfy \e N &ge; &minus;1.
153  * @exception GeographicErr if \e C or \e S is not big enough to hold the
154  * coefficients.
155  * @exception std::bad_alloc if the memory for the square root table
156  * can't be allocated.
157  **********************************************************************/
158  coeff(const std::vector<real>& C,
159  const std::vector<real>& S,
160  int N)
161  : _Nx(N)
162  , _nmx(N)
163  , _mmx(N)
164  , _Cnm(C.begin())
165  , _Snm(S.begin())
166  {
167  if (!(_Nx >= -1))
168  throw GeographicErr("Bad indices for coeff");
169  if (!(index(_nmx, _mmx) < int(C.size()) &&
170  index(_nmx, _mmx) < int(S.size()) + (_Nx + 1)))
171  throw GeographicErr("Arrays too small in coeff");
173  }
174  /**
175  * @return \e N the degree giving storage layout for \e C and \e S.
176  **********************************************************************/
177  inline int N() const { return _Nx; }
178  /**
179  * @return \e nmx the maximum degree to be used.
180  **********************************************************************/
181  inline int nmx() const { return _nmx; }
182  /**
183  * @return \e mmx the maximum order to be used.
184  **********************************************************************/
185  inline int mmx() const { return _mmx; }
186  /**
187  * The one-dimensional index into \e C and \e S.
188  *
189  * @param[in] n the degree.
190  * @param[in] m the order.
191  * @return the one-dimensional index.
192  **********************************************************************/
193  inline int index(int n, int m) const
194  { return m * _Nx - m * (m - 1) / 2 + n; }
195  /**
196  * An element of \e C.
197  *
198  * @param[in] k the one-dimensional index.
199  * @return the value of the \e C coefficient.
200  **********************************************************************/
201  inline Math::real Cv(int k) const { return *(_Cnm + k); }
202  /**
203  * An element of \e S.
204  *
205  * @param[in] k the one-dimensional index.
206  * @return the value of the \e S coefficient.
207  **********************************************************************/
208  inline Math::real Sv(int k) const { return *(_Snm + (k - (_Nx + 1))); }
209  /**
210  * An element of \e C with checking.
211  *
212  * @param[in] k the one-dimensional index.
213  * @param[in] n the requested degree.
214  * @param[in] m the requested order.
215  * @param[in] f a multiplier.
216  * @return the value of the \e C coefficient multiplied by \e f in \e n
217  * and \e m are in range else 0.
218  **********************************************************************/
219  inline Math::real Cv(int k, int n, int m, real f) const
220  { return m > _mmx || n > _nmx ? 0 : *(_Cnm + k) * f; }
221  /**
222  * An element of \e S with checking.
223  *
224  * @param[in] k the one-dimensional index.
225  * @param[in] n the requested degree.
226  * @param[in] m the requested order.
227  * @param[in] f a multiplier.
228  * @return the value of the \e S coefficient multiplied by \e f in \e n
229  * and \e m are in range else 0.
230  **********************************************************************/
231  inline Math::real Sv(int k, int n, int m, real f) const
232  { return m > _mmx || n > _nmx ? 0 : *(_Snm + (k - (_Nx + 1))) * f; }
233 
234  /**
235  * The size of the coefficient vector for the cosine terms.
236  *
237  * @param[in] N the maximum degree.
238  * @param[in] M the maximum order.
239  * @return the size of the vector of cosine terms as stored in column
240  * major order.
241  **********************************************************************/
242  static inline int Csize(int N, int M)
243  { return (M + 1) * (2 * N - M + 2) / 2; }
244 
245  /**
246  * The size of the coefficient vector for the sine terms.
247  *
248  * @param[in] N the maximum degree.
249  * @param[in] M the maximum order.
250  * @return the size of the vector of cosine terms as stored in column
251  * major order.
252  **********************************************************************/
253  static inline int Ssize(int N, int M)
254  { return Csize(N, M) - (N + 1); }
255 
256  /**
257  * Load coefficients from a binary stream.
258  *
259  * @param[in] stream the input stream.
260  * @param[out] N The maximum degree of the coefficients.
261  * @param[out] M The maximum order of the coefficients.
262  * @param[out] C The vector of cosine coefficients.
263  * @param[out] S The vector of sine coefficients.
264  * @exception GeographicErr if \e N and \e M do not satisfy \e N &ge;
265  * \e M &ge; &minus;1.
266  * @exception GeographicErr if there's an error reading the data.
267  * @exception std::bad_alloc if the memory for \e C or \e S can't be
268  * allocated.
269  *
270  * \e N and \e M are read as 4-byte ints. \e C and \e S are resized to
271  * accommodate all the coefficients (with the \e m = 0 coefficients for
272  * \e S excluded) and the data for these coefficients read as 8-byte
273  * doubles. The coefficients are stored in column major order. The
274  * bytes in the stream should use little-endian ordering. IEEE floating
275  * point is assumed for the coefficients.
276  **********************************************************************/
277  static void readcoeffs(std::istream& stream, int& N, int& M,
278  std::vector<real>& C, std::vector<real>& S);
279  };
280 
281  /**
282  * Evaluate a spherical harmonic sum and its gradient.
283  *
284  * @tparam gradp should the gradient be calculated.
285  * @tparam norm the normalization for the associated Legendre polynomials.
286  * @tparam L the number of terms in the coefficients.
287  * @param[in] c an array of coeff objects.
288  * @param[in] f array of coefficient multipliers. f[0] should be 1.
289  * @param[in] x the \e x component of the cartesian position.
290  * @param[in] y the \e y component of the cartesian position.
291  * @param[in] z the \e z component of the cartesian position.
292  * @param[in] a the normalizing radius.
293  * @param[out] gradx the \e x component of the gradient.
294  * @param[out] grady the \e y component of the gradient.
295  * @param[out] gradz the \e z component of the gradient.
296  * @result the spherical harmonic sum.
297  *
298  * See the SphericalHarmonic class for the definition of the sum. The
299  * coefficients used by this function are, for example, c[0].Cv + f[1] *
300  * c[1].Cv + ... + f[L&minus;1] * c[L&minus;1].Cv. (Note that f[0] is \e
301  * not used.) The upper limits on the sum are determined by c[0].nmx() and
302  * c[0].mmx(); these limits apply to \e all the components of the
303  * coefficients. The parameters \e gradp, \e norm, and \e L are template
304  * parameters, to allow more optimization to be done at compile time.
305  *
306  * Clenshaw summation is used which permits the evaluation of the sum
307  * without the need to allocate temporary arrays. Thus this function never
308  * throws an exception.
309  **********************************************************************/
310  template<bool gradp, normalization norm, int L>
311  static Math::real Value(const coeff c[], const real f[],
312  real x, real y, real z, real a,
313  real& gradx, real& grady, real& gradz);
314 
315  /**
316  * Create a CircularEngine object
317  *
318  * @tparam gradp should the gradient be calculated.
319  * @tparam norm the normalization for the associated Legendre polynomials.
320  * @tparam L the number of terms in the coefficients.
321  * @param[in] c an array of coeff objects.
322  * @param[in] f array of coefficient multipliers. f[0] should be 1.
323  * @param[in] p the radius of the circle = sqrt(<i>x</i><sup>2</sup> +
324  * <i>y</i><sup>2</sup>).
325  * @param[in] z the height of the circle.
326  * @param[in] a the normalizing radius.
327  * @exception std::bad_alloc if the memory for the CircularEngine can't be
328  * allocated.
329  * @result the CircularEngine object.
330  *
331  * If you need to evaluate the spherical harmonic sum for several points
332  * with constant \e f, \e p = sqrt(<i>x</i><sup>2</sup> +
333  * <i>y</i><sup>2</sup>), \e z, and \e a, it is more efficient to construct
334  * call SphericalEngine::Circle to give a CircularEngine object and then
335  * call CircularEngine::operator()() with arguments <i>x</i>/\e p and
336  * <i>y</i>/\e p.
337  **********************************************************************/
338  template<bool gradp, normalization norm, int L>
339  static CircularEngine Circle(const coeff c[], const real f[],
340  real p, real z, real a);
341  /**
342  * Check that the static table of square roots is big enough and enlarge it
343  * if necessary.
344  *
345  * @param[in] N the maximum degree to be used in SphericalEngine.
346  * @exception std::bad_alloc if the memory for the square root table can't
347  * be allocated.
348  *
349  * Typically, there's no need for an end-user to call this routine, because
350  * the constructors for SphericalEngine::coeff do so. However, since this
351  * updates a static table, there's a possible race condition in a
352  * multi-threaded environment. Because this routine does nothing if the
353  * table is already large enough, one way to avoid race conditions is to
354  * call this routine at program start up (when it's still single threaded),
355  * supplying the largest degree that your program will use. E.g., \code
356  GeographicLib::SphericalEngine::RootTable(2190);
357  \endcode
358  * suffices to accommodate extant magnetic and gravity models.
359  **********************************************************************/
360  static void RootTable(int N);
361 
362  /**
363  * Clear the static table of square roots and release the memory. Call
364  * this only when you are sure you no longer will be using SphericalEngine.
365  * Your program will crash if you call SphericalEngine after calling this
366  * routine. <b>It's safest not to call this routine at all.</b> (The space
367  * used by the table is modest.)
368  **********************************************************************/
369  static void ClearRootTable() {
370  std::vector<real> temp(0);
371  root_.swap(temp);
372  }
373  };
374 
375 } // namespace GeographicLib
376 
377 #if defined(_MSC_VER)
378 # pragma warning (pop)
379 #endif
380 
381 #endif // GEOGRAPHICLIB_SPHERICALENGINE_HPP