GeographicLib  1.21
SphericalHarmonic2.hpp
Go to the documentation of this file.
00001 /**
00002  * \file SphericalHarmonic2.hpp
00003  * \brief Header for GeographicLib::SphericalHarmonic2 class
00004  *
00005  * Copyright (c) Charles Karney (2011, 2012) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP)
00011 #define GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP \
00012   "$Id: ce4cda614c1966dea65610bc73bc4db562677fa8 $"
00013 
00014 #include <vector>
00015 #include <GeographicLib/Constants.hpp>
00016 #include <GeographicLib/SphericalEngine.hpp>
00017 #include <GeographicLib/CircularEngine.hpp>
00018 
00019 namespace GeographicLib {
00020 
00021   /**
00022    * \brief Spherical Harmonic series with two corrections to the coefficients.
00023    *
00024    * This classes is similar to SphericalHarmonic, except that the coefficients
00025    * \e C<sub>\e nm</sub> are replaced by \e C<sub>\e nm</sub> + \e tau'
00026    * C'<sub>\e nm</sub> + \e tau'' C''<sub>\e nm</sub> (and similarly for \e
00027    * S<sub>\e nm</sub>).
00028    *
00029    * Example of use:
00030    * \include example-SphericalHarmonic2.cpp
00031    **********************************************************************/
00032 
00033   class GEOGRAPHIC_EXPORT SphericalHarmonic2 {
00034   public:
00035     /**
00036      * Supported normalizations for associate Legendre polynomials.
00037      **********************************************************************/
00038     enum normalization {
00039       /**
00040        * Fully normalized associated Legendre polynomials.  See
00041        * SphericalHarmonic::FULL for documentation.
00042        *
00043        * @hideinitializer
00044        **********************************************************************/
00045       FULL = SphericalEngine::FULL,
00046       /**
00047        * Schmidt semi-normalized associated Legendre polynomials.  See
00048        * SphericalHarmonic::SCHMIDT for documentation.
00049        *
00050        * @hideinitializer
00051        **********************************************************************/
00052       SCHMIDT = SphericalEngine::SCHMIDT,
00053       /// \cond SKIP
00054       // These are deprecated...
00055       full = FULL,
00056       schmidt = SCHMIDT,
00057       /// \endcond
00058     };
00059 
00060   private:
00061     typedef Math::real real;
00062     SphericalEngine::coeff _c[3];
00063     real _a;
00064     unsigned _norm;
00065 
00066   public:
00067     /**
00068      * Constructor with a full set of coefficients specified.
00069      *
00070      * @param[in] C the coefficients \e C<sub>\e nm</sub>.
00071      * @param[in] S the coefficients \e S<sub>\e nm</sub>.
00072      * @param[in] N the maximum degree and order of the sum
00073      * @param[in] C1 the coefficients \e C'<sub>\e nm</sub>.
00074      * @param[in] S1 the coefficients \e S'<sub>\e nm</sub>.
00075      * @param[in] N1 the maximum degree and order of the first correction
00076      *   coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>.
00077      * @param[in] C2 the coefficients \e C''<sub>\e nm</sub>.
00078      * @param[in] S2 the coefficients \e S''<sub>\e nm</sub>.
00079      * @param[in] N2 the maximum degree and order of the second correction
00080      *   coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>.
00081      * @param[in] a the reference radius appearing in the definition of the
00082      *   sum.
00083      * @param[in] norm the normalization for the associated Legendre
00084      *   polynomials, either SphericalHarmonic2::FULL (the default) or
00085      *   SphericalHarmonic2::SCHMIDT.
00086      *
00087      * See SphericalHarmonic for the way the coefficients should be stored.  \e
00088      * N1 and \e N2 should satisfy \e N1 <= \e N and \e N2 <= \e N.
00089      *
00090      * The class stores <i>pointers</i> to the first elements of \e C, \e S, \e
00091      * C', \e S', \e C'', and \e S''.  These arrays should not be altered or
00092      * destroyed during the lifetime of a SphericalHarmonic object.
00093      **********************************************************************/
00094     SphericalHarmonic2(const std::vector<real>& C,
00095                        const std::vector<real>& S,
00096                        int N,
00097                        const std::vector<real>& C1,
00098                        const std::vector<real>& S1,
00099                        int N1,
00100                        const std::vector<real>& C2,
00101                        const std::vector<real>& S2,
00102                        int N2,
00103                        real a, unsigned norm = FULL)
00104       : _a(a)
00105       , _norm(norm) {
00106       if (!(N1 <= N && N2 <= N))
00107         throw GeographicErr("N1 and N2 cannot be larger that N");
00108       _c[0] = SphericalEngine::coeff(C, S, N);
00109       _c[1] = SphericalEngine::coeff(C1, S1, N1);
00110       _c[2] = SphericalEngine::coeff(C2, S2, N2);
00111     }
00112 
00113     /**
00114      * Constructor with a subset of coefficients specified.
00115      *
00116      * @param[in] C the coefficients \e C<sub>\e nm</sub>.
00117      * @param[in] S the coefficients \e S<sub>\e nm</sub>.
00118      * @param[in] N the degree used to determine the layout of \e C and \e S.
00119      * @param[in] nmx the maximum degree used in the sum.  The sum over \e n is
00120      *   from 0 thru \e nmx.
00121      * @param[in] mmx the maximum order used in the sum.  The sum over \e m is
00122      *   from 0 thru min(\e n, \e mmx).
00123      * @param[in] C1 the coefficients \e C'<sub>\e nm</sub>.
00124      * @param[in] S1 the coefficients \e S'<sub>\e nm</sub>.
00125      * @param[in] N1 the degree used to determine the layout of \e C' and \e
00126      *   S'.
00127      * @param[in] nmx1 the maximum degree used for \e C' and \e S'.
00128      * @param[in] mmx1 the maximum order used for \e C' and \e S'.
00129      * @param[in] C2 the coefficients \e C''<sub>\e nm</sub>.
00130      * @param[in] S2 the coefficients \e S''<sub>\e nm</sub>.
00131      * @param[in] N2 the degree used to determine the layout of \e C'' and \e
00132      *   S''.
00133      * @param[in] nmx2 the maximum degree used for \e C'' and \e S''.
00134      * @param[in] mmx2 the maximum order used for \e C'' and \e S''.
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 SphericalHarmonic2::FULL (the default) or
00139      *   SphericalHarmonic2::SCHMIDT.
00140      *
00141      * The class stores <i>pointers</i> to the first elements of \e C, \e S, \e
00142      * C', \e S', \e C'', and \e S''.  These arrays should not be altered or
00143      * destroyed during the lifetime of a SphericalHarmonic object.
00144      **********************************************************************/
00145     SphericalHarmonic2(const std::vector<real>& C,
00146                        const std::vector<real>& S,
00147                        int N, int nmx, int mmx,
00148                        const std::vector<real>& C1,
00149                        const std::vector<real>& S1,
00150                        int N1, int nmx1, int mmx1,
00151                        const std::vector<real>& C2,
00152                        const std::vector<real>& S2,
00153                        int N2, int nmx2, int mmx2,
00154                        real a, unsigned norm = FULL)
00155       : _a(a)
00156       , _norm(norm) {
00157       if (!(nmx1 <= nmx && nmx2 <= nmx))
00158         throw GeographicErr("nmx1 and nmx2 cannot be larger that nmx");
00159       if (!(mmx1 <= mmx && mmx2 <= mmx))
00160         throw GeographicErr("mmx1 and mmx2cannot be larger that mmx");
00161       _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx);
00162       _c[1] = SphericalEngine::coeff(C1, S1, N1, nmx1, mmx1);
00163       _c[2] = SphericalEngine::coeff(C2, S2, N2, nmx2, mmx2);
00164     }
00165 
00166     /**
00167      * A default constructor so that the object can be created when the
00168      * constructor for another object is initialized.  This default object can
00169      * then be reset with the default copy assignment operator.
00170      **********************************************************************/
00171     SphericalHarmonic2() {}
00172 
00173     /**
00174      * Compute a spherical harmonic sum with two correction terms.
00175      *
00176      * @param[in] tau1 multiplier for correction coefficients \e C' and \e S'.
00177      * @param[in] tau2 multiplier for correction coefficients \e C'' and \e S''.
00178      * @param[in] x cartesian coordinate.
00179      * @param[in] y cartesian coordinate.
00180      * @param[in] z cartesian coordinate.
00181      * @return \e V the spherical harmonic sum.
00182      *
00183      * This routine requires constant memory and thus never throws an
00184      * exception.
00185      **********************************************************************/
00186     Math::real operator()(real tau1, real tau2, real x, real y, real z)
00187       const throw() {
00188       real f[] = {1, tau1, tau2};
00189       real v = 0;
00190       real dummy;
00191       switch (_norm) {
00192       case FULL:
00193         v = SphericalEngine::Value<false, SphericalEngine::FULL, 3>
00194           (_c, f, x, y, z, _a, dummy, dummy, dummy);
00195         break;
00196       case SCHMIDT:
00197         v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
00198           (_c, f, x, y, z, _a, dummy, dummy, dummy);
00199         break;
00200       }
00201       return v;
00202     }
00203 
00204     /**
00205      * Compute a spherical harmonic sum with two correction terms and its
00206      * gradient.
00207      *
00208      * @param[in] tau1 multiplier for correction coefficients \e C' and \e S'.
00209      * @param[in] tau2 multiplier for correction coefficients \e C'' and \e S''.
00210      * @param[in] x cartesian coordinate.
00211      * @param[in] y cartesian coordinate.
00212      * @param[in] z cartesian coordinate.
00213      * @param[out] gradx \e x component of the gradient
00214      * @param[out] grady \e y component of the gradient
00215      * @param[out] gradz \e z component of the gradient
00216      * @return \e V the spherical harmonic sum.
00217      *
00218      * This is the same as the previous function, except that the components of
00219      * the gradients of the sum in the \e x, \e y, and \e z directions are
00220      * computed.  This routine requires constant memory and thus never throws
00221      * an exception.
00222      **********************************************************************/
00223     Math::real operator()(real tau1, real tau2, real x, real y, real z,
00224                           real& gradx, real& grady, real& gradz) const throw() {
00225       real f[] = {1, tau1, tau2};
00226       real v = 0;
00227       switch (_norm) {
00228       case FULL:
00229         v = SphericalEngine::Value<true, SphericalEngine::FULL, 3>
00230           (_c, f, x, y, z, _a, gradx, grady, gradz);
00231         break;
00232       case SCHMIDT:
00233         v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
00234           (_c, f, x, y, z, _a, gradx, grady, gradz);
00235         break;
00236       }
00237       return v;
00238     }
00239 
00240     /**
00241      * Create a CircularEngine to allow the efficient evaluation of several
00242      * points on a circle of latitude at fixed values of \e tau1 and \e tau2.
00243      *
00244      * @param[in] tau1 multiplier for correction coefficients \e C' and \e S'.
00245      * @param[in] tau2 multiplier for correction coefficients \e C'' and \e S''.
00246      * @param[in] p the radius of the circle.
00247      * @param[in] z the height of the circle above the equatorial plane.
00248      * @param[in] gradp if true the returned object will be able to compute the
00249      *   gradient of the sum.
00250      * @return the CircularEngine object.
00251      *
00252      * SphericalHarmonic2::operator()() exchanges the order of the sums in the
00253      * definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m =
00254      * 0..N)[sum(n = m..N)[...]].  SphericalHarmonic2::Circle performs the
00255      * inner sum over degree \e n (which entails about <i>N</i><sup>2</sup>
00256      * operations).  Calling CircularEngine::operator()() on the returned
00257      * object performs the outer sum over the order \e m (about \e N
00258      * operations).  This routine may throw a bad_alloc exception in the
00259      * CircularEngine constructor.
00260      *
00261      * See SphericalHarmonic::Circle for an example of its use.
00262      **********************************************************************/
00263     CircularEngine Circle(real tau1, real tau2, real p, real z, bool gradp)
00264       const {
00265       real f[] = {1, tau1, tau2};
00266       switch (_norm) {
00267       case FULL:
00268         return gradp ?
00269           SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
00270           (_c, f, p, z, _a) :
00271           SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
00272           (_c, f, p, z, _a);
00273         break;
00274       case SCHMIDT:
00275       default:                  // To avoid compiler warnings
00276         return gradp ?
00277           SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
00278           (_c, f, p, z, _a) :
00279           SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
00280           (_c, f, p, z, _a);
00281         break;
00282       }
00283     }
00284 
00285     /**
00286      * @return the zeroth SphericalEngine::coeff object.
00287      **********************************************************************/
00288     const SphericalEngine::coeff& Coefficients() const throw()
00289     { return _c[0]; }
00290     /**
00291      * @return the first SphericalEngine::coeff object.
00292      **********************************************************************/
00293     const SphericalEngine::coeff& Coefficients1() const throw()
00294     { return _c[1]; }
00295     /**
00296      * @return the second SphericalEngine::coeff object.
00297      **********************************************************************/
00298     const SphericalEngine::coeff& Coefficients2() const throw()
00299     { return _c[2]; }
00300   };
00301 
00302 } // namespace GeographicLib
00303 
00304 #endif  // GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP