1 // 2 // TMatrix.h 3 // pbsam_xcode 4 // 5 // Created by David Brookes on 6/15/16. 6 // Copyright © 2016 David Brookes. All rights reserved. 7 // 8 9 #ifndef TMatrix_h 10 #define TMatrix_h 11 12 #include <map> 13 #include <memory> 14 #include <vector> 15 #include <stdio.h> 16 #include "SHCalc.h" 17 #include "ReExpCalc.h" 18 #include "SystemSAM.h" 19 20 using namespace std; 21 22 /* 23 Re-expansion coefficients 24 */ 25 class TMatrix 26 { 27 protected: 28 29 /* 30 enum for telling the ReExpCoeffs which values to retrieve 31 in the below methods 32 */ 33 enum WhichReEx { BASE, FBASE, DDR, DDPHI, DDTHETA }; 34 35 int p_; 36 double kappa_; 37 vector<double> lam_scl_; // S factors, 0=kpio and 1=kpoo, for PB-SAM 38 39 vector<shared_ptr<ReExpCoeffs> > T_; 40 41 // maps (I,k), (J,l) indices to a single index in T. If this returns -1 42 // then the spheres are overlapping or less than 5A away and numerical 43 // re-expansion is required 44 map<vector<int>, int> idxMap_; 45 shared_ptr<SHCalc> _shCalc_; 46 shared_ptr<BesselCalc> _besselCalc_; 47 shared_ptr<SystemSAM> _system_; 48 49 int Nmol_; 50 vector<int> Nsi_; // number of spheres in each Molecule 51 52 53 // inner functions for re-expansion 54 MyMatrix<cmplx> expand_RX(MyMatrix<cmplx> X, 55 int I, int k, int J, int l, 56 WhichReEx whichR); 57 58 MyMatrix<cmplx> expand_SX(MyMatrix<cmplx> x1, 59 int I, int k, int J, int l, 60 WhichReEx whichS); 61 62 MyMatrix<cmplx> expand_RHX(MyMatrix<cmplx> x2, 63 int I, int k, int J, int l, 64 WhichReEx whichRH); 65 66 MyMatrix<cmplx> expand_dRdtheta_sing(MyMatrix<cmplx> mat, 67 int I, int k, int J, int l, 68 double theta, bool ham); 69 70 MyMatrix<cmplx> expand_dRdphi_sing(MyMatrix<cmplx> mat, 71 int I, int k, int J, int l, 72 double theta, bool ham); 73 74 // returns True if (J,l) > (I,k) (i.e. J > I or (I==J and k<l)) 75 bool is_Jl_greater(int I, int k, int J, int l); 76 77 /* 78 Convert derivatives from spherical to cartesian coords 79 */ 80 VecOfMats<cmplx>::type conv_to_cart(VecOfMats<cmplx>::type dZ, 81 int I, int k, int J, int l); 82 83 84 public: 85 TMatrix()86 TMatrix() { } 87 88 TMatrix(int p, shared_ptr<SystemSAM> _sys, shared_ptr<SHCalc> _shcalc, 89 shared_ptr<Constants> _consts, shared_ptr<BesselCalc> _besselcalc, 90 shared_ptr<ReExpCoeffsConstants> _reexpconsts); 91 92 // if these spheres can be re-expanded analytically, return true is_analytic(int I,int k,int J,int l)93 bool is_analytic(int I, int k, int J, int l) 94 { 95 if (idxMap_[{I, k, J, l}] == -1 ) return false; 96 else return true; 97 } 98 99 // get re-expansion of sphere (I, k) with respect to (J, l) get_T_Ik_Jl(int I,int k,int J,int l)100 shared_ptr<ReExpCoeffs> get_T_Ik_Jl(int I, int k, int J, int l) 101 { 102 return T_[idxMap_[{I, k, J, l}]]; 103 } 104 105 void update_vals(shared_ptr<SystemSAM> _sys, shared_ptr<SHCalc> _shcalc, 106 shared_ptr<BesselCalc> _besselcalc, 107 shared_ptr<ReExpCoeffsConstants> _reexpconsts); 108 109 /* 110 Re-expand a matrix X with respect to T(I,k)(J,l) 111 */ 112 MyMatrix<cmplx> re_expandX(MyMatrix<cmplx> X, int I, int k, int J, int l, 113 bool isF = false); 114 115 /* 116 Re-expand a numerical surface with respect to T(I,k)(J,l) (Equation 27b [1]) 117 */ 118 // MyMatrix<cmplx> re_expandX_numeric(vector<vector<double> > X, int I, int k, 119 // int J, int l, double kappa); 120 121 MyMatrix<cmplx> re_expandX_numeric(vector<vector<double> > X, int I, int k, 122 int J, int l, double kappa, 123 shared_ptr<PreCalcSH> pre_sh, 124 bool no_pre_sh=false); 125 126 /* 127 re-expand element j of grad(X) with element (I,k,J l) of T. REquires 128 the three components of grad(X) 129 */ 130 MyMatrix<Ptx> re_expand_gradX(MyMatrix<Ptx> dX, 131 int I, int k, int J, int l, bool isF = false); 132 133 134 135 /* 136 Re-expand X with element (I, k, J, l) of grad(T) and return 137 a matrix of Point objects containing each element of the gradient 138 */ 139 MyMatrix<Ptx> re_expandX_gradT(MyMatrix<cmplx> X, 140 int I, int k, 141 int J, int l); 142 143 /* 144 Locally re-expand X with element (I, k, J, l) of grad(T) and return 145 a matrix of Point objects containing each element of the gradient 146 */ 147 MyMatrix<Ptx> re_expandgradX_numeric(vector<vector<Pt> > X, 148 int I, int k, 149 int J, int l, double kappa, 150 shared_ptr<PreCalcSH> pre_sh, 151 bool no_pre_sh=false); 152 153 get_nmol()154 int get_nmol() const { return Nmol_; } get_nsi(int i)155 int get_nsi(int i) { return Nsi_[i]; } get_T_ct()156 int get_T_ct() { return (int) T_.size();} 157 compute_derivatives_i(int i)158 void compute_derivatives_i(int i) { T_[i]->calc_derivatives();} 159 160 161 // convert a matrix of Pts into a vector of 3 matrices 162 VecOfMats<cmplx>::type convert_from_ptx(MyMatrix<Ptx> X); 163 //and do the opposite of the above 164 MyMatrix<Ptx> convert_to_ptx(VecOfMats<cmplx>::type X); 165 166 }; 167 168 169 #endif /* TMatrix_h */ 170