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 }