GeographicLib  1.21
SphericalEngine.hpp
Go to the documentation of this file.
00001 /**
00002  * \file SphericalEngine.hpp
00003  * \brief Header for GeographicLib::SphericalEngine class
00004  *
00005  * Copyright (c) Charles Karney (2011, 2012) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #if !defined(GEOGRAPHICLIB_SPHERICALENGINE_HPP)
00011 #define GEOGRAPHICLIB_SPHERICALENGINE_HPP \
00012   "$Id: f48320a694ecf901d997b23d32ea625e589f9534 $"
00013 
00014 #include <vector>
00015 #include <istream>
00016 #include <GeographicLib/Constants.hpp>
00017 
00018 #if defined(_MSC_VER)
00019 // Squelch warnings about dll vs vector
00020 #pragma warning (push)
00021 #pragma warning (disable: 4251)
00022 #endif
00023 
00024 namespace GeographicLib {
00025 
00026   class CircularEngine;
00027 
00028   /**
00029    * \brief The evaluation engine for SphericalHarmonic
00030    *
00031    * This serves as the backend to SphericalHarmonic, SphericalHarmonic1, and
00032    * SphericalHarmonic2.  Typically end-users will not have to access this
00033    * class directly.
00034    *
00035    * See SphericalEngine.cpp for more information on the implementation.
00036    *
00037    * Example of use:
00038    * \include example-SphericalEngine.cpp
00039    **********************************************************************/
00040 
00041   class GEOGRAPHIC_EXPORT SphericalEngine {
00042   private:
00043     typedef Math::real real;
00044     // A table of the square roots of integers
00045     static std::vector<real> root_;
00046     friend class CircularEngine; // CircularEngine needs access to root_, scale_
00047     // An internal scaling of the coefficients to avoid overflow in
00048     // intermediate calculations.
00049     static const real scale_;
00050     // Move latitudes near the pole off the axis by this amount.
00051     static const real eps_;
00052     static const std::vector<real> Z_;
00053     SphericalEngine();          // Disable constructor
00054   public:
00055     /**
00056      * Supported normalizations for associated Legendre polynomials.
00057      **********************************************************************/
00058     enum normalization {
00059       /**
00060        * Fully normalized associated Legendre polynomials.  See
00061        * SphericalHarmonic::FULL for documentation.
00062        *
00063        * @hideinitializer
00064        **********************************************************************/
00065       FULL = 0,
00066       /**
00067        * Schmidt semi-normalized associated Legendre polynomials.  See
00068        * SphericalHarmonic::SCHMIDT for documentation.
00069        *
00070        * @hideinitializer
00071        **********************************************************************/
00072       SCHMIDT = 1,
00073       /// \cond SKIP
00074       // These are deprecated...
00075       full = FULL,
00076       schmidt = SCHMIDT,
00077       /// \endcond
00078     };
00079 
00080     /**
00081      * \brief Package up coefficients for SphericalEngine
00082      *
00083      * This packages up the \e C, \e S coefficients and information about how
00084      * the coefficients are stored into a single structure.  This allows a
00085      * vector of type SphericalEngine::coeff to be passed to
00086      * SphericalEngine::Value.  This class also includes functions to aid
00087      * indexing into \e C and \e S.
00088      *
00089      * The storage layout of the coefficients is documented in
00090      * SphericalHarmonic and SphericalHarmonic::SphericalHarmonic.
00091      **********************************************************************/
00092     class GEOGRAPHIC_EXPORT coeff {
00093     private:
00094       int _N, _nmx, _mmx;
00095       std::vector<real>::const_iterator _Cnm;
00096       std::vector<real>::const_iterator _Snm;
00097     public:
00098       /**
00099        * A default constructor
00100        **********************************************************************/
00101       coeff()
00102         : _N(-1)
00103         , _nmx(-1)
00104         , _mmx(-1)
00105         , _Cnm(Z_.begin())
00106         , _Snm(Z_.begin()) {}
00107       /**
00108        * The general constructor.
00109        *
00110        * @param[in] C a vector of coefficients for the cosine terms.
00111        * @param[in] S a vector of coefficients for the sine terms.
00112        * @param[in] N the degree giving storage layout for \e C and \e S.
00113        * @param[in] nmx the maximum degree to be used.
00114        * @param[in] mmx the maximum order to be used.
00115        *
00116        * This requires \e N >= \e nmx >= \e mmx >= -1.  \e C and \e S must also
00117        * be large enough to hold the coefficients.  Otherwise an exception is
00118        * thrown.
00119        **********************************************************************/
00120       coeff(const std::vector<real>& C,
00121             const std::vector<real>& S,
00122             int N, int nmx, int mmx)
00123         : _N(N)
00124         , _nmx(nmx)
00125         , _mmx(mmx)
00126         , _Cnm(C.begin())
00127         , _Snm(S.begin())
00128       {
00129         if (!(_N >= _nmx && _nmx >= _mmx && _mmx >= -1))
00130           throw GeographicErr("Bad indices for coeff");
00131         if (!(index(_nmx, _mmx) < int(C.size()) &&
00132               index(_nmx, _mmx) < int(S.size()) + (_N + 1)))
00133           throw GeographicErr("Arrays too small in coeff");
00134         SphericalEngine::RootTable(_nmx);
00135       }
00136       /**
00137        * The constructor for full coefficient vectors.
00138        *
00139        * @param[in] C a vector of coefficients for the cosine terms.
00140        * @param[in] S a vector of coefficients for the sine terms.
00141        * @param[in] N the maximum degree and order.
00142        *
00143        * This requires \e N >= -1.  \e C and \e S must also be large enough to
00144        * hold the coefficients.  Otherwise an exception is thrown.
00145        **********************************************************************/
00146       coeff(const std::vector<real>& C,
00147             const std::vector<real>& S,
00148             int N)
00149         : _N(N)
00150         , _nmx(N)
00151         , _mmx(N)
00152         , _Cnm(C.begin())
00153         , _Snm(S.begin())
00154       {
00155         if (!(_N >= -1))
00156           throw GeographicErr("Bad indices for coeff");
00157         if (!(index(_nmx, _mmx) < int(C.size()) &&
00158               index(_nmx, _mmx) < int(S.size()) + (_N + 1)))
00159           throw GeographicErr("Arrays too small in coeff");
00160         SphericalEngine::RootTable(_nmx);
00161       }
00162       /**
00163        * @return \e N the degree giving storage layout for \e C and \e S.
00164        **********************************************************************/
00165       inline int N() const throw() { return _N; }
00166       /**
00167        * @return \e nmx the maximum degree to be used.
00168        **********************************************************************/
00169       inline int nmx() const throw() { return _nmx; }
00170       /**
00171        * @return \e mmx the maximum order to be used.
00172        **********************************************************************/
00173       inline int mmx() const throw() { return _mmx; }
00174       /**
00175        * The one-dimensional index into \e C and \e S.
00176        *
00177        * @param[in] n the degree.
00178        * @param[in] m the order.
00179        * @return the one-dimensional index.
00180        **********************************************************************/
00181       inline int index(int n, int m) const throw()
00182       { return m * _N - m * (m - 1) / 2 + n; }
00183       /**
00184        * An element of \e C.
00185        *
00186        * @param[in] k the one-dimensional index.
00187        * @return the value of the \e C coefficient.
00188        **********************************************************************/
00189       inline Math::real Cv(int k) const { return *(_Cnm + k); }
00190       /**
00191        * An element of \e S.
00192        *
00193        * @param[in] k the one-dimensional index.
00194        * @return the value of the \e S coefficient.
00195        **********************************************************************/
00196       inline Math::real Sv(int k) const { return *(_Snm + (k - (_N + 1))); }
00197       /**
00198        * An element of \e C with checking.
00199        *
00200        * @param[in] k the one-dimensional index.
00201        * @param[in] n the requested degree.
00202        * @param[in] m the requested order.
00203        * @param[in] f a multiplier.
00204        * @return the value of the \e C coefficient multiplied by \e f in \e n
00205        *   and \e m are in range else 0.
00206        **********************************************************************/
00207       inline Math::real Cv(int k, int n, int m, real f) const
00208       { return m > _mmx || n > _nmx ? 0 : *(_Cnm + k) * f; }
00209       /**
00210        * An element of \e S with checking.
00211        *
00212        * @param[in] k the one-dimensional index.
00213        * @param[in] n the requested degree.
00214        * @param[in] m the requested order.
00215        * @param[in] f a multiplier.
00216        * @return the value of the \e S coefficient multiplied by \e f in \e n
00217        *   and \e m are in range else 0.
00218        **********************************************************************/
00219       inline Math::real Sv(int k, int n, int m, real f) const
00220       { return m > _mmx || n > _nmx ? 0 : *(_Snm + (k - (_N + 1))) * f; }
00221 
00222       /**
00223        * The size of the coefficient vector for the cosine terms.
00224        *
00225        * @param[in] N the maximum degree.
00226        * @param[in] M the maximum order.
00227        * @return the size of the vector of cosine terms as stored in column
00228        *   major order.
00229        **********************************************************************/
00230       static inline int Csize(int N, int M)
00231       { return (M + 1) * (2 * N - M + 2) / 2; }
00232 
00233       /**
00234        * The size of the coefficient vector for the sine terms.
00235        *
00236        * @param[in] N the maximum degree.
00237        * @param[in] M the maximum order.
00238        * @return the size of the vector of cosine terms as stored in column
00239        *   major order.
00240        **********************************************************************/
00241       static inline int Ssize(int N, int M)
00242       { return Csize(N, M) - (N + 1); }
00243 
00244       /**
00245        * Load coefficients from a binary stream.
00246        *
00247        * @param[in] stream the input stream.
00248        * @param[out] N The maximum degree of the coefficients.
00249        * @param[out] M The maximum order of the coefficients.
00250        * @param[out] C The vector of cosine coefficients.
00251        * @param[out] S The vector of sine coefficients.
00252        *
00253        * \e N and \e M are read as 4-byte ints.  \e C and \e S are resized to
00254        * accommodate all the coefficients (with the \e m = 0 coefficients for
00255        * \e S excluded) and the data for these coefficients read as 8-byte
00256        * doubles.  The coefficients are stored in column major order.  The
00257        * bytes in the stream should use little-endian ordering.  IEEE floating
00258        * point is assumed for the coefficients.
00259        **********************************************************************/
00260       static void readcoeffs(std::istream& stream, int& N, int& M,
00261                              std::vector<real>& C, std::vector<real>& S);
00262     };
00263 
00264     /**
00265      * Evaluate a spherical harmonic sum and its gradient.
00266      *
00267      * @tparam gradp should the gradient be calculated.
00268      * @tparam norm the normalization for the associated Legendre polynomials.
00269      * @tparam L the number of terms in the coefficients.
00270      * @param[in] c an array of coeff objects.
00271      * @param[in] f array of coefficient multipliers.  f[0] should be 1.
00272      * @param[in] x the \e x component of the cartesian position.
00273      * @param[in] y the \e y component of the cartesian position.
00274      * @param[in] z the \e z component of the cartesian position.
00275      * @param[in] a the normalizing radius.
00276      * @param[out] gradx the \e x component of the gradient.
00277      * @param[out] grady the \e y component of the gradient.
00278      * @param[out] gradz the \e z component of the gradient.
00279      * @result the spherical harmonic sum.
00280      *
00281      * See the SphericalHarmonic class for the definition of the sum.
00282      * The coefficients used by this function are, for example,
00283      * c[0].Cv + f[1] * c[1].Cv + ... + f[L-1] * c[L-1].Cv.  (Note
00284      * that f[0] is \e not used.)  The upper limits on the sum are
00285      * determined by c[0].nmx() and c[0].mmx(); these limits apply to
00286      * \e all the components of the coefficients.  The parameters \e
00287      * gradp, \e norm, and \e L are template parameters, to allow more
00288      * optimization to be done at compile time.
00289      *
00290      * Clenshaw summation is used which permits the evaluation of the sum
00291      * without the need to allocate temporary arrays.  Thus this function never
00292      * throws an exception.
00293      **********************************************************************/
00294     template<bool gradp, normalization norm, int L>
00295       static Math::real Value(const coeff c[], const real f[],
00296                               real x, real y, real z, real a,
00297                               real& gradx, real& grady, real& gradz) throw();
00298 
00299     /**
00300      * Create a CircularEngine object
00301      *
00302      * @tparam gradp should the gradient be calculated.
00303      * @tparam norm the normalization for the associated Legendre polynomials.
00304      * @tparam L the number of terms in the coefficients.
00305      * @param[in] c an array of coeff objects.
00306      * @param[in] f array of coefficient multipliers.  f[0] should be 1.
00307      * @param[in] p the radius of the circle = sqrt(<i>x</i><sup>2</sup> + 
00308      *   <i>y</i><sup>2</sup>).
00309      * @param[in] z the height of the circle.
00310      * @param[in] a the normalizing radius.
00311      * @result the CircularEngine object.
00312      *
00313      * If you need to evaluate the spherical harmonic sum for several points
00314      * with constant \e f, \e p = sqrt(<i>x</i><sup>2</sup> +
00315      * <i>y</i><sup>2</sup>), \e z, and \e a, it is more efficient to construct
00316      * call SphericalEngine::Circle to give a CircularEngine object and then
00317      * call CircularEngine::operator()() with arguments <i>x</i>/\e p and
00318      * <i>y</i>/\e p.
00319      **********************************************************************/
00320     template<bool gradp, normalization norm, int L>
00321       static CircularEngine Circle(const coeff c[], const real f[],
00322                                    real p, real z, real a);
00323     /**
00324      * Check that the static table of square roots is big enough and enlarge it
00325      * if necessary.
00326      *
00327      * @param[in] N the maximum degree to be used in SphericalEngine.
00328      *
00329      * Typically, there's no need for an end-user to call this routine, because
00330      * the constructors for SphericalEngine::coeff do so.  However, since this
00331      * updates a static table, there's a possible race condition in a
00332      * multi-threaded environment.  Because this routine does nothing if the
00333      * table is already large enough, one way to avoid race conditions is to
00334      * call this routine at program start up (when it's still single threaded),
00335      * supplying the largest degree that your program will use.  E.g.,
00336      \code
00337   GeographicLib::SphericalEngine::RootTable(2190);
00338      \endcode
00339      * suffices to accommodate extant magnetic and gravity models.
00340      **********************************************************************/
00341     static void RootTable(int N);
00342 
00343     /**
00344      * Clear the static table of square roots and release the memory.  Call
00345      * this only when you are sure you no longer will be using SphericalEngine.
00346      * Your program will crash if you call SphericalEngine after calling this
00347      * routine.  <b>It's safest not to call this routine at all.</b> (The space
00348      * used by the table is modest.)
00349      **********************************************************************/
00350     static void ClearRootTable() {
00351       std::vector<real> temp(0);
00352       root_.swap(temp);
00353     }
00354   };
00355 
00356 } // namespace GeographicLib
00357 
00358 #if defined(_MSC_VER)
00359 #pragma warning (pop)
00360 #endif
00361 
00362 #endif  // GEOGRAPHICLIB_SPHERICALENGINE_HPP