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