00001 /** 00002 * \file SphericalHarmonic.hpp 00003 * \brief Header for GeographicLib::SphericalHarmonic class 00004 * 00005 * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under 00006 * the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC_HPP) 00011 #define GEOGRAPHICLIB_SPHERICALHARMONIC_HPP 1 00012 00013 #include <vector> 00014 #include <GeographicLib/Constants.hpp> 00015 #include <GeographicLib/SphericalEngine.hpp> 00016 #include <GeographicLib/CircularEngine.hpp> 00017 #include <GeographicLib/Geocentric.hpp> 00018 00019 namespace GeographicLib { 00020 00021 /** 00022 * \brief Spherical harmonic series 00023 * 00024 * This class evaluates the spherical harmonic sum \verbatim 00025 V(x, y, z) = sum(n = 0..N)[ q^(n+1) * sum(m = 0..n)[ 00026 (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) * 00027 P[n,m](cos(theta)) ] ] 00028 \endverbatim 00029 * where 00030 * - <i>p</i><sup>2</sup> = <i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>, 00031 * - <i>r</i><sup>2</sup> = <i>p</i><sup>2</sup> + <i>z</i><sup>2</sup>, 00032 * - \e q = <i>a</i>/<i>r</i>, 00033 * - θ = atan2(\e p, \e z) = the spherical \e colatitude, 00034 * - λ = atan2(\e y, \e x) = the longitude. 00035 * - P<sub><i>nm</i></sub>(\e t) is the associated Legendre polynomial of 00036 * degree \e n and order \e m. 00037 * 00038 * Two normalizations are supported for P<sub><i>nm</i></sub> 00039 * - fully normalized denoted by SphericalHarmonic::FULL. 00040 * - Schmidt semi-normalized denoted by SphericalHarmonic::SCHMIDT. 00041 * 00042 * Clenshaw summation is used for the sums over both \e n and \e m. This 00043 * allows the computation to be carried out without the need for any 00044 * temporary arrays. See SphericalEngine.cpp for more information on the 00045 * implementation. 00046 * 00047 * References: 00048 * - C. W. Clenshaw, A note on the summation of Chebyshev series, 00049 * %Math. Tables Aids Comput. 9(51), 118--120 (1955). 00050 * - R. E. Deakin, Derivatives of the earth's potentials, Geomatics 00051 * Research Australasia 68, 31--60, (June 1998). 00052 * - W. A. Heiskanen and H. Moritz, Physical Geodesy, (Freeman, San 00053 * Francisco, 1967). (See Sec. 1-14, for a definition of Pbar.) 00054 * - S. A. Holmes and W. E. Featherstone, A unified approach to the Clenshaw 00055 * summation and the recursive computation of very high degree and order 00056 * normalised associated Legendre functions, J. Geodesy 76(5), 00057 * 279--299 (2002). 00058 * - C. C. Tscherning and K. Poder, Some geodetic applications of Clenshaw 00059 * summation, Boll. Geod. Sci. Aff. 41(4), 349--375 (1982). 00060 * 00061 * Example of use: 00062 * \include example-SphericalHarmonic.cpp 00063 **********************************************************************/ 00064 00065 class GEOGRAPHICLIB_EXPORT SphericalHarmonic { 00066 public: 00067 /** 00068 * Supported normalizations for the associated Legendre polynomials. 00069 **********************************************************************/ 00070 enum normalization { 00071 /** 00072 * Fully normalized associated Legendre polynomials. 00073 * 00074 * These are defined by 00075 * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(\e z) 00076 * = (−1)<sup><i>m</i></sup> 00077 * sqrt(\e k (2\e n + 1) (\e n − \e m)! / (\e n + \e m)!) 00078 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where 00079 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers 00080 * function (also known as the Legendre function on the cut or the 00081 * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k 00082 * = 1 for \e m = 0 and \e k = 2 otherwise. 00083 * 00084 * The mean squared value of 00085 * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cosθ) 00086 * cos(<i>m</i>λ) and 00087 * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cosθ) 00088 * sin(<i>m</i>λ) over the sphere is 1. 00089 * 00090 * @hideinitializer 00091 **********************************************************************/ 00092 FULL = SphericalEngine::FULL, 00093 /** 00094 * Schmidt semi-normalized associated Legendre polynomials. 00095 * 00096 * These are defined by 00097 * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(\e z) 00098 * = (−1)<sup><i>m</i></sup> 00099 * sqrt(\e k (\e n − \e m)! / (\e n + \e m)!) 00100 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where 00101 * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers 00102 * function (also known as the Legendre function on the cut or the 00103 * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k 00104 * = 1 for \e m = 0 and \e k = 2 otherwise. 00105 * 00106 * The mean squared value of 00107 * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cosθ) 00108 * cos(<i>m</i>λ) and 00109 * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cosθ) 00110 * sin(<i>m</i>λ) over the sphere is 1/(2\e n + 1). 00111 * 00112 * @hideinitializer 00113 **********************************************************************/ 00114 SCHMIDT = SphericalEngine::SCHMIDT, 00115 /// \cond SKIP 00116 // These are deprecated... 00117 full = FULL, 00118 schmidt = SCHMIDT, 00119 /// \endcond 00120 }; 00121 00122 private: 00123 typedef Math::real real; 00124 SphericalEngine::coeff _c[1]; 00125 real _a; 00126 unsigned _norm; 00127 00128 public: 00129 /** 00130 * Constructor with a full set of coefficients specified. 00131 * 00132 * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>. 00133 * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>. 00134 * @param[in] N the maximum degree and order of the sum 00135 * @param[in] a the reference radius appearing in the definition of the 00136 * sum. 00137 * @param[in] norm the normalization for the associated Legendre 00138 * polynomials, either SphericalHarmonic::full (the default) or 00139 * SphericalHarmonic::schmidt. 00140 * @exception GeographicErr if \e N does not satisfy \e N ≥ −1. 00141 * @exception GeographicErr if \e C or \e S is not big enough to hold the 00142 * coefficients. 00143 * 00144 * The coefficients <i>C</i><sub><i>nm</i></sub> and 00145 * <i>S</i><sub><i>nm</i></sub> are stored in the one-dimensional vectors 00146 * \e C and \e S which must contain (\e N + 1)(\e N + 2)/2 and \e N (\e N + 00147 * 1)/2 elements, respectively, stored in "column-major" order. Thus for 00148 * \e N = 3, the order would be: 00149 * <i>C</i><sub>00</sub>, 00150 * <i>C</i><sub>10</sub>, 00151 * <i>C</i><sub>20</sub>, 00152 * <i>C</i><sub>30</sub>, 00153 * <i>C</i><sub>11</sub>, 00154 * <i>C</i><sub>21</sub>, 00155 * <i>C</i><sub>31</sub>, 00156 * <i>C</i><sub>22</sub>, 00157 * <i>C</i><sub>32</sub>, 00158 * <i>C</i><sub>33</sub>. 00159 * In general the (\e n,\e m) element is at index \e m \e N − \e m 00160 * (\e m − 1)/2 + \e n. The layout of \e S is the same except that 00161 * the first column is omitted (since the \e m = 0 terms never contribute 00162 * to the sum) and the 0th element is <i>S</i><sub>11</sub> 00163 * 00164 * The class stores <i>pointers</i> to the first elements of \e C and \e S. 00165 * These arrays should not be altered or destroyed during the lifetime of a 00166 * SphericalHarmonic object. 00167 **********************************************************************/ 00168 SphericalHarmonic(const std::vector<real>& C, 00169 const std::vector<real>& S, 00170 int N, real a, unsigned norm = FULL) 00171 : _a(a) 00172 , _norm(norm) 00173 { _c[0] = SphericalEngine::coeff(C, S, N); } 00174 00175 /** 00176 * Constructor with a subset of coefficients specified. 00177 * 00178 * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>. 00179 * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>. 00180 * @param[in] N the degree used to determine the layout of \e C and \e S. 00181 * @param[in] nmx the maximum degree used in the sum. The sum over \e n is 00182 * from 0 thru \e nmx. 00183 * @param[in] mmx the maximum order used in the sum. The sum over \e m is 00184 * from 0 thru min(\e n, \e mmx). 00185 * @param[in] a the reference radius appearing in the definition of the 00186 * sum. 00187 * @param[in] norm the normalization for the associated Legendre 00188 * polynomials, either SphericalHarmonic::FULL (the default) or 00189 * SphericalHarmonic::SCHMIDT. 00190 * @exception GeographicErr if \e N, \e nmx, and \e mmx do not satisfy 00191 * \e N ≥ \e nmx ≥ \e mmx ≥ −1. 00192 * @exception GeographicErr if \e C or \e S is not big enough to hold the 00193 * coefficients. 00194 * 00195 * The class stores <i>pointers</i> to the first elements of \e C and \e S. 00196 * These arrays should not be altered or destroyed during the lifetime of a 00197 * SphericalHarmonic object. 00198 **********************************************************************/ 00199 SphericalHarmonic(const std::vector<real>& C, 00200 const std::vector<real>& S, 00201 int N, int nmx, int mmx, 00202 real a, unsigned norm = FULL) 00203 : _a(a) 00204 , _norm(norm) 00205 { _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); } 00206 00207 /** 00208 * A default constructor so that the object can be created when the 00209 * constructor for another object is initialized. This default object can 00210 * then be reset with the default copy assignment operator. 00211 **********************************************************************/ 00212 SphericalHarmonic() {} 00213 00214 /** 00215 * Compute the spherical harmonic sum. 00216 * 00217 * @param[in] x cartesian coordinate. 00218 * @param[in] y cartesian coordinate. 00219 * @param[in] z cartesian coordinate. 00220 * @return \e V the spherical harmonic sum. 00221 * 00222 * This routine requires constant memory and thus never throws an 00223 * exception. 00224 **********************************************************************/ 00225 Math::real operator()(real x, real y, real z) const { 00226 real f[] = {1}; 00227 real v = 0; 00228 real dummy; 00229 switch (_norm) { 00230 case FULL: 00231 v = SphericalEngine::Value<false, SphericalEngine::FULL, 1> 00232 (_c, f, x, y, z, _a, dummy, dummy, dummy); 00233 break; 00234 case SCHMIDT: 00235 v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1> 00236 (_c, f, x, y, z, _a, dummy, dummy, dummy); 00237 break; 00238 } 00239 return v; 00240 } 00241 00242 /** 00243 * Compute a spherical harmonic sum and its gradient. 00244 * 00245 * @param[in] x cartesian coordinate. 00246 * @param[in] y cartesian coordinate. 00247 * @param[in] z cartesian coordinate. 00248 * @param[out] gradx \e x component of the gradient 00249 * @param[out] grady \e y component of the gradient 00250 * @param[out] gradz \e z component of the gradient 00251 * @return \e V the spherical harmonic sum. 00252 * 00253 * This is the same as the previous function, except that the components of 00254 * the gradients of the sum in the \e x, \e y, and \e z directions are 00255 * computed. This routine requires constant memory and thus never throws 00256 * an exception. 00257 **********************************************************************/ 00258 Math::real operator()(real x, real y, real z, 00259 real& gradx, real& grady, real& gradz) const { 00260 real f[] = {1}; 00261 real v = 0; 00262 switch (_norm) { 00263 case FULL: 00264 v = SphericalEngine::Value<true, SphericalEngine::FULL, 1> 00265 (_c, f, x, y, z, _a, gradx, grady, gradz); 00266 break; 00267 case SCHMIDT: 00268 v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1> 00269 (_c, f, x, y, z, _a, gradx, grady, gradz); 00270 break; 00271 } 00272 return v; 00273 } 00274 00275 /** 00276 * Create a CircularEngine to allow the efficient evaluation of several 00277 * points on a circle of latitude. 00278 * 00279 * @param[in] p the radius of the circle. 00280 * @param[in] z the height of the circle above the equatorial plane. 00281 * @param[in] gradp if true the returned object will be able to compute the 00282 * gradient of the sum. 00283 * @exception std::bad_alloc if the memory for the CircularEngine can't be 00284 * allocated. 00285 * @return the CircularEngine object. 00286 * 00287 * SphericalHarmonic::operator()() exchanges the order of the sums in the 00288 * definition, i.e., ∑<sub>n = 0..N</sub> ∑<sub>m = 0..n</sub> 00289 * becomes ∑<sub>m = 0..N</sub> ∑<sub>n = m..N</sub>. 00290 * SphericalHarmonic::Circle performs the inner sum over degree \e n (which 00291 * entails about <i>N</i><sup>2</sup> operations). Calling 00292 * CircularEngine::operator()() on the returned object performs the outer 00293 * sum over the order \e m (about \e N operations). 00294 * 00295 * Here's an example of computing the spherical sum at a sequence of 00296 * longitudes without using a CircularEngine object \code 00297 SphericalHarmonic h(...); // Create the SphericalHarmonic object 00298 double r = 2, lat = 33, lon0 = 44, dlon = 0.01; 00299 double 00300 phi = lat * Math::degree<double>(), 00301 z = r * sin(phi), p = r * cos(phi); 00302 for (int i = 0; i <= 100; ++i) { 00303 real 00304 lon = lon0 + i * dlon, 00305 lam = lon * Math::degree<double>(); 00306 std::cout << lon << " " << h(p * cos(lam), p * sin(lam), z) << "\n"; 00307 } 00308 \endcode 00309 * Here is the same calculation done using a CircularEngine object. This 00310 * will be about <i>N</i>/2 times faster. \code 00311 SphericalHarmonic h(...); // Create the SphericalHarmonic object 00312 double r = 2, lat = 33, lon0 = 44, dlon = 0.01; 00313 double 00314 phi = lat * Math::degree<double>(), 00315 z = r * sin(phi), p = r * cos(phi); 00316 CircularEngine c(h(p, z, false)); // Create the CircularEngine object 00317 for (int i = 0; i <= 100; ++i) { 00318 real 00319 lon = lon0 + i * dlon; 00320 std::cout << lon << " " << c(lon) << "\n"; 00321 } 00322 \endcode 00323 **********************************************************************/ 00324 CircularEngine Circle(real p, real z, bool gradp) const { 00325 real f[] = {1}; 00326 switch (_norm) { 00327 case FULL: 00328 return gradp ? 00329 SphericalEngine::Circle<true, SphericalEngine::FULL, 1> 00330 (_c, f, p, z, _a) : 00331 SphericalEngine::Circle<false, SphericalEngine::FULL, 1> 00332 (_c, f, p, z, _a); 00333 break; 00334 case SCHMIDT: 00335 default: // To avoid compiler warnings 00336 return gradp ? 00337 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1> 00338 (_c, f, p, z, _a) : 00339 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1> 00340 (_c, f, p, z, _a); 00341 break; 00342 } 00343 } 00344 00345 /** 00346 * @return the zeroth SphericalEngine::coeff object. 00347 **********************************************************************/ 00348 const SphericalEngine::coeff& Coefficients() const 00349 { return _c[0]; } 00350 }; 00351 00352 } // namespace GeographicLib 00353 00354 #endif // GEOGRAPHICLIB_SPHERICALHARMONIC_HPP