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