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