GeographicLib  2.0
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-2019) <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 _nNx, _nmx, _mmx;
102  std::vector<real>::const_iterator _cCnm;
103  std::vector<real>::const_iterator _sSnm;
104  public:
105  /**
106  * A default constructor
107  **********************************************************************/
108  coeff() : _nNx(-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  : _nNx(N)
128  , _nmx(nmx)
129  , _mmx(mmx)
130  , _cCnm(C.begin())
131  , _sSnm(S.begin())
132  {
133  if (!((_nNx >= _nmx && _nmx >= _mmx && _mmx >= 0) ||
134  // If mmx = -1 then the sums are empty so require nmx = -1 also.
135  (_nmx == -1 && _mmx == -1)))
136  throw GeographicErr("Bad indices for coeff");
137  if (!(index(_nmx, _mmx) < int(C.size()) &&
138  index(_nmx, _mmx) < int(S.size()) + (_nNx + 1)))
139  throw GeographicErr("Arrays too small in coeff");
141  }
142  /**
143  * The constructor for full coefficient vectors.
144  *
145  * @param[in] C a vector of coefficients for the cosine terms.
146  * @param[in] S a vector of coefficients for the sine terms.
147  * @param[in] N the maximum degree and order.
148  * @exception GeographicErr if \e N does not satisfy \e N &ge; &minus;1.
149  * @exception GeographicErr if \e C or \e S is not big enough to hold the
150  * coefficients.
151  * @exception std::bad_alloc if the memory for the square root table
152  * can't be allocated.
153  **********************************************************************/
154  coeff(const std::vector<real>& C,
155  const std::vector<real>& S,
156  int N)
157  : _nNx(N)
158  , _nmx(N)
159  , _mmx(N)
160  , _cCnm(C.begin())
161  , _sSnm(S.begin())
162  {
163  if (!(_nNx >= -1))
164  throw GeographicErr("Bad indices for coeff");
165  if (!(index(_nmx, _mmx) < int(C.size()) &&
166  index(_nmx, _mmx) < int(S.size()) + (_nNx + 1)))
167  throw GeographicErr("Arrays too small in coeff");
169  }
170  /**
171  * @return \e N the degree giving storage layout for \e C and \e S.
172  **********************************************************************/
173  int N() const { return _nNx; }
174  /**
175  * @return \e nmx the maximum degree to be used.
176  **********************************************************************/
177  int nmx() const { return _nmx; }
178  /**
179  * @return \e mmx the maximum order to be used.
180  **********************************************************************/
181  int mmx() const { return _mmx; }
182  /**
183  * The one-dimensional index into \e C and \e S.
184  *
185  * @param[in] n the degree.
186  * @param[in] m the order.
187  * @return the one-dimensional index.
188  **********************************************************************/
189  int index(int n, int m) const
190  { return m * _nNx - m * (m - 1) / 2 + n; }
191  /**
192  * An element of \e C.
193  *
194  * @param[in] k the one-dimensional index.
195  * @return the value of the \e C coefficient.
196  **********************************************************************/
197  Math::real Cv(int k) const { return *(_cCnm + k); }
198  /**
199  * An element of \e S.
200  *
201  * @param[in] k the one-dimensional index.
202  * @return the value of the \e S coefficient.
203  **********************************************************************/
204  Math::real Sv(int k) const { return *(_sSnm + (k - (_nNx + 1))); }
205  /**
206  * An element of \e C with checking.
207  *
208  * @param[in] k the one-dimensional index.
209  * @param[in] n the requested degree.
210  * @param[in] m the requested order.
211  * @param[in] f a multiplier.
212  * @return the value of the \e C coefficient multiplied by \e f in \e n
213  * and \e m are in range else 0.
214  **********************************************************************/
215  Math::real Cv(int k, int n, int m, real f) const
216  { return m > _mmx || n > _nmx ? 0 : *(_cCnm + k) * f; }
217  /**
218  * An element of \e S with checking.
219  *
220  * @param[in] k the one-dimensional index.
221  * @param[in] n the requested degree.
222  * @param[in] m the requested order.
223  * @param[in] f a multiplier.
224  * @return the value of the \e S coefficient multiplied by \e f in \e n
225  * and \e m are in range else 0.
226  **********************************************************************/
227  Math::real Sv(int k, int n, int m, real f) const
228  { return m > _mmx || n > _nmx ? 0 : *(_sSnm + (k - (_nNx + 1))) * f; }
229 
230  /**
231  * The size of the coefficient vector for the cosine terms.
232  *
233  * @param[in] N the maximum degree.
234  * @param[in] M the maximum order.
235  * @return the size of the vector of cosine terms as stored in column
236  * major order.
237  **********************************************************************/
238  static int Csize(int N, int M)
239  { return (M + 1) * (2 * N - M + 2) / 2; }
240 
241  /**
242  * The size of the coefficient vector for the sine terms.
243  *
244  * @param[in] N the maximum degree.
245  * @param[in] M the maximum order.
246  * @return the size of the vector of cosine terms as stored in column
247  * major order.
248  **********************************************************************/
249  static int Ssize(int N, int M)
250  { return Csize(N, M) - (N + 1); }
251 
252  /**
253  * Load coefficients from a binary stream.
254  *
255  * @param[in] stream the input stream.
256  * @param[in,out] N The maximum degree of the coefficients.
257  * @param[in,out] M The maximum order of the coefficients.
258  * @param[out] C The vector of cosine coefficients.
259  * @param[out] S The vector of sine coefficients.
260  * @param[in] truncate if false (the default) then \e N and \e M are
261  * determined by the values in the binary stream; otherwise, the input
262  * values of \e N and \e M are used to truncate the coefficients read
263  * from the stream at the given degree and order.
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  bool truncate = false);
280  };
281 
282  /**
283  * Evaluate a spherical harmonic sum and its gradient.
284  *
285  * @tparam gradp should the gradient be calculated.
286  * @tparam norm the normalization for the associated Legendre polynomials.
287  * @tparam L the number of terms in the coefficients.
288  * @param[in] c an array of coeff objects.
289  * @param[in] f array of coefficient multipliers. f[0] should be 1.
290  * @param[in] x the \e x component of the cartesian position.
291  * @param[in] y the \e y component of the cartesian position.
292  * @param[in] z the \e z component of the cartesian position.
293  * @param[in] a the normalizing radius.
294  * @param[out] gradx the \e x component of the gradient.
295  * @param[out] grady the \e y component of the gradient.
296  * @param[out] gradz the \e z component of the gradient.
297  * @result the spherical harmonic sum.
298  *
299  * See the SphericalHarmonic class for the definition of the sum. The
300  * coefficients used by this function are, for example, c[0].Cv + f[1] *
301  * c[1].Cv + ... + f[L&minus;1] * c[L&minus;1].Cv. (Note that f[0] is \e
302  * not used.) The upper limits on the sum are determined by c[0].nmx() and
303  * c[0].mmx(); these limits apply to \e all the components of the
304  * coefficients. The parameters \e gradp, \e norm, and \e L are template
305  * parameters, to allow more optimization to be done at compile time.
306  *
307  * Clenshaw summation is used which permits the evaluation of the sum
308  * without the need to allocate temporary arrays. Thus this function never
309  * throws an exception.
310  **********************************************************************/
311  template<bool gradp, normalization norm, int L>
312  static Math::real Value(const coeff c[], const real f[],
313  real x, real y, real z, real a,
314  real& gradx, real& grady, real& gradz);
315 
316  /**
317  * Create a CircularEngine object
318  *
319  * @tparam gradp should the gradient be calculated.
320  * @tparam norm the normalization for the associated Legendre polynomials.
321  * @tparam L the number of terms in the coefficients.
322  * @param[in] c an array of coeff objects.
323  * @param[in] f array of coefficient multipliers. f[0] should be 1.
324  * @param[in] p the radius of the circle = sqrt(<i>x</i><sup>2</sup> +
325  * <i>y</i><sup>2</sup>).
326  * @param[in] z the height of the circle.
327  * @param[in] a the normalizing radius.
328  * @exception std::bad_alloc if the memory for the CircularEngine can't be
329  * allocated.
330  * @result the CircularEngine object.
331  *
332  * If you need to evaluate the spherical harmonic sum for several points
333  * with constant \e f, \e p = sqrt(<i>x</i><sup>2</sup> +
334  * <i>y</i><sup>2</sup>), \e z, and \e a, it is more efficient to construct
335  * call SphericalEngine::Circle to give a CircularEngine object and then
336  * call CircularEngine::operator()() with arguments <i>x</i>/\e p and
337  * <i>y</i>/\e p.
338  **********************************************************************/
339  template<bool gradp, normalization norm, int L>
340  static CircularEngine Circle(const coeff c[], const real f[],
341  real p, real z, real a);
342  /**
343  * Check that the static table of square roots is big enough and enlarge it
344  * if necessary.
345  *
346  * @param[in] N the maximum degree to be used in SphericalEngine.
347  * @exception std::bad_alloc if the memory for the square root table can't
348  * be allocated.
349  *
350  * Typically, there's no need for an end-user to call this routine, because
351  * the constructors for SphericalEngine::coeff do so. However, since this
352  * updates a static table, there's a possible race condition in a
353  * multi-threaded environment. Because this routine does nothing if the
354  * table is already large enough, one way to avoid race conditions is to
355  * call this routine at program start up (when it's still single threaded),
356  * supplying the largest degree that your program will use. E.g., \code
357  GeographicLib::SphericalEngine::RootTable(2190);
358  \endcode
359  * suffices to accommodate extant magnetic and gravity models.
360  **********************************************************************/
361  static void RootTable(int N);
362 
363  /**
364  * Clear the static table of square roots and release the memory. Call
365  * this only when you are sure you no longer will be using SphericalEngine.
366  * Your program will crash if you call SphericalEngine after calling this
367  * routine.
368  *
369  * \warning It's safest not to call this routine at all. (The space used
370  * by the table is modest.)
371  **********************************************************************/
372  static void ClearRootTable() {
373  std::vector<real> temp(0);
374  sqrttable().swap(temp);
375  }
376  };
377 
378 } // namespace GeographicLib
379 
380 #if defined(_MSC_VER)
381 # pragma warning (pop)
382 #endif
383 
384 #endif // GEOGRAPHICLIB_SPHERICALENGINE_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Spherical harmonic sums for a circle.
Exception handling for GeographicLib.
Definition: Constants.hpp:316
Package up coefficients for SphericalEngine.
coeff(const std::vector< real > &C, const std::vector< real > &S, int N)
coeff(const std::vector< real > &C, const std::vector< real > &S, int N, int nmx, int mmx)
Math::real Sv(int k, int n, int m, real f) const
Math::real Cv(int k, int n, int m, real f) const
The evaluation engine for SphericalHarmonic.
Namespace for GeographicLib.
Definition: Accumulator.cpp:12