GeographicLib
1.21
|
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