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>
15 #include <GeographicLib/Constants.hpp>
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 
40   class GEOGRAPHICLIB_EXPORT SphericalEngine {
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.
scale()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.
eps()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      **********************************************************************/
70     enum normalization {
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      **********************************************************************/
99     class GEOGRAPHICLIB_EXPORT coeff {
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        **********************************************************************/
coeff()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        **********************************************************************/
coeff(const std::vector<real> & C,const std::vector<real> & S,int N,int nmx,int mmx)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 >= 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()) + (_Nx + 1)))
139           throw GeographicErr("Arrays too small in coeff");
140         SphericalEngine::RootTable(_nmx);
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        **********************************************************************/
coeff(const std::vector<real> & C,const std::vector<real> & S,int N)154       coeff(const std::vector<real>& C,
155             const std::vector<real>& S,
156             int N)
157         : _Nx(N)
158         , _nmx(N)
159         , _mmx(N)
160         , _Cnm(C.begin())
161         , _Snm(S.begin())
162       {
163         if (!(_Nx >= -1))
164           throw GeographicErr("Bad indices for coeff");
165         if (!(index(_nmx, _mmx) < int(C.size()) &&
166               index(_nmx, _mmx) < int(S.size()) + (_Nx + 1)))
167           throw GeographicErr("Arrays too small in coeff");
168         SphericalEngine::RootTable(_nmx);
169       }
170       /**
171        * @return \e N the degree giving storage layout for \e C and \e S.
172        **********************************************************************/
N() const173       int N() const { return _Nx; }
174       /**
175        * @return \e nmx the maximum degree to be used.
176        **********************************************************************/
nmx() const177       int nmx() const { return _nmx; }
178       /**
179        * @return \e mmx the maximum order to be used.
180        **********************************************************************/
mmx() const181       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        **********************************************************************/
index(int n,int m) const189       int index(int n, int m) const
190       { return m * _Nx - 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        **********************************************************************/
Cv(int k) const197       Math::real Cv(int k) const { return *(_Cnm + 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        **********************************************************************/
Sv(int k) const204       Math::real Sv(int k) const { return *(_Snm + (k - (_Nx + 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        **********************************************************************/
Cv(int k,int n,int m,real f) const215       Math::real Cv(int k, int n, int m, real f) const
216       { return m > _mmx || n > _nmx ? 0 : *(_Cnm + 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        **********************************************************************/
Sv(int k,int n,int m,real f) const227       Math::real Sv(int k, int n, int m, real f) const
228       { return m > _mmx || n > _nmx ? 0 : *(_Snm + (k - (_Nx + 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        **********************************************************************/
Csize(int N,int M)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        **********************************************************************/
Ssize(int N,int M)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      **********************************************************************/
ClearRootTable()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
385