1 /** 2 * \file SphericalHarmonic2.hpp 3 * \brief Header for GeographicLib::SphericalHarmonic2 class 4 * 5 * Copyright (c) Charles Karney (2011-2012) <charles@karney.com> and licensed 6 * under the MIT/X11 License. For more information, see 7 * https://geographiclib.sourceforge.io/ 8 **********************************************************************/ 9 10 #if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP) 11 #define GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP 1 12 13 #include <vector> 14 #include <GeographicLib/Constants.hpp> 15 #include <GeographicLib/SphericalEngine.hpp> 16 #include <GeographicLib/CircularEngine.hpp> 17 18 namespace GeographicLib { 19 20 /** 21 * \brief Spherical harmonic series with two corrections to the coefficients 22 * 23 * This classes is similar to SphericalHarmonic, except that the coefficients 24 * <i>C</i><sub><i>nm</i></sub> are replaced by 25 * <i>C</i><sub><i>nm</i></sub> + \e tau' <i>C'</i><sub><i>nm</i></sub> + \e 26 * tau'' <i>C''</i><sub><i>nm</i></sub> (and similarly for 27 * <i>S</i><sub><i>nm</i></sub>). 28 * 29 * Example of use: 30 * \include example-SphericalHarmonic2.cpp 31 **********************************************************************/ 32 33 // Don't include the GEOGRPAHIC_EXPORT because this header-only class isn't 34 // used by any other classes in the library. 35 class /*GEOGRAPHICLIB_EXPORT*/ SphericalHarmonic2 { 36 public: 37 /** 38 * Supported normalizations for associate Legendre polynomials. 39 **********************************************************************/ 40 enum normalization { 41 /** 42 * Fully normalized associated Legendre polynomials. See 43 * SphericalHarmonic::FULL for documentation. 44 * 45 * @hideinitializer 46 **********************************************************************/ 47 FULL = SphericalEngine::FULL, 48 /** 49 * Schmidt semi-normalized associated Legendre polynomials. See 50 * SphericalHarmonic::SCHMIDT for documentation. 51 * 52 * @hideinitializer 53 **********************************************************************/ 54 SCHMIDT = SphericalEngine::SCHMIDT, 55 }; 56 57 private: 58 typedef Math::real real; 59 SphericalEngine::coeff _c[3]; 60 real _a; 61 unsigned _norm; 62 63 public: 64 /** 65 * Constructor with a full set of coefficients specified. 66 * 67 * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>. 68 * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>. 69 * @param[in] N the maximum degree and order of the sum 70 * @param[in] C1 the coefficients <i>C'</i><sub><i>nm</i></sub>. 71 * @param[in] S1 the coefficients <i>S'</i><sub><i>nm</i></sub>. 72 * @param[in] N1 the maximum degree and order of the first correction 73 * coefficients <i>C'</i><sub><i>nm</i></sub> and 74 * <i>S'</i><sub><i>nm</i></sub>. 75 * @param[in] C2 the coefficients <i>C''</i><sub><i>nm</i></sub>. 76 * @param[in] S2 the coefficients <i>S''</i><sub><i>nm</i></sub>. 77 * @param[in] N2 the maximum degree and order of the second correction 78 * coefficients <i>C'</i><sub><i>nm</i></sub> and 79 * <i>S'</i><sub><i>nm</i></sub>. 80 * @param[in] a the reference radius appearing in the definition of the 81 * sum. 82 * @param[in] norm the normalization for the associated Legendre 83 * polynomials, either SphericalHarmonic2::FULL (the default) or 84 * SphericalHarmonic2::SCHMIDT. 85 * @exception GeographicErr if \e N and \e N1 do not satisfy \e N ≥ 86 * \e N1 ≥ −1, and similarly for \e N2. 87 * @exception GeographicErr if any of the vectors of coefficients is not 88 * large enough. 89 * 90 * See SphericalHarmonic for the way the coefficients should be stored. \e 91 * N1 and \e N2 should satisfy \e N1 ≤ \e N and \e N2 ≤ \e N. 92 * 93 * The class stores <i>pointers</i> to the first elements of \e C, \e S, \e 94 * C', \e S', \e C'', and \e S''. These arrays should not be altered or 95 * destroyed during the lifetime of a SphericalHarmonic object. 96 **********************************************************************/ SphericalHarmonic2(const std::vector<real> & C,const std::vector<real> & S,int N,const std::vector<real> & C1,const std::vector<real> & S1,int N1,const std::vector<real> & C2,const std::vector<real> & S2,int N2,real a,unsigned norm=FULL)97 SphericalHarmonic2(const std::vector<real>& C, 98 const std::vector<real>& S, 99 int N, 100 const std::vector<real>& C1, 101 const std::vector<real>& S1, 102 int N1, 103 const std::vector<real>& C2, 104 const std::vector<real>& S2, 105 int N2, 106 real a, unsigned norm = FULL) 107 : _a(a) 108 , _norm(norm) { 109 if (!(N1 <= N && N2 <= N)) 110 throw GeographicErr("N1 and N2 cannot be larger that N"); 111 _c[0] = SphericalEngine::coeff(C, S, N); 112 _c[1] = SphericalEngine::coeff(C1, S1, N1); 113 _c[2] = SphericalEngine::coeff(C2, S2, N2); 114 } 115 116 /** 117 * Constructor with a subset of coefficients specified. 118 * 119 * @param[in] C the coefficients <i>C</i><sub><i>nm</i></sub>. 120 * @param[in] S the coefficients <i>S</i><sub><i>nm</i></sub>. 121 * @param[in] N the degree used to determine the layout of \e C and \e S. 122 * @param[in] nmx the maximum degree used in the sum. The sum over \e n is 123 * from 0 thru \e nmx. 124 * @param[in] mmx the maximum order used in the sum. The sum over \e m is 125 * from 0 thru min(\e n, \e mmx). 126 * @param[in] C1 the coefficients <i>C'</i><sub><i>nm</i></sub>. 127 * @param[in] S1 the coefficients <i>S'</i><sub><i>nm</i></sub>. 128 * @param[in] N1 the degree used to determine the layout of \e C' and \e 129 * S'. 130 * @param[in] nmx1 the maximum degree used for \e C' and \e S'. 131 * @param[in] mmx1 the maximum order used for \e C' and \e S'. 132 * @param[in] C2 the coefficients <i>C''</i><sub><i>nm</i></sub>. 133 * @param[in] S2 the coefficients <i>S''</i><sub><i>nm</i></sub>. 134 * @param[in] N2 the degree used to determine the layout of \e C'' and \e 135 * S''. 136 * @param[in] nmx2 the maximum degree used for \e C'' and \e S''. 137 * @param[in] mmx2 the maximum order used for \e C'' and \e S''. 138 * @param[in] a the reference radius appearing in the definition of the 139 * sum. 140 * @param[in] norm the normalization for the associated Legendre 141 * polynomials, either SphericalHarmonic2::FULL (the default) or 142 * SphericalHarmonic2::SCHMIDT. 143 * @exception GeographicErr if the parameters do not satisfy \e N ≥ \e 144 * nmx ≥ \e mmx ≥ −1; \e N1 ≥ \e nmx1 ≥ \e mmx1 ≥ 145 * −1; \e N ≥ \e N1; \e nmx ≥ \e nmx1; \e mmx ≥ \e mmx1; 146 * and similarly for \e N2, \e nmx2, and \e mmx2. 147 * @exception GeographicErr if any of the vectors of coefficients is not 148 * large enough. 149 * 150 * The class stores <i>pointers</i> to the first elements of \e C, \e S, \e 151 * C', \e S', \e C'', and \e S''. These arrays should not be altered or 152 * destroyed during the lifetime of a SphericalHarmonic object. 153 **********************************************************************/ SphericalHarmonic2(const std::vector<real> & C,const std::vector<real> & S,int N,int nmx,int mmx,const std::vector<real> & C1,const std::vector<real> & S1,int N1,int nmx1,int mmx1,const std::vector<real> & C2,const std::vector<real> & S2,int N2,int nmx2,int mmx2,real a,unsigned norm=FULL)154 SphericalHarmonic2(const std::vector<real>& C, 155 const std::vector<real>& S, 156 int N, int nmx, int mmx, 157 const std::vector<real>& C1, 158 const std::vector<real>& S1, 159 int N1, int nmx1, int mmx1, 160 const std::vector<real>& C2, 161 const std::vector<real>& S2, 162 int N2, int nmx2, int mmx2, 163 real a, unsigned norm = FULL) 164 : _a(a) 165 , _norm(norm) { 166 if (!(nmx1 <= nmx && nmx2 <= nmx)) 167 throw GeographicErr("nmx1 and nmx2 cannot be larger that nmx"); 168 if (!(mmx1 <= mmx && mmx2 <= mmx)) 169 throw GeographicErr("mmx1 and mmx2 cannot be larger that mmx"); 170 _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); 171 _c[1] = SphericalEngine::coeff(C1, S1, N1, nmx1, mmx1); 172 _c[2] = SphericalEngine::coeff(C2, S2, N2, nmx2, mmx2); 173 } 174 175 /** 176 * A default constructor so that the object can be created when the 177 * constructor for another object is initialized. This default object can 178 * then be reset with the default copy assignment operator. 179 **********************************************************************/ SphericalHarmonic2()180 SphericalHarmonic2() {} 181 182 /** 183 * Compute a spherical harmonic sum with two correction terms. 184 * 185 * @param[in] tau1 multiplier for correction coefficients \e C' and \e S'. 186 * @param[in] tau2 multiplier for correction coefficients \e C'' and \e 187 * S''. 188 * @param[in] x cartesian coordinate. 189 * @param[in] y cartesian coordinate. 190 * @param[in] z cartesian coordinate. 191 * @return \e V the spherical harmonic sum. 192 * 193 * This routine requires constant memory and thus never throws an 194 * exception. 195 **********************************************************************/ operator ()(real tau1,real tau2,real x,real y,real z) const196 Math::real operator()(real tau1, real tau2, real x, real y, real z) 197 const { 198 real f[] = {1, tau1, tau2}; 199 real v = 0; 200 real dummy; 201 switch (_norm) { 202 case FULL: 203 v = SphericalEngine::Value<false, SphericalEngine::FULL, 3> 204 (_c, f, x, y, z, _a, dummy, dummy, dummy); 205 break; 206 case SCHMIDT: 207 default: // To avoid compiler warnings 208 v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3> 209 (_c, f, x, y, z, _a, dummy, dummy, dummy); 210 break; 211 } 212 return v; 213 } 214 215 /** 216 * Compute a spherical harmonic sum with two correction terms and its 217 * gradient. 218 * 219 * @param[in] tau1 multiplier for correction coefficients \e C' and \e S'. 220 * @param[in] tau2 multiplier for correction coefficients \e C'' and \e 221 * S''. 222 * @param[in] x cartesian coordinate. 223 * @param[in] y cartesian coordinate. 224 * @param[in] z cartesian coordinate. 225 * @param[out] gradx \e x component of the gradient 226 * @param[out] grady \e y component of the gradient 227 * @param[out] gradz \e z component of the gradient 228 * @return \e V the spherical harmonic sum. 229 * 230 * This is the same as the previous function, except that the components of 231 * the gradients of the sum in the \e x, \e y, and \e z directions are 232 * computed. This routine requires constant memory and thus never throws 233 * an exception. 234 **********************************************************************/ operator ()(real tau1,real tau2,real x,real y,real z,real & gradx,real & grady,real & gradz) const235 Math::real operator()(real tau1, real tau2, real x, real y, real z, 236 real& gradx, real& grady, real& gradz) const { 237 real f[] = {1, tau1, tau2}; 238 real v = 0; 239 switch (_norm) { 240 case FULL: 241 v = SphericalEngine::Value<true, SphericalEngine::FULL, 3> 242 (_c, f, x, y, z, _a, gradx, grady, gradz); 243 break; 244 case SCHMIDT: 245 default: // To avoid compiler warnings 246 v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3> 247 (_c, f, x, y, z, _a, gradx, grady, gradz); 248 break; 249 } 250 return v; 251 } 252 253 /** 254 * Create a CircularEngine to allow the efficient evaluation of several 255 * points on a circle of latitude at fixed values of \e tau1 and \e tau2. 256 * 257 * @param[in] tau1 multiplier for correction coefficients \e C' and \e S'. 258 * @param[in] tau2 multiplier for correction coefficients \e C'' and \e 259 * S''. 260 * @param[in] p the radius of the circle. 261 * @param[in] z the height of the circle above the equatorial plane. 262 * @param[in] gradp if true the returned object will be able to compute the 263 * gradient of the sum. 264 * @exception std::bad_alloc if the memory for the CircularEngine can't be 265 * allocated. 266 * @return the CircularEngine object. 267 * 268 * SphericalHarmonic2::operator()() exchanges the order of the sums in the 269 * definition, i.e., ∑<sub><i>n</i> = 0..<i>N</i></sub> 270 * ∑<sub><i>m</i> = 0..<i>n</i></sub> becomes ∑<sub><i>m</i> = 271 * 0..<i>N</i></sub> ∑<sub><i>n</i> = <i>m</i>..<i>N</i></sub>.. 272 * SphericalHarmonic2::Circle performs the inner sum over degree \e n 273 * (which entails about <i>N</i><sup>2</sup> operations). Calling 274 * CircularEngine::operator()() on the returned object performs the outer 275 * sum over the order \e m (about \e N operations). 276 * 277 * See SphericalHarmonic::Circle for an example of its use. 278 **********************************************************************/ Circle(real tau1,real tau2,real p,real z,bool gradp) const279 CircularEngine Circle(real tau1, real tau2, real p, real z, bool gradp) 280 const { 281 real f[] = {1, tau1, tau2}; 282 switch (_norm) { 283 case FULL: 284 return gradp ? 285 SphericalEngine::Circle<true, SphericalEngine::FULL, 3> 286 (_c, f, p, z, _a) : 287 SphericalEngine::Circle<false, SphericalEngine::FULL, 3> 288 (_c, f, p, z, _a); 289 break; 290 case SCHMIDT: 291 default: // To avoid compiler warnings 292 return gradp ? 293 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3> 294 (_c, f, p, z, _a) : 295 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3> 296 (_c, f, p, z, _a); 297 break; 298 } 299 } 300 301 /** 302 * @return the zeroth SphericalEngine::coeff object. 303 **********************************************************************/ Coefficients() const304 const SphericalEngine::coeff& Coefficients() const 305 { return _c[0]; } 306 /** 307 * @return the first SphericalEngine::coeff object. 308 **********************************************************************/ Coefficients1() const309 const SphericalEngine::coeff& Coefficients1() const 310 { return _c[1]; } 311 /** 312 * @return the second SphericalEngine::coeff object. 313 **********************************************************************/ Coefficients2() const314 const SphericalEngine::coeff& Coefficients2() const 315 { return _c[2]; } 316 }; 317 318 } // namespace GeographicLib 319 320 #endif // GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP 321