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 1 00012 00013 #include <vector> 00014 #include <GeographicLib/Constants.hpp> 00015 #include <GeographicLib/SphericalEngine.hpp> 00016 00017 #if defined(_MSC_VER) 00018 // Squelch warnings about dll vs vector 00019 # pragma warning (push) 00020 # pragma warning (disable: 4251) 00021 #endif 00022 00023 namespace GeographicLib { 00024 00025 /** 00026 * \brief Spherical harmonic sums for a circle 00027 * 00028 * The class is a companion to SphericalEngine. If the results of a 00029 * spherical harmonic sum are needed for several points on a circle of 00030 * constant latitude \e lat and height \e h, then SphericalEngine::Circle can 00031 * compute the inner sum, which is independent of longitude \e lon, and 00032 * produce a CircularEngine object. CircularEngine::operator()() can 00033 * then be used to perform the outer sum for particular vales of \e lon. 00034 * This can lead to substantial improvements in computational speed for high 00035 * degree sum (approximately by a factor of \e N / 2 where \e N is the 00036 * maximum degree). 00037 * 00038 * CircularEngine is tightly linked to the internals of SphericalEngine. For 00039 * that reason, the constructor for this class is private. Use 00040 * SphericalHarmonic::Circle, SphericalHarmonic1::Circle, and 00041 * SphericalHarmonic2::Circle to create instances of this class. 00042 * 00043 * CircularEngine stores the coefficients needed to allow the summation over 00044 * order to be performed in 2 or 6 vectors of length \e M + 1 (depending on 00045 * whether gradients are to be calculated). For this reason the constructor 00046 * may throw a std::bad_alloc exception. 00047 * 00048 * Example of use: 00049 * \include example-CircularEngine.cpp 00050 **********************************************************************/ 00051 00052 class GEOGRAPHICLIB_EXPORT CircularEngine { 00053 private: 00054 typedef Math::real real; 00055 enum normalization { 00056 FULL = SphericalEngine::FULL, 00057 SCHMIDT = SphericalEngine::SCHMIDT, 00058 }; 00059 int _M; 00060 bool _gradp; 00061 unsigned _norm; 00062 real _a, _r, _u, _t; 00063 std::vector<real> _wc, _ws, _wrc, _wrs, _wtc, _wts; 00064 real _q, _uq, _uq2; 00065 00066 Math::real Value(bool gradp, real cl, real sl, 00067 real& gradx, real& grady, real& gradz) const; 00068 00069 static inline void cossin(real x, real& cosx, real& sinx) { 00070 using std::abs; using std::cos; using std::sin; 00071 x = x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); 00072 real xi = x * Math::degree(); 00073 cosx = 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 { 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 { 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 { 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 { 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