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 ≥ \e nmx ≥ \e mmx ≥ −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 ≥ −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 ≥ 265 * \e M ≥ −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−1] * c[L−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