1 // 2 // SHCalc.hpp 3 // pb_solvers_code 4 // 5 /* 6 Copyright (c) 2015, Teresa Head-Gordon, Lisa Felberg, Enghui Yap, David Brookes 7 All rights reserved. 8 9 Redistribution and use in source and binary forms, with or without 10 modification, are permitted provided that the following conditions are met: 11 * Redistributions of source code must retain the above copyright 12 notice, this list of conditions and the following disclaimer. 13 * Redistributions in binary form must reproduce the above copyright 14 notice, this list of conditions and the following disclaimer in the 15 documentation and/or other materials provided with the distribution. 16 * Neither the name of UC Berkeley nor the 17 names of its contributors may be used to endorse or promote products 18 derived from this software without specific prior written permission. 19 20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 21 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 22 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 23 DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY 24 DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 25 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 26 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 27 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 28 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 29 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 30 */ 31 32 #ifndef SHCalc_h 33 #define SHCalc_h 34 35 #include <complex> 36 #include <vector> 37 #include <assert.h> 38 #include <memory> 39 #include "BesselCalc.h" 40 #include <unordered_map> 41 42 43 using namespace std; 44 45 46 /* 47 Class for storing the constants that may be used in 48 multiple spherical harmonics calculations 49 */ 50 class SHCalcConstants 51 { 52 protected: 53 int numVals_; // number of poles 54 MyMatrix<double> legConsts1_; // (2l-1)/(l-m) for use in legdre comp 55 MyMatrix<double> legConsts2_; // (l+m-1)/(l-m) for use in legdre compn 56 MyMatrix<double> shConsts_; // sqrt((n-m)!/(n+m)!) in EQ1, Lotan 2006 57 vector<double> dubFac_; // (2l-1)!! double factrl, in legendre recursion 58 59 public: 60 SHCalcConstants(const int num_vals=Constants::MAX_NUM_POLES); 61 get_leg_consts1_val(const int n,const int m)62 const double get_leg_consts1_val(const int n, const int m) 63 { return legConsts1_(n, m); } get_leg_consts2_val(const int n,const int m)64 const double get_leg_consts2_val(const int n, const int m) 65 { return legConsts2_(n, m); } get_sh_consts_val(const int n,const int m)66 const double get_sh_consts_val(const int n, const int m) 67 { return shConsts_(n, m); } get_dub_fac_val(const int i)68 const double get_dub_fac_val(const int i) const 69 { return dubFac_[i]; } get_n()70 const int get_n() const 71 { return numVals_; } 72 }; 73 74 75 /* 76 Class for computing spherical harmonics. This includes 77 calcualtion of the associated Legendre Polynomials 78 79 The spherical harmonics in this case are defined by the equation: 80 Y_(n,m)(theta, phi) = (-1)^m * sqrt((n-m)! / (n + m)!) * 81 P_(n,m)(cos(theta)) * exp(i*m*phi) 82 where P_(n, m) are the associated Legendre polynomials. 83 84 These are constructed dynamically and returned as a matrix of 85 values for every n,m 86 */ 87 class SHCalc 88 { 89 protected: 90 91 int numVals_; //# of poles (output matrix will be NxN) 92 shared_ptr<SHCalcConstants> _consts_; 93 MyMatrix<double> P_; // legendre polynomials 94 MyMatrix<cmplx> Y_; // spherical harmonics calcualted by this class 95 96 void calc_legendre(const double theta); // calculate the legendre polynomial 97 // at every n, m (store in this.P_) 98 99 public: SHCalc()100 SHCalc() {} 101 102 SHCalc(const int num_vals, shared_ptr<SHCalcConstants> _consts); 103 104 // calculate the spherical harmonics at every n, m (store in this.Y_) 105 void calc_sh(const double theta, const double phi); 106 107 // retrieve the result for n, m values 108 cmplx get_result(const int n, const int m); 109 110 // retrieve the full calculated Y_ matrix get_full_result()111 MyMatrix<cmplx> get_full_result() { return Y_; } 112 113 double get_legendre_result( int n, int m ); get_num_vals()114 int get_num_vals() { return numVals_; } 115 }; 116 117 118 // class for pre-calculated spherical harmonics. Map 119 class PreCalcSH 120 { 121 protected: 122 123 unordered_map<Pt, MyMatrix<cmplx> > map_; 124 125 public: 126 PreCalcSH()127 PreCalcSH() { } 128 129 // calculates spherical harmonics at this point and add to hashmap 130 void calc_and_add(Pt p, shared_ptr<SHCalc> sh); 131 132 cmplx get_sh(Pt p, int n, int m); 133 clear_sh()134 void clear_sh() 135 { 136 map_.clear(); 137 } 138 // MyMatrix<cmplx> get_full_sh(Pt p) { return map_[p]; } 139 140 // bool is_calculated(Pt p) { return map_.find(p) != map_.end(); } 141 142 143 144 }; 145 146 147 #endif /* SHCalc_hpp */