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