GeographicLib  1.21
CircularEngine.hpp
Go to the documentation of this file.
00001 /**
00002  * \file CircularEngine.hpp
00003  * \brief Header for GeographicLib::CircularEngine 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_CIRCULARENGINE_HPP)
00011 #define GEOGRAPHICLIB_CIRCULARENGINE_HPP \
00012   "$Id: d0528f468369dbc1b7af11e02278ad7a361d398b $"
00013 
00014 #include <vector>
00015 #include <GeographicLib/Constants.hpp>
00016 #include <GeographicLib/SphericalEngine.hpp>
00017 
00018 #if defined(_MSC_VER)
00019 // Squelch warnings about dll vs vector
00020 #pragma warning (push)
00021 #pragma warning (disable: 4251)
00022 #endif
00023 
00024 namespace GeographicLib {
00025 
00026   /**
00027    * \brief Spherical Harmonic sums for a circle.
00028    *
00029    * The class is a companion to SphericalEngine.  If the results of a
00030    * spherical harmonic sum are needed for several points on a circle of
00031    * constant latitude \e lat and height \e h, then SphericalEngine::Circle can
00032    * compute the inner sum, which is independent of longitude \e lon, and
00033    * produce a CircularEngine object.  CircularEngine::operator()() can
00034    * then be used to perform the outer sum for particular vales of \e lon.
00035    * This can lead to substantial improvements in computational speed for high
00036    * degree sum (approximately by a factor of \e N / 2 where \e N is the
00037    * maximum degree).
00038    *
00039    * CircularEngine is tightly linked to the internals of SphericalEngine.  For
00040    * that reason, the constructor for this class is private.  Use
00041    * SphericalHarmonic::Circle, SphericalHarmonic1::Circle, and
00042    * SphericalHarmonic2::Circle to create instances of this class.
00043    *
00044    * CircularEngine stores the coefficients needed to allow the summation over
00045    * order to be performed in 2 or 6 vectors of length \e M + 1 (depending on
00046    * whether gradients are to be calculated).  For this reason the constructor
00047    * may throw a bad_alloc exception.
00048    *
00049    * Example of use:
00050    * \include example-CircularEngine.cpp
00051    **********************************************************************/
00052 
00053   class GEOGRAPHIC_EXPORT CircularEngine {
00054   private:
00055     typedef Math::real real;
00056     enum normalization {
00057       FULL = SphericalEngine::FULL,
00058       SCHMIDT = SphericalEngine::SCHMIDT,
00059     };
00060     int _M;
00061     bool _gradp;
00062     unsigned _norm;
00063     real _a, _r, _u, _t;
00064     std::vector<real> _wc, _ws, _wrc, _wrs, _wtc, _wts;
00065     real _q, _uq, _uq2;
00066 
00067     Math::real Value(bool gradp, real cl, real sl,
00068                      real& gradx, real& grady, real& gradz) const throw();
00069 
00070     static inline void cossin(real x, real& cosx, real& sinx) {
00071       x = x >= 180 ? x - 360 : (x < -180 ? x + 360 : x);
00072       real xi = x * Math::degree<real>();
00073       cosx = std::abs(x) ==   90 ? 0 : cos(xi);
00074       sinx =          x  == -180 ? 0 : sin(xi);
00075     }
00076 
00077     friend class SphericalEngine;
00078     friend class GravityCircle;  // Access to cossin
00079     friend class MagneticCircle; // Access to cossin
00080     CircularEngine(int M, bool gradp, unsigned norm,
00081                    real a, real r, real u, real t)
00082       : _M(M)
00083       , _gradp(gradp)
00084       , _norm(norm)
00085       , _a(a)
00086       , _r(r)
00087       , _u(u)
00088       , _t(t)
00089       , _wc(std::vector<real>(_M + 1, 0))
00090       , _ws(std::vector<real>(_M + 1, 0))
00091       , _wrc(std::vector<real>(_gradp ? _M + 1 : 0, 0))
00092       , _wrs(std::vector<real>(_gradp ? _M + 1 : 0, 0))
00093       , _wtc(std::vector<real>(_gradp ? _M + 1 : 0, 0))
00094       , _wts(std::vector<real>(_gradp ? _M + 1 : 0, 0))
00095       {
00096         _q = _a / _r;
00097         _uq = _u * _q;
00098         _uq2 = Math::sq(_uq);
00099       }
00100 
00101     void SetCoeff(int m, real wc, real ws)
00102     { _wc[m] = wc; _ws[m] = ws; }
00103 
00104     void SetCoeff(int m, real wc, real ws,
00105                   real wrc, real wrs, real wtc, real wts) {
00106       _wc[m] = wc; _ws[m] = ws;
00107       if (_gradp) {
00108         _wrc[m] = wrc; _wrs[m] = wrs;
00109         _wtc[m] = wtc; _wts[m] = wts;
00110       }
00111     }
00112 
00113   public:
00114 
00115     /**
00116      * A default constructor.  CircularEngine::operator()() on the resulting
00117      * object returns zero.  The resulting object can be assigned to the result
00118      * of SphericalHarmonic::Circle.
00119      **********************************************************************/
00120     CircularEngine()
00121       : _M(-1)
00122       , _gradp(true)
00123       , _u(0)
00124       , _t(1)
00125       {}
00126 
00127     /**
00128      * Evaluate the sum for a particular longitude given in terms of its
00129      * cosine and sine.
00130      *
00131      * @param[in] coslon the cosine of the longitude.
00132      * @param[in] sinlon the sine of the longitude.
00133      * @return \e V the value of the sum.
00134      *
00135      * The arguments must satisfy <i>coslon</i><sup>2</sup> +
00136      * <i>sinlon</i><sup>2</sup> = 1.
00137      **********************************************************************/
00138     Math::real operator()(real coslon, real sinlon) const throw() {
00139       real dummy;
00140       return Value(false, coslon, sinlon, dummy, dummy, dummy);
00141     }
00142 
00143     /**
00144      * Evaluate the sum for a particular longitude.
00145      *
00146      * @param[in] lon the longitude (degrees).
00147      * @return \e V the value of the sum.
00148      **********************************************************************/
00149     Math::real operator()(real lon) const throw() {
00150       real coslon, sinlon;
00151       cossin(lon, coslon, sinlon);
00152       return (*this)(coslon, sinlon);
00153     }
00154 
00155     /**
00156      * Evaluate the sum and its gradient for a particular longitude given in
00157      * terms of its cosine and sine.
00158      *
00159      * @param[in] coslon the cosine of the longitude.
00160      * @param[in] sinlon the sine of the longitude.
00161      * @param[out] gradx \e x component of the gradient.
00162      * @param[out] grady \e y component of the gradient.
00163      * @param[out] gradz \e z component of the gradient.
00164      * @return \e V the value of the sum.
00165      *
00166      * The gradients will only be computed if the CircularEngine object was
00167      * created with this capability (e.g., via \e gradp = true in
00168      * SphericalHarmonic::Circle).  If not, \e gradx, etc., will not be
00169      * touched.  The arguments must satisfy <i>coslon</i><sup>2</sup> +
00170      * <i>sinlon</i><sup>2</sup> = 1.
00171      **********************************************************************/
00172     Math::real operator()(real coslon, real sinlon,
00173                           real& gradx, real& grady, real& gradz) const throw() {
00174       return Value(true, coslon, sinlon, gradx, grady, gradz);
00175     }
00176 
00177     /**
00178      * Evaluate the sum and its gradient for a particular longitude.
00179      *
00180      * @param[in] lon the longitude (degrees).
00181      * @param[out] gradx \e x component of the gradient.
00182      * @param[out] grady \e y component of the gradient.
00183      * @param[out] gradz \e z component of the gradient.
00184      * @return \e V the value of the sum.
00185      *
00186      * The gradients will only be computed if the CircularEngine object was
00187      * created with this capability (e.g., via \e gradp = true in
00188      * SphericalHarmonic::Circle).  If not, \e gradx, etc., will not be
00189      * touched.
00190      **********************************************************************/
00191     Math::real operator()(real lon,
00192                           real& gradx, real& grady, real& gradz) const throw() {
00193       real coslon, sinlon;
00194       cossin(lon, coslon, sinlon);
00195       return (*this)(coslon, sinlon, gradx, grady, gradz);
00196     }
00197   };
00198 
00199 } // namespace GeographicLib
00200 
00201 #if defined(_MSC_VER)
00202 #pragma warning (pop)
00203 #endif
00204 
00205 #endif  // GEOGRAPHICLIB_CIRCULARENGINE_HPP