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