GeographicLib  1.21
SphericalHarmonic.hpp
Go to the documentation of this file.
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