1 //
2 //  TMatrix.cpp
3 //  pbsam_xcode
4 //
5 //  Created by David Brookes on 6/15/16.
6 //  Copyright © 2016 David Brookes. All rights reserved.
7 //
8 
9 #include "TMatrix.h"
10 
TMatrix(int p,shared_ptr<SystemSAM> _sys,shared_ptr<SHCalc> _shcalc,shared_ptr<Constants> _consts,shared_ptr<BesselCalc> _besselcalc,shared_ptr<ReExpCoeffsConstants> _reexpconsts)11 TMatrix::TMatrix(int p, shared_ptr<SystemSAM> _sys,
12                  shared_ptr<SHCalc> _shcalc,
13                  shared_ptr<Constants> _consts,
14                  shared_ptr<BesselCalc> _besselcalc,
15                  shared_ptr<ReExpCoeffsConstants> _reexpconsts)
16 :p_(p), kappa_(_consts->get_kappa()), Nmol_(_sys->get_n()), _system_(_sys),
17 _besselCalc_(_besselcalc), _shCalc_(_shcalc)
18 {
19   int total_spheres=0;
20   Nsi_ = vector<int> (Nmol_);
21   for (int I = 0; I < Nmol_; I++)
22   {
23     total_spheres += _sys->get_Ns_i(I);
24     Nsi_[I] = _sys->get_Ns_i(I);
25   }
26   T_.reserve(total_spheres*total_spheres);
27 
28   update_vals(_sys, _shcalc, _besselcalc, _reexpconsts);
29 }
30 
31 
update_vals(shared_ptr<SystemSAM> _sys,shared_ptr<SHCalc> _shcalc,shared_ptr<BesselCalc> _besselcalc,shared_ptr<ReExpCoeffsConstants> _reexpconsts)32 void TMatrix::update_vals(shared_ptr<SystemSAM> _sys, shared_ptr<SHCalc> _shcalc,
33                           shared_ptr<BesselCalc> _besselcalc,
34                           shared_ptr<ReExpCoeffsConstants> _reexpconsts)
35 {
36   T_.clear();
37   idxMap_.clear();
38 
39   int idx = 0;
40   double cutoff = 5.0;
41   Pt c_Ik, c_Jl, v;
42   double ak, al;
43   vector<int> idx_vec;
44   vector<double> kapVal(2);
45   for (int I = 0; I < _sys->get_n(); I++)
46   {
47     for (int J = 0; J < _sys->get_n(); J++)
48     {
49       for (int k = 0; k < _sys->get_Ns_i(I); k++)
50       {
51         c_Ik = _sys->get_centerik(I, k);
52         for (int l = 0; l < _sys->get_Ns_i(J); l++)
53         {
54           if (I==J && k==l) continue;
55 
56           idx_vec = {I, k, J, l};
57           c_Jl = _sys->get_centerik(J, l);
58           ak = _sys->get_aik(I, k);
59           al = _sys->get_aik(J, l);
60 
61           if ( (I==J) && ((c_Ik.dist(c_Jl)<cutoff)||
62                           (c_Ik.dist(c_Jl)<ak+al+cutoff)))
63           {
64             idxMap_[idx_vec] = -1;
65             continue;
66           }
67 
68 //          cout << "This is Ik, Jl pair : "  << I << ", "<< k
69 //              << " & "  << J << ", "<< l << " IK " << c_Ik.x() << ", "
70 //          << c_Ik.y() << ", "<< c_Ik.z() << " and " << c_Jl.x() << ", "
71 //          << c_Jl.y() << ", "<< c_Jl.z() << endl;
72 //          cout << " dist: " << c_Ik.dist(c_Jl) << " and aIk "
73 //              << ak << " and ajl "
74 //              << al << " dis1 : " <<  c_Ik.dist(c_Jl) - ak - al << endl;
75           if ( I == J ) kapVal = {0.0, kappa_};
76           kapVal = {kappa_, kappa_};
77           v = _sys->get_pbc_dist_vec_base(c_Ik, c_Jl);
78           vector<double> besselK = _besselcalc->calc_mbfK(2*p_,kapVal[1]*v.r());
79           _shcalc->calc_sh(v.theta(), v.phi());
80 
81           vector<double> lambdas = {_sys->get_aik(J, l), _sys->get_aik(I, k)};
82           auto re_exp = make_shared<ReExpCoeffs>(p_, v,
83                                                  _shcalc->get_full_result(),
84                                                  besselK, _reexpconsts,
85                                                  kapVal,
86                                                  lambdas, false);
87           T_.push_back(re_exp);
88           idxMap_[idx_vec] = idx;
89           idx++;
90         }
91       }
92     }
93   }
94 }
95 
96 //// Perform local expansion from J, l onto I, k
97 //MyMatrix<cmplx> TMatrix::re_expandX_numeric(vector<vector<double> > X,
98 //                                          int I, int k,
99 //                                          int J, int l, double kappa)
100 //{
101 //  int h, n, m;
102 //  cmplx val, sh;
103 //  double chgscl, rscl, ekr;
104 //  MyMatrix<cmplx> Z(p_, 2*p_+1);
105 //  vector<int> exp_pts = _system_->get_gdpt_expij(J, l);
106 //
107 //
108 //  for (h = 0; h < X[l].size(); h++)
109 //  {
110 //    Pt sph_dist = _system_->get_centerik(I, k) - _system_->get_centerik(J, l);
111 //    Pt loc = _system_->get_gridijh(J, l, exp_pts[h]) - sph_dist;
112 //
113 //
114 //    _shCalc_->calc_sh(loc.theta(), loc.phi());
115 //
116 //    vector<double> bessI = _besselCalc_->calc_mbfK(p_+1, kappa*loc.r());
117 //    rscl = _system_->get_aik(I, k) / loc.r();
118 //    chgscl = X[l][h] / loc.r();
119 //    ekr = exp(-kappa*loc.r());
120 //    for (n = 0; n < p_; n++)
121 //    {
122 //      for (m = -n; m <= n; m++)
123 //      {
124 //        sh = _shCalc_->get_result(n, m);
125 //        val = bessI[n] * ekr * chgscl * sh + Z(n, m+p_);
126 //        Z.set_val(n, m+p_, val);
127 //      }
128 //      chgscl *= rscl;
129 //    }
130 //  }
131 //  return Z;
132 //}
133 
134 
re_expandX_numeric(vector<vector<double>> X,int I,int k,int J,int l,double kappa,shared_ptr<PreCalcSH> pre_sh,bool no_pre_sh)135 MyMatrix<cmplx> TMatrix::re_expandX_numeric(vector<vector<double> > X,
136                                             int I, int k,
137                                             int J, int l, double kappa,
138                                             shared_ptr<PreCalcSH> pre_sh,
139                                             bool no_pre_sh)
140 {
141   int h, n, m;
142   cmplx val, sh;
143   double chgscl, rscl, ekr;
144   MyMatrix<cmplx> Z(p_, 2*p_+1);
145   vector<int> exp_pts = _system_->get_gdpt_expij(J, l);
146   for (h = 0; h < X[l].size(); h++)
147   {
148     Pt sph_dist = _system_->get_centerik(I, k) - _system_->get_centerik(J, l);
149     Pt loc = _system_->get_gridijh(J, l, exp_pts[h]) - sph_dist;
150 
151     if (no_pre_sh) pre_sh->calc_and_add(loc, _shCalc_);
152 
153     vector<double> bessI = _besselCalc_->calc_mbfK(p_+1, kappa*loc.r());
154     rscl = _system_->get_aik(I, k) / loc.r();
155     chgscl = X[l][h] / loc.r();
156     ekr = exp(-kappa*loc.r());
157 
158     for (n = 0; n < p_; n++)
159     {
160       for (m = -n; m <= n; m++)
161       {
162         sh = pre_sh->get_sh(loc, n, m);
163         val = bessI[n] * ekr * chgscl * sh + Z(n, m+p_);
164         Z.set_val(n, m+p_, val);
165       }
166       chgscl *= rscl;
167     }
168   }
169   return Z;
170 }
171 
172 
173 
re_expandX(MyMatrix<cmplx> X,int I,int k,int J,int l,bool isF)174 MyMatrix<cmplx> TMatrix::re_expandX(MyMatrix<cmplx> X,
175                                     int I, int k,
176                                    int J, int l, bool isF)
177 {
178   MyMatrix<cmplx> X1, X2, Z;
179   WhichReEx whichR=BASE, whichS=BASE, whichRH=BASE;
180 
181   if (isF) whichS = FBASE;
182   X1 = expand_RX(X, I, k, J, l, whichR);
183   X2 = expand_SX(X1, I, k, J, l, whichS);
184   Z  = expand_RHX(X2, I, k, J, l, whichRH);
185   return Z;
186 }
187 
188 
189 /*
190  re-expand element j of grad(X) with element (i, j) of T. Requires
191  the three components of grad(X)
192  */
re_expand_gradX(MyMatrix<Ptx> dX,int I,int k,int J,int l,bool isF)193 MyMatrix<Ptx> TMatrix::re_expand_gradX(MyMatrix<Ptx> dX,
194                                         int I, int k,
195                                         int J, int l, bool isF)
196 
197 {
198   VecOfMats<cmplx>::type dX_comps = convert_from_ptx(dX);
199 
200   MyMatrix<cmplx> x1, x2, z;
201   VecOfMats<cmplx>::type Z (3);
202   WhichReEx whichR=BASE, whichS=BASE, whichRH=BASE;
203 
204   if (isF) whichS = FBASE;
205   // dA/dR = 0, dA/dtheta = 1, dA/dphi = 2
206   for (int d = 0; d < 3; d++)
207   {
208     x1 = expand_RX(dX_comps[d], I, k, J, l, whichR);
209     x2 = expand_SX(x1, I, k, J, l, whichS);
210     z  = expand_RHX(x2, I, k, J, l, whichRH);
211     Z.set_val(d, z);
212   }
213 
214 //  MyMatrix<Ptx> Zpt = convert_to_ptx(Z);
215 //
216 //  cout << "after re_exp_grad" << endl;
217 //  for (int d = 0; d < 3; d++)
218 //  {
219 //    cout << " Dim: " << d <<  endl;
220 //    for (int n = 0; n < p_; n++)
221 //    {
222 //      for (int m = 0; m <= n; m++)
223 //      {
224 //        cout << Zpt(n,m+p_).get_cart(d) << ", ";
225 //      }
226 //      cout << endl;
227 //    }
228 //  }
229 //  cout << endl;
230 //
231 //  return Zpt;
232 
233   return convert_to_ptx(Z);
234 }
235 
236 
re_expandX_gradT(MyMatrix<cmplx> X,int I,int k,int J,int l)237 MyMatrix<Ptx> TMatrix::re_expandX_gradT(MyMatrix<cmplx> X,
238                                          int I, int k,
239                                          int J, int l)
240 {
241   MyMatrix<cmplx> x1, x2, z, z1, z2;
242   VecOfMats<cmplx>::type Z (3); // output vector
243   WhichReEx whichR=BASE, whichS=BASE, whichRH=BASE;
244 
245 
246   // first find with respect to dT/dr:
247   whichS = DDR;
248   x1 = expand_RX(X, I, k, J, l, whichR);
249   x2 = expand_SX(x1, I, k, J, l, whichS);
250   z  = expand_RHX(x2, I, k, J, l, whichRH);
251   Z.set_val(0, z);
252 
253   // dT/dtheta:
254   whichR=BASE;
255   whichS = BASE;
256   whichRH = DDTHETA;
257   x1 = expand_RX(X, I, k, J, l, whichR);
258   x2 = expand_SX(x1, I, k, J, l, whichS);
259   z1 = expand_RHX(x2, I, k, J, l, whichRH);
260 
261   whichRH = BASE;
262   whichR = DDTHETA;
263   x1 = expand_RX(X, I, k, J, l, whichR);
264   x2 = expand_SX(x1, I, k, J, l, whichS);
265   z2 = expand_RHX(x2, I, k, J, l, whichRH);
266   Z.set_val(1, z1 + z2);
267 
268   // dT/dphi:
269   whichR = BASE;
270   whichRH = DDPHI;
271   x1 = expand_RX(X, I, k, J, l, whichR);
272   x2 = expand_SX(x1, I, k, J, l, whichS);
273   z1 = expand_RHX(x2, I, k, J, l, whichRH);
274 
275   whichRH = BASE;
276   whichR = DDPHI;
277   x1 = expand_RX(X, I, k, J, l, whichR);
278   x2 = expand_SX(x1, I, k, J, l, whichS);
279   z2 = expand_RHX(x2, I, k, J, l, whichRH);
280   Z.set_val(2, z1 + z2);
281 
282   Z = conv_to_cart(Z, I, k, J, l);
283   return convert_to_ptx(Z);
284 }
285 
286 // Perform local expansion from J, l onto I, k
re_expandgradX_numeric(vector<vector<Pt>> X,int I,int k,int J,int l,double kappa,shared_ptr<PreCalcSH> pre_sh,bool no_pre_sh)287 MyMatrix<Ptx> TMatrix::re_expandgradX_numeric(vector<vector<Pt> > X,
288                                               int I, int k,
289                                               int J, int l, double kappa,
290                                               shared_ptr<PreCalcSH> pre_sh,
291                                               bool no_pre_sh)
292 {
293   int h, n, m;
294   cmplx val, sh;
295   double chgscl, rscl, ekr, xval;
296   VecOfMats<cmplx>::type Z (3, MyMatrix<cmplx> (p_, 2*p_+1));
297   vector<int> exp_pts = _system_->get_gdpt_expij(J, l);
298 
299   for (h = 0; h < X[l].size(); h++)
300   {
301     if (X[l][h].norm2() < 1e-15) continue;
302     for (int d = 0; d < 3; d++)
303     {
304       xval = X[l][h].get_cart(d);
305 
306       Pt sph_dist = _system_->get_centerik(I, k) - _system_->get_centerik(J, l);
307       Pt loc = _system_->get_gridijh(J, l, exp_pts[h]) - sph_dist;
308 //      _shCalc_->calc_sh(loc.theta(), loc.phi());
309       if (no_pre_sh) pre_sh->calc_and_add(loc, _shCalc_);
310 
311       vector<double> bessI = _besselCalc_->calc_mbfK(p_+1, kappa*loc.r());
312       rscl = _system_->get_aik(I, k) / loc.r();
313       chgscl = xval / loc.r();
314       ekr = exp(-kappa*loc.r());
315       for (n = 0; n < p_; n++)
316       {
317         for (m = -n; m <= n; m++)
318         {
319           sh = pre_sh->get_sh(loc, n, m);
320           val = bessI[n]*ekr*chgscl*sh + Z[d](n, m+p_);
321           Z[d](n, m+p_) = val;
322         }
323         chgscl *= rscl;
324       }
325     }
326   }
327   return convert_to_ptx(Z);
328 }
329 
is_Jl_greater(int I,int k,int J,int l)330 bool TMatrix::is_Jl_greater(int I, int k, int J, int l)
331 {
332   if (J > I || (I==J && l > k)) return true;
333   else return false;
334 }
335 
336 // perform first part of T*A and return results
expand_RX(MyMatrix<cmplx> X,int I,int k,int J,int l,WhichReEx whichR)337 MyMatrix<cmplx> TMatrix::expand_RX(MyMatrix<cmplx> X,
338                                    int I, int k, int J, int l,
339                                    WhichReEx whichR)
340 {
341 
342   int n, m, s, map_idx;
343   map_idx = idxMap_[{I, k, J, l}];
344 
345   MyMatrix<cmplx> x1(p_, 2*p_ + 1);
346   cmplx inter, rval, aval;
347 
348   // fill X1:
349   for (n = 0; n < p_; n++)
350   {
351     for (m = -n; m <= n; m++)
352     {
353       inter  = 0;
354       if (T_[map_idx]->isSingular())
355       {
356         Pt vec = T_[map_idx]->get_TVec();
357         if (whichR == BASE)
358         {
359           aval = X(n, -m+p_);
360           if (vec.theta() > M_PI/2.0)
361             x1.set_val(n, m+p_, (n%2 == 0 ? aval : -aval));
362           else x1.set_val(n, m+p_, X(n, m+p_));
363         } else if (whichR == DDTHETA)
364         {
365           x1 = expand_dRdtheta_sing(X, I, k, J, l, vec.theta(), false);
366           return x1;
367         }
368         else
369         {
370           x1 = expand_dRdphi_sing(X, I, k, J, l, vec.theta(), false);
371           return x1;
372         }
373       } else
374       {
375         for (s = -n; s <= n; s++)
376         {
377           if (whichR == DDPHI)
378             rval = T_[map_idx]->get_dr_dphi_val(n, m, s);
379           else if (whichR == DDTHETA)
380             rval = T_[map_idx]->get_dr_dtheta_val(n, m, s);
381           else
382             rval = T_[map_idx]->get_rval(n, m, s);
383 
384           aval = X(n, s+p_);
385 
386           inter += rval * aval;
387         } // end s
388         x1.set_val(n, m+p_, inter);
389       }
390     } // end m
391   } //end n
392   return x1;
393 }
394 
395 // perform second part of T*A and return results
expand_SX(MyMatrix<cmplx> x1,int I,int k,int J,int l,WhichReEx whichS)396 MyMatrix<cmplx> TMatrix::expand_SX(MyMatrix<cmplx> x1,
397                                    int I, int k, int J, int l,
398                                    WhichReEx whichS)
399 {
400   double fac;
401   cmplx inter, sval;
402   MyMatrix<cmplx> x2(p_, 2*p_ + 1);
403 
404   int n, m, s, map_idx = idxMap_[{I, k, J, l}];
405   vector<double> lam = T_[map_idx]->get_lambdas();
406   vector<double> lamScl = T_[map_idx]->get_lam_scale();
407 
408   // fill x2:
409   for (n = 0; n < p_; n++)
410   {
411     for (m = -n; m <= n; m++)
412     {
413       inter  = 0;
414       for (s = abs(m); s < p_; s++)
415       {
416         fac = ( s <= n ) ? lamScl[n-s] : 1.0;
417         if (whichS == DDR) sval = T_[map_idx]->get_dsdr_val(n, s, m);
418         else if (whichS == FBASE) sval = T_[map_idx]->get_s_fval(n, s, m);
419         else sval = T_[map_idx]->get_sval(n, s, m);
420         inter += fac * sval * x1(s, m+p_);
421       } // end l
422       x2.set_val(n, m+p_, inter);
423     } // end m
424   } //end n
425   return x2;
426 }
427 //
428 // perform third part of T*A and return results
expand_RHX(MyMatrix<cmplx> x2,int I,int k,int J,int l,WhichReEx whichRH)429 MyMatrix<cmplx> TMatrix::expand_RHX( MyMatrix<cmplx> x2,
430                                     int I, int k, int J, int l,
431                                     WhichReEx whichRH)
432 {
433   int n, m, s, map_idx = idxMap_[{I, k, J, l}];
434   cmplx inter, rval;
435   MyMatrix<cmplx> z(p_, 2*p_ + 1);
436 
437 //  bool jl_greater = is_Jl_greater(I, k, J, l);
438 //  if (jl_greater) map_idx = idxMap_[{I, k, J, l}];
439 //  else            map_idx = idxMap_[{J, l, I, k}];
440 
441   //fill zj:
442   for (n = 0; n < p_; n++)
443   {
444     for (m = -n; m <= n; m++)
445     {
446       inter  = 0;
447       if (T_[map_idx]->isSingular())
448       {
449         Pt vec = T_[map_idx]->get_TVec();
450         if (whichRH == BASE)
451         {
452           if (vec.theta() > M_PI/2.0)
453             z.set_val(n, m+p_, (n%2 == 0 ? x2(n,-m+p_) : -x2(n,-m+p_)));
454           else
455             z.set_val(n, m+p_, x2(n,m+p_));
456         } else if (whichRH == DDTHETA)
457         {
458           z = expand_dRdtheta_sing(x2, I, k, J, l, vec.theta(), true);
459           return z;
460         }
461         else
462         {
463           z = expand_dRdphi_sing(x2, I, k, J, l, vec.theta(), true);
464           return z;
465         }
466       } else
467       {
468         for (s = -n; s <= n; s++)
469         {
470           if (whichRH == DDPHI)
471             rval = T_[map_idx]->get_dr_dphi_val(n, s, m);
472           else if (whichRH == DDTHETA)
473             rval = T_[map_idx]->get_dr_dtheta_val(n, s, m);
474           else
475             rval = T_[map_idx]->get_rval(n, s, m);
476 
477           inter += conj(rval) * x2(n, s+p_);
478         } // end s
479         z.set_val(n, m+p_, inter);
480       }
481     } // end m
482   } //end n
483   return z;
484 }
485 
expand_dRdtheta_sing(MyMatrix<cmplx> mat,int I,int k,int J,int l,double theta,bool ham)486 MyMatrix<cmplx> TMatrix::expand_dRdtheta_sing(MyMatrix<cmplx> mat,
487                                               int I, int k, int J, int l,
488                                               double theta, bool ham)
489 {
490   MyMatrix<cmplx> x(p_, 2*p_ + 1);
491   x.set_val( 0, p_, cmplx(0.0, 0.0));
492   double rec = (ham ? -1.0 : 1.0);
493   int map_idx;
494 
495   bool jl_greater = is_Jl_greater(I, k, J, l);
496   if (jl_greater) map_idx = idxMap_[{I, k, J, l}];
497   else            map_idx = idxMap_[{J, l, I, k}];
498 
499   if (theta < M_PI/2)
500   {
501     for (int n = 1; n < p_; n++)
502     {
503       x.set_val( n, p_, rec*cmplx(2.0*T_[map_idx]->get_prefac_dR_val(n,0,1)*
504                                   mat(n, 1+p_).real(), 0.0)); // m = 0
505       for (int m = 1; m < n; m++)
506       {
507         x.set_val( n, m+p_, rec*T_[map_idx]->get_prefac_dR_val(n,m,0)*mat(n,m-1+p_)
508                   + rec*T_[map_idx]->get_prefac_dR_val(n,m,1)*mat(n,m+1+p_));
509         x.set_val( n,-m+p_, conj( x( n, m+p_)));
510       }
511 
512       x.set_val(n, n+p_,
513                 rec*T_[map_idx]->get_prefac_dR_val( n, n, 0)*mat(n, n-1+p_));
514       x.set_val(n,-n+p_, conj( x( n, n+p_)));
515     }
516   }
517   else
518   {
519     double s = -1.0;
520     for (int n = 1; n < p_; n++, s = -s)
521     {
522       x.set_val( n, p_, rec*cmplx(2.0*s*T_[map_idx]->get_prefac_dR_val(n,0,1)*
523                                   mat(n,1+p_).real(), 0.0)); // m = 0
524       for (int m = 1; m < n; m++)
525       {
526         x.set_val( n, m+p_, rec*s*
527                   (T_[map_idx]->get_prefac_dR_val(n,m,0)*mat(n,-m+1+p_)
528                    + T_[map_idx]->get_prefac_dR_val(n,m,1)*mat(n,-m-1+p_)));
529         x.set_val( n,-m+p_, conj( x( n, m+p_)));
530       }
531 
532       x.set_val(n, n+p_,rec*s*T_[map_idx]->get_prefac_dR_val(n, n,0)*mat(n,-n+1+p_));
533       x.set_val(n,-n+p_, conj( x( n, n+p_)));
534     }
535   }
536   return x;
537 }
538 
expand_dRdphi_sing(MyMatrix<cmplx> mat,int I,int k,int J,int l,double theta,bool ham)539 MyMatrix<cmplx> TMatrix::expand_dRdphi_sing(MyMatrix<cmplx> mat,
540                                             int I, int k, int J, int l,
541                                             double theta, bool ham)
542 {
543   MyMatrix<cmplx> x(p_, 2*p_ + 1);
544   x.set_val( 0, p_, cmplx(0.0, 0.0));
545   double rec = ((ham && (theta < M_PI/2)) ? -1.0 : 1.0);
546   int map_idx;
547 
548   bool jl_greater = is_Jl_greater(I, k, J, l);
549   if (jl_greater) map_idx = idxMap_[{I, k, J, l}];
550   else            map_idx = idxMap_[{J, l, I, k}];
551 
552 
553   if (theta < M_PI/2)
554   {
555     for (int n = 1; n < p_; n++)
556     {
557       x.set_val( n, p_, rec*cmplx(2.0*T_[map_idx]->get_prefac_dR_val(n,0,1)*
558                                   mat(n, 1+p_).imag(),0.0));
559       for (int m = 1; m < n; m++)
560       {
561         x.set_val(n, m+p_,
562                   rec*(cmplx( 0.0, T_[map_idx]->get_prefac_dR_val(n,m, 0))*mat(n,m-1+p_)
563                        - cmplx( 0.0, T_[map_idx]->get_prefac_dR_val(n,m,1))*mat(n,m+1+p_)));
564         x.set_val( n, -m+p_, conj(x(n, m+p_)));
565       }
566 
567       x.set_val( n, n+p_,
568                 rec*cmplx( 0.0, T_[map_idx]->get_prefac_dR_val(n,n,0))*mat(n,n-1+p_));
569       x.set_val( n, -n+p_, conj(x( n, n+p_)));
570     }
571   }
572   else
573   {
574     double s = 1.0;
575     for (int n = 1; n < p_; n++, s = -s)
576     {
577       x.set_val(n, p_, rec*cmplx(2.0*s*T_[map_idx]->get_prefac_dR_val(n,0,1)*
578                                  mat(n,1+p_).imag(),0.0));
579       for (int m = 1; m < n; m++)
580       {
581         x.set_val(n, m+p_,
582                   rec*s*(-cmplx(0.0,T_[map_idx]->get_prefac_dR_val(n,m,0))*mat(n,-m+1+p_)
583                          + cmplx(0.0, T_[map_idx]->get_prefac_dR_val(n,m,1))*mat(n,-m-1+p_)));
584         x.set_val( n, -m+p_, conj(x(n, m+p_)));
585       }
586 
587       x.set_val( n,  n+p_, rec*cmplx(0.0,
588                                      -s*T_[map_idx]->get_prefac_dR_val(n,n,0))*mat(n,-n+1+p_));
589       x.set_val( n, -n+p_, conj(x(n, n+p_)));
590     }
591   }
592   return x;
593 }
594 
595 
conv_to_cart(VecOfMats<cmplx>::type dZ,int I,int k,int J,int l)596 VecOfMats<cmplx>::type TMatrix::conv_to_cart(VecOfMats<cmplx>::type dZ,
597                                              int I, int k, int J, int l)
598 {
599   int m, n, map_idx;
600   double the, r, phi;
601   VecOfMats<cmplx>::type Zcart (3); vector<double> con1(3), con2(3), con3(3);
602 
603   map_idx = idxMap_[{I, k, J, l}];
604   Pt v = T_[map_idx]->get_TVec();
605   the = v.theta(); r = v.r(); phi = v.phi();
606 
607   if (T_[map_idx]->isSingular())
608   {
609     double cost = (the < M_PI/2) ? 1.0 : -1.0;
610     con1 = {  0.0, cost/r,   0.0};
611     con2 = {  0.0,    0.0, 1.0/r};
612     con3 = { cost,    0.0,   0.0};
613   } else
614   {
615     con1 = {sin(the)*cos(phi),  cos(the)*cos(phi)/r, -sin(phi)/sin(the)/r};
616     con2 = {sin(the)*sin(phi),  cos(the)*sin(phi)/r,  cos(phi)/sin(the)/r};
617     con3 = {         cos(the),          -sin(the)/r,                  0.0};
618   }
619 
620   for (n = 0; n < 3; n++)
621     Zcart[n] = MyMatrix<cmplx> (p_, 2*p_ + 1);
622 
623   for ( n = 0; n < p_; n++ )
624     for ( m = -n; m <= n; m++ )
625     {
626       (&Zcart[0])->set_val( n, m+p_, con1[0]*dZ[0](n, m+p_)
627                            + con1[1]*dZ[1](n, m+p_)+con1[2]*dZ[2](n, m+p_));
628       (&Zcart[1])->set_val( n, m+p_, con2[0]*dZ[0](n, m+p_)
629                            + con2[1]*dZ[1](n, m+p_)+con2[2]*dZ[2](n, m+p_));
630       (&Zcart[2])->set_val( n, m+p_, con3[0]*dZ[0](n, m+p_)
631                            + con3[1]*dZ[1](n, m+p_)+con3[2]*dZ[2](n, m+p_));
632     }
633   return Zcart;
634 }
635 
636 
637 // convert a matrix of Point objects into a vector of 3 matrices
convert_from_ptx(MyMatrix<Ptx> X)638 VecOfMats<cmplx>::type TMatrix::convert_from_ptx(MyMatrix<Ptx> X)
639 {
640   VecOfMats<cmplx>::type result (3, MyMatrix<cmplx> (X.get_nrows(),
641                                                      X.get_ncols()));
642   for (int i = 0 ; i < 3; i++)
643     for (int j = 0; j < X.get_nrows(); j++)
644       for (int k = 0; k < X.get_ncols(); k++)
645         result[i].set_val(j, k, X(j, k)[i]);
646 
647   return result;
648 
649 }
650 
651 
convert_to_ptx(VecOfMats<cmplx>::type X)652 MyMatrix<Ptx> TMatrix::convert_to_ptx(VecOfMats<cmplx>::type X)
653 {
654   MyMatrix<Ptx> result (X[0].get_nrows(), X[0].get_ncols());
655   for (int j = 0; j < X[0].get_nrows(); j++)
656     for (int k = 0; k < X[0].get_ncols(); k++)
657     {
658       result(j, k).set_x(X[0](j, k));
659       result(j, k).set_y(X[1](j, k));
660       result(j, k).set_z(X[2](j, k));
661     }
662   return result;
663 }