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