1 /** 2 * \file SphericalHarmonic.hpp 3 * \brief Header for GeographicLib::SphericalHarmonic 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_SPHERICALHARMONIC_HPP) 11 #define GEOGRAPHICLIB_SPHERICALHARMONIC_HPP 1 12 13 #include <vector> 14 #include <GeographicLib/Constants.hpp> 15 #include <GeographicLib/SphericalEngine.hpp> 16 #include <GeographicLib/CircularEngine.hpp> 17 18 namespace GeographicLib { 19 20 /** 21 * \brief Spherical harmonic series 22 * 23 * This class evaluates the spherical harmonic sum \verbatim 24 V(x, y, z) = sum(n = 0..N)[ q^(n+1) * sum(m = 0..n)[ 25 (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) * 26 P[n,m](cos(theta)) ] ] 27 \endverbatim 28 * where 29 * - <i>p</i><sup>2</sup> = <i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>, 30 * - <i>r</i><sup>2</sup> = <i>p</i><sup>2</sup> + <i>z</i><sup>2</sup>, 31 * - \e q = <i>a</i>/<i>r</i>, 32 * - θ = atan2(\e p, \e z) = the spherical \e colatitude, 33 * - λ = atan2(\e y, \e x) = the longitude. 34 * - P<sub><i>nm</i></sub>(\e t) is the associated Legendre polynomial of 35 * degree \e n and order \e m. 36 * 37 * Two normalizations are supported for P<sub><i>nm</i></sub> 38 * - fully normalized denoted by SphericalHarmonic::FULL. 39 * - Schmidt semi-normalized denoted by SphericalHarmonic::SCHMIDT. 40 * 41 * Clenshaw summation is used for the sums over both \e n and \e m. This 42 * allows the computation to be carried out without the need for any 43 * temporary arrays. See SphericalEngine.cpp for more information on the 44 * implementation. 45 * 46 * References: 47 * - C. W. Clenshaw, 48 * <a href="https://doi.org/10.1090/S0025-5718-1955-0071856-0"> 49 * A note on the summation of Chebyshev series</a>, 50 * %Math. Tables Aids Comput. 9(51), 118--120 (1955). 51 * - R. E. Deakin, Derivatives of the earth's potentials, Geomatics 52 * Research Australasia 68, 31--60, (June 1998). 53 * - W. A. Heiskanen and H. Moritz, Physical Geodesy, (Freeman, San 54 * Francisco, 1967). (See Sec. 1-14, for a definition of Pbar.) 55 * - S. A. Holmes and W. E. Featherstone, 56 * <a href="https://doi.org/10.1007/s00190-002-0216-2"> 57 * A unified approach to the Clenshaw summation and the recursive 58 * computation of very high degree and order normalised associated Legendre 59 * functions</a>, J. Geodesy 76(5), 279--299 (2002). 60 * - C. C. Tscherning and K. Poder, 61 * <a href="http://cct.gfy.ku.dk/publ_cct/cct80.pdf"> 62 * Some geodetic applications of Clenshaw summation</a>, 63 * Boll. Geod. Sci. Aff. 41(4), 349--375 (1982). 64 * 65 * Example of use: 66 * \include example-SphericalHarmonic.cpp 67 **********************************************************************/ 68 69 class GEOGRAPHICLIB_EXPORT SphericalHarmonic { 70 public: 71 /** 72 * Supported normalizations for the associated Legendre polynomials. 73 **********************************************************************/ 74 enum normalization { 75 /** 76 * Fully normalized associated Legendre polynomials. 77 * 78 * These are defined by 79 * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(\e z) 80 * = (−1)<sup><i>m</i></sup> 81 * sqrt(\e k (2\e n + 1) (\e n − \e m)! / (\e n + \e m)!) 82 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where 83 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers 84 * function (also known as the Legendre function on the cut or the 85 * associated Legendre polynomial) https://dlmf.nist.gov/14.7.E10 and 86 * \e k = 1 for \e m = 0 and \e k = 2 otherwise. 87 * 88 * The mean squared value of 89 * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cosθ) 90 * cos(<i>m</i>λ) and 91 * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cosθ) 92 * sin(<i>m</i>λ) over the sphere is 1. 93 * 94 * @hideinitializer 95 **********************************************************************/ 96 FULL = SphericalEngine::FULL, 97 /** 98 * Schmidt semi-normalized associated Legendre polynomials. 99 * 100 * These are defined by 101 * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(\e z) 102 * = (−1)<sup><i>m</i></sup> 103 * sqrt(\e k (\e n − \e m)! / (\e n + \e m)!) 104 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where 105 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers 106 * function (also known as the Legendre function on the cut or the 107 * associated Legendre polynomial) https://dlmf.nist.gov/14.7.E10 and 108 * \e k = 1 for \e m = 0 and \e k = 2 otherwise. 109 * 110 * The mean squared value of 111 * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cosθ) 112 * cos(<i>m</i>λ) and 113 * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cosθ) 114 * sin(<i>m</i>λ) over the sphere is 1/(2\e n + 1). 115 * 116 * @hideinitializer 117 **********************************************************************/ 118 SCHMIDT = SphericalEngine::SCHMIDT, 119 }; 120 121 private: 122 typedef Math::real real; 123 SphericalEngine::coeff _c[1]; 124 real _a; 125 unsigned _norm; 126 127 public: 128 /** 129 * Constructor with a full set of coefficients specified. 130 * 131 * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>. 132 * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>. 133 * @param[in] N the maximum degree and order of the sum 134 * @param[in] a the reference radius appearing in the definition of the 135 * sum. 136 * @param[in] norm the normalization for the associated Legendre 137 * polynomials, either SphericalHarmonic::FULL (the default) or 138 * SphericalHarmonic::SCHMIDT. 139 * @exception GeographicErr if \e N does not satisfy \e N ≥ −1. 140 * @exception GeographicErr if \e C or \e S is not big enough to hold the 141 * coefficients. 142 * 143 * The coefficients <i>C</i><sub><i>nm</i></sub> and 144 * <i>S</i><sub><i>nm</i></sub> are stored in the one-dimensional vectors 145 * \e C and \e S which must contain (\e N + 1)(\e N + 2)/2 and \e N (\e N + 146 * 1)/2 elements, respectively, stored in "column-major" order. Thus for 147 * \e N = 3, the order would be: 148 * <i>C</i><sub>00</sub>, 149 * <i>C</i><sub>10</sub>, 150 * <i>C</i><sub>20</sub>, 151 * <i>C</i><sub>30</sub>, 152 * <i>C</i><sub>11</sub>, 153 * <i>C</i><sub>21</sub>, 154 * <i>C</i><sub>31</sub>, 155 * <i>C</i><sub>22</sub>, 156 * <i>C</i><sub>32</sub>, 157 * <i>C</i><sub>33</sub>. 158 * In general the (\e n,\e m) element is at index \e m \e N − \e m 159 * (\e m − 1)/2 + \e n. The layout of \e S is the same except that 160 * the first column is omitted (since the \e m = 0 terms never contribute 161 * to the sum) and the 0th element is <i>S</i><sub>11</sub> 162 * 163 * The class stores <i>pointers</i> to the first elements of \e C and \e S. 164 * These arrays should not be altered or destroyed during the lifetime of a 165 * SphericalHarmonic object. 166 **********************************************************************/ SphericalHarmonic(const std::vector<real> & C,const std::vector<real> & S,int N,real a,unsigned norm=FULL)167 SphericalHarmonic(const std::vector<real>& C, 168 const std::vector<real>& S, 169 int N, real a, unsigned norm = FULL) 170 : _a(a) 171 , _norm(norm) 172 { _c[0] = SphericalEngine::coeff(C, S, N); } 173 174 /** 175 * Constructor with a subset of coefficients specified. 176 * 177 * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>. 178 * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>. 179 * @param[in] N the degree used to determine the layout of \e C and \e S. 180 * @param[in] nmx the maximum degree used in the sum. The sum over \e n is 181 * from 0 thru \e nmx. 182 * @param[in] mmx the maximum order used in the sum. The sum over \e m is 183 * from 0 thru min(\e n, \e mmx). 184 * @param[in] a the reference radius appearing in the definition of the 185 * sum. 186 * @param[in] norm the normalization for the associated Legendre 187 * polynomials, either SphericalHarmonic::FULL (the default) or 188 * SphericalHarmonic::SCHMIDT. 189 * @exception GeographicErr if \e N, \e nmx, and \e mmx do not satisfy 190 * \e N ≥ \e nmx ≥ \e mmx ≥ −1. 191 * @exception GeographicErr if \e C or \e S is not big enough to hold the 192 * coefficients. 193 * 194 * The class stores <i>pointers</i> to the first elements of \e C and \e S. 195 * These arrays should not be altered or destroyed during the lifetime of a 196 * SphericalHarmonic object. 197 **********************************************************************/ SphericalHarmonic(const std::vector<real> & C,const std::vector<real> & S,int N,int nmx,int mmx,real a,unsigned norm=FULL)198 SphericalHarmonic(const std::vector<real>& C, 199 const std::vector<real>& S, 200 int N, int nmx, int mmx, 201 real a, unsigned norm = FULL) 202 : _a(a) 203 , _norm(norm) 204 { _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); } 205 206 /** 207 * A default constructor so that the object can be created when the 208 * constructor for another object is initialized. This default object can 209 * then be reset with the default copy assignment operator. 210 **********************************************************************/ SphericalHarmonic()211 SphericalHarmonic() {} 212 213 /** 214 * Compute the spherical harmonic sum. 215 * 216 * @param[in] x cartesian coordinate. 217 * @param[in] y cartesian coordinate. 218 * @param[in] z cartesian coordinate. 219 * @return \e V the spherical harmonic sum. 220 * 221 * This routine requires constant memory and thus never throws an 222 * exception. 223 **********************************************************************/ operator ()(real x,real y,real z) const224 Math::real operator()(real x, real y, real z) const { 225 real f[] = {1}; 226 real v = 0; 227 real dummy; 228 switch (_norm) { 229 case FULL: 230 v = SphericalEngine::Value<false, SphericalEngine::FULL, 1> 231 (_c, f, x, y, z, _a, dummy, dummy, dummy); 232 break; 233 case SCHMIDT: 234 default: // To avoid compiler warnings 235 v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1> 236 (_c, f, x, y, z, _a, dummy, dummy, dummy); 237 break; 238 } 239 return v; 240 } 241 242 /** 243 * Compute a spherical harmonic sum and its gradient. 244 * 245 * @param[in] x cartesian coordinate. 246 * @param[in] y cartesian coordinate. 247 * @param[in] z cartesian coordinate. 248 * @param[out] gradx \e x component of the gradient 249 * @param[out] grady \e y component of the gradient 250 * @param[out] gradz \e z component of the gradient 251 * @return \e V the spherical harmonic sum. 252 * 253 * This is the same as the previous function, except that the components of 254 * the gradients of the sum in the \e x, \e y, and \e z directions are 255 * computed. This routine requires constant memory and thus never throws 256 * an exception. 257 **********************************************************************/ operator ()(real x,real y,real z,real & gradx,real & grady,real & gradz) const258 Math::real operator()(real x, real y, real z, 259 real& gradx, real& grady, real& gradz) const { 260 real f[] = {1}; 261 real v = 0; 262 switch (_norm) { 263 case FULL: 264 v = SphericalEngine::Value<true, SphericalEngine::FULL, 1> 265 (_c, f, x, y, z, _a, gradx, grady, gradz); 266 break; 267 case SCHMIDT: 268 default: // To avoid compiler warnings 269 v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1> 270 (_c, f, x, y, z, _a, gradx, grady, gradz); 271 break; 272 } 273 return v; 274 } 275 276 /** 277 * Create a CircularEngine to allow the efficient evaluation of several 278 * points on a circle of latitude. 279 * 280 * @param[in] p the radius of the circle. 281 * @param[in] z the height of the circle above the equatorial plane. 282 * @param[in] gradp if true the returned object will be able to compute the 283 * gradient of the sum. 284 * @exception std::bad_alloc if the memory for the CircularEngine can't be 285 * allocated. 286 * @return the CircularEngine object. 287 * 288 * SphericalHarmonic::operator()() exchanges the order of the sums in the 289 * definition, i.e., ∑<sub><i>n</i> = 0..<i>N</i></sub> 290 * ∑<sub><i>m</i> = 0..<i>n</i></sub> becomes ∑<sub><i>m</i> = 291 * 0..<i>N</i></sub> ∑<sub><i>n</i> = <i>m</i>..<i>N</i></sub>. 292 * SphericalHarmonic::Circle performs the inner sum over degree \e n (which 293 * entails about <i>N</i><sup>2</sup> operations). Calling 294 * CircularEngine::operator()() on the returned object performs the outer 295 * sum over the order \e m (about \e N operations). 296 * 297 * Here's an example of computing the spherical sum at a sequence of 298 * longitudes without using a CircularEngine object \code 299 SphericalHarmonic h(...); // Create the SphericalHarmonic object 300 double r = 2, lat = 33, lon0 = 44, dlon = 0.01; 301 double 302 phi = lat * Math::degree<double>(), 303 z = r * sin(phi), p = r * cos(phi); 304 for (int i = 0; i <= 100; ++i) { 305 real 306 lon = lon0 + i * dlon, 307 lam = lon * Math::degree<double>(); 308 std::cout << lon << " " << h(p * cos(lam), p * sin(lam), z) << "\n"; 309 } 310 \endcode 311 * Here is the same calculation done using a CircularEngine object. This 312 * will be about <i>N</i>/2 times faster. \code 313 SphericalHarmonic h(...); // Create the SphericalHarmonic object 314 double r = 2, lat = 33, lon0 = 44, dlon = 0.01; 315 double 316 phi = lat * Math::degree<double>(), 317 z = r * sin(phi), p = r * cos(phi); 318 CircularEngine c(h(p, z, false)); // Create the CircularEngine object 319 for (int i = 0; i <= 100; ++i) { 320 real 321 lon = lon0 + i * dlon; 322 std::cout << lon << " " << c(lon) << "\n"; 323 } 324 \endcode 325 **********************************************************************/ Circle(real p,real z,bool gradp) const326 CircularEngine Circle(real p, real z, bool gradp) const { 327 real f[] = {1}; 328 switch (_norm) { 329 case FULL: 330 return gradp ? 331 SphericalEngine::Circle<true, SphericalEngine::FULL, 1> 332 (_c, f, p, z, _a) : 333 SphericalEngine::Circle<false, SphericalEngine::FULL, 1> 334 (_c, f, p, z, _a); 335 break; 336 case SCHMIDT: 337 default: // To avoid compiler warnings 338 return gradp ? 339 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1> 340 (_c, f, p, z, _a) : 341 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1> 342 (_c, f, p, z, _a); 343 break; 344 } 345 } 346 347 /** 348 * @return the zeroth SphericalEngine::coeff object. 349 **********************************************************************/ Coefficients() const350 const SphericalEngine::coeff& Coefficients() const 351 { return _c[0]; } 352 }; 353 354 } // namespace GeographicLib 355 356 #endif // GEOGRAPHICLIB_SPHERICALHARMONIC_HPP 357