laguerre2bessel.cc 18 KB
Newer Older
1 2 3
#include "laguerre2bessel.h"
#include "lagsht_bessel.h"
#include "laguerreTransform.h"
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
4 5 6
#include "lagsht_exceptions.h"


7

8
#include "walltimer.h"    //timing
9
#include <fstream>       //dump
10
#include <sstream>
11
#include <numeric>
12 13


14
#define DEBUG 0
15
#define LIMBER 0
16

17 18
namespace LagSHT {

19 20
  /* Version 0: test de calcul des J_{lnp} via une quadrature de Laguerre 
     mais ca n'a pas l'air de converger (28 Sept 15)
21
     Version 1: on fait une integration a la Clenshaw-Curtis avec startegie Localeo ou Globale
22
     recursive
23
     Version 2: Jlnp calcule via une Discrete Fourier-Bessel Transform and a  Clenshaw-Curtis quadrature and global integration strategy. 
24 25
     Version 3: calcul via une quadrature de  Clenshaw-Curtis quadrature and global integration en controlant les bornes
                d'integration
26

27
  */
28 29 30 31 32 33 34 35
Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
				 LaguerreTransform* lagTrans,
				 int Nmax, int Pmax, r_8 Rmax,
				 string clenshawDir,
				 string jlnpDir,
				 bool recomputeJlnp
				 ):
  sphere_(sphere),lagTrans_(lagTrans),
36 37
  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
  //JEC 1/2/16  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmaax), PPmax_(Pmax) {
38 39 40 41 42 43 44 45 46 47 48
  
  if(sphere_) {
    //store localy usefull parameters for the algorithms Analysis/Synthesis
    Npix_ = sphere_->Npix();
    Nalm_ = sphere_->Nalm();
    Lmax_ = sphere_->Lmax();
    Mmax_ = sphere_->Mmax();
  } else {
    throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
  }
  
49
  //  jnu_ = new Bessel(Lmax_,std::max(std::max(Pmax_,PPmax_),Nmax_));  //the max is here just in case PPmax =/= Pmax
50 51 52
  
  cout << "(JEC) lag2bes Ctor  Compute jnu_ START" << endl;

53
  jnu_ = new Bessel(Lmax_,std::max(Pmax_,Nmax_));  //the max is here just in case Nmax =/= Pmax
54
  
55 56 57
  cout << "(JEC) lag2bes Ctor  Compute jnu_ END" << endl;

  
58 59
  if(!recomputeJlnp){
    
60 61
    cout << "READ mode......." << endl; 

62
    //read the Rmax-independant values of the jlnp and rescale
63 64
    //    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
    stringstream ss; ss << "limber0-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
65
    
66 67
    string fname(jlnpDir+"/"+ss.str());
    std::ifstream ifs(fname.c_str(), std::ifstream::in);
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
    if(!ifs.is_open())
      throw LagSHTError("ERROR: jlnp file open failed.:"+ss.str());
    
    r_8 tau = lagTrans_->GetTau();
    r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
    int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
    Jlnp_.resize(Nmax_*Jstride);
    
    int pcur, lcur, ncur;
    r_8 val;
    for(int p=0; p<Pmax_; p++){
      for(int l=0; l<Lmax_; l++){
	int ploff7 = p + l*Pmax_;
	for (int n=0; n<Nmax_; n++){
	  ifs >> lcur >> ncur >> pcur >> val;         //read
83 84
 	  if( (l!=lcur) || (n!=ncur) || (p!=pcur))                         //check
 	    throw LagSHTError("ERROR: jlnp file read failed.:"+ss.str());
85 86 87 88 89 90
	  Jlnp_[ploff7 + n*Jstride] =  tau3sqrt * val;                     //rescaling
	}//n
      }//l
    }//p
    
  } else {
91 92 93

    cout << "Compute mode......." << endl; 

94 95 96
    ComputeJInteg(clenshawDir,jlnpDir);
  }//jlnp omputation/load
  
97 98
}//Ctor

99

100

101

102 103 104 105 106 107
/*!
  Computation of the Jlnp = Jln(k_{lp}) = Sqrt(2/pi) tau^3/2 Int_0^infty jl(k_{lp} tau x) K_n(x) x^2 dx
  with
  K_n(x) = Exp[-x/2] Ln(x) Sqrt(n!/(n+alpha)!) ->cf. LaguerreFuncQuad
  k_{lp}*tau = q_{lp}/r_{Nmax-1}
  
108 109 110 111 112 113 114 115 116 117 118 119 120
  Use a mixed version:
  1- compute Limber approx (order 0)
  2- compute higher order expensiuon (order 3)
  If |order0-order3|<tol Then 
     result = order3
  Else
  Clenshaw-Curtis quadrature with 2*40-1 pts between LowBound and UppBound
      - LowBound is (50% * l_value) or( 0 if l_value<10)
      - UppBound is the Max between : 
      - for n>=2 nodes[n-1] + 5*(nodes[n-1]-nodes[n-2]) this is an estimate of Kn cut-off overwise use the highest Laguerre Node
      - q(l,1)/(k_{lp}*tau) this is an estimate of first oscillation of Bessel function  
  Endif

121
*/
122

123 124
void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){

125 126
  cout << "(JEC) ComputeJInteg START " << endl;

127
#if LIMBER>0
128 129
  r_8 h2=1.0; //put to 0 to switch off the 2nd order of Limber expension approx
  r_8 h3=1.0; //put to 0 to switch off the 3rd order of Limber expension approxxs
130
#endif
131 132 133 134 135 136 137 138 139

  if(Nmax_<2)throw LagSHTError("Laguerre2Bessel::ComputeJInteg Nmax<2");

  tstack_push("J Integ compute ");


  r_8 tau = lagTrans_->GetTau();
  r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
  r_8 rNm1 = Rmax_ / tau; // R = r[N-1] with N the original transform value. Todo (lagTrans_->GetNodes()).back()
140
  r_8 tol4Integral = 1e-10;
141 142

#if LIMBER>0
143
  r_8 tol4Limber   = 1e-6;
144
#endif
145

146 147 148 149 150
  typedef  ClenshawCurtisQuadrature<double> Integrator_t;
  string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-40.txt";
  Integrator_t theQuad(40,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
  
  Quadrature<r_8,Integrator_t>::values_t integ_val;
151 152


153
  int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
154

155 156 157
  Jlnp_.resize(Nmax_*Jstride);  

  r_8 alpha = lagTrans_->GetAlpha();
158 159 160 161 162 163 164
  r_8 alphaOv2M1 = alpha/2.0-1.0;
  r_8 alpha2  = alpha*alpha;
  r_8 alphap1 = alpha+1.0;


  cout << "(JEC) ComputeJInteg Facts START " << endl;

165 166 167 168
  vector<r_8> facts(Nmax_); //sqrt(n!/(n+alpha)!) the normalization
  facts[0] = 1.0/sqrt(boost::math::tgamma(alpha+1.0)); //1/sqrt(alpha)!       
  for(int n=1;n<Nmax_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) );

169 170
  cout << "(JEC) ComputeJInteg Facts END " << endl;

171
 
172 173
  cout << "(JEC) ComputeJInteg  integUpperBound START " << endl;

174 175 176 177 178 179 180 181 182 183 184
  vector<r_8> integUpperBound(Nmax_);
  integUpperBound[0] = rNm1; //this is the case of L_0^{\alpha}(x) =1 with no root
  integUpperBound[1] = rNm1; //this is the case of L_1^{\alpha}(x) =1+alpha-x with single root = 1/(1+alpha)
  for (int n=2;n<Nmax_;n++){
    LaguerreFuncQuad lag(n,alpha);
    vector<r_8> nodes;
    vector<r_8> weights;
    lag.QuadWeightNodes(nodes,weights);
    integUpperBound[n] = nodes[n-1] + 5*(nodes[n-1]-nodes[n-2]); //upper bound from Laguerre function end point, may be modified by Bessel part later
  }

185 186
  cout << "(JEC) ComputeJInteg  integUpperBound END " << endl;

187

188
  stringstream ss; ss << "NR-NoLimber-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
189 190 191
  std::ofstream ofs;
  string fname(jlnpDir+"/"+ss.str());
  ofs.open (fname.c_str(), std::ofstream::out);
192

193
  cout << "(JEC) ComputeJInteg write on <" <<fname << ">"<< endl;
194

195 196
  r_8 sqrt2Ovpi= sqrt(2.0/M_PI); 
  
197
  for(int p=0; p<Pmax_; p++){
198
    //    cout << "p = " << p;
199

200
    for(int l=0; l<Lmax_; l++){
201
      //      cout << " l = " << l;
202 203

      int ploff7 = p + l*Pmax_;
204

205
      r_8 qlp = (*jnu_)(l,p);                //ql,p
206 207 208 209 210 211 212 213 214 215 216 217 218
      r_8 klptau = qlp/rNm1;                 //klp * tau
      r_8 klptau2 = klptau*klptau;
      r_8 klptau3 = klptau*klptau2;
      r_8 klptau6 = klptau3*klptau3;

      r_8 r8l = (r_8)(l);
      
      r_8 fl = 1.0+2.0*r8l;
      r_8 fl2 = fl*fl;
      r_8 fl3 = fl2*fl;
      
      

219
      r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau; //JEC 4/2/16 for low l just start at 0   
220

221 222 223
    
      r_8 ql2scaled = (*jnu_)(l,1)/klptau;   //  the 2nd BesselRoot rescaled

224 225

#if LIMBER>0      
226
      r_8 arg = (r8l+0.5)/klptau;            // (L+1/2)/(klp*tau) for Limber development
227

228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
      r_8 limb0Norm =  (1./klptau3) * pow((r8l+0.5),(r_8)1.5) * pow(arg,alphaOv2M1);
      
      r_8 limb3C0  = -(h3*(8*(-3 + alpha)*(-1 + alpha)*alphap1*klptau3 
		     + 12*(-1 + alpha2)*klptau2*fl + 
		     6*alphap1*klptau*fl2 + fl3)) 
	       - 12*klptau*(-8*klptau2*fl2 + h2*(4*(-1 + alpha2)*klptau2 + 4*alphap1*klptau*fl + fl2) );

      r_8 limb3C1 = 6*fl * ( 8*h2*klptau*(1 + 2*alphap1*klptau + 2*r8l) + 
			     h3*(4*(-1 + alpha2)*klptau2 + 4*alphap1*klptau*fl + fl2)
			     );

      r_8 limb3C2 = - 12*fl2 * ( 4*h2*klptau + h3*(1 + 2*alphap1*klptau + 2*r8l));

      r_8 limb3C3 = 8*h3*fl3;

      r_8 limb3Norm = pow(arg,alphaOv2M1)/(384 * klptau6 * pow(0.5+r8l,0.5));
244
#endif
245

246
      for (int n=0; n<Nmax_; n++){
247
	//	cout << " n = " << n << endl;
248 249 250

	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);

251 252 253
	r_8 JlnpApprox = 0; 
#if LIMBER>0
	//1. compute the Limber (Order 0) approximation
254
	r_8 limber0 = limb0Norm * Ln->Value(arg)*facts[n];
255
	
256 257
	//2.  compute a higher order Limber (Order 3) approximation
	//h2: term of 2nd order; h3: term of 3rd order
258

259 260 261 262 263 264 265 266
	LaguerreFuncQuad* Ln_ap1 = new LaguerreFuncQuad(n, alpha+1);
	LaguerreFuncQuad* Ln_ap2 = new LaguerreFuncQuad(n, alpha+2);
	LaguerreFuncQuad* Ln_ap3 = new LaguerreFuncQuad(n, alpha+3);
	
	r_8 lag0 = Ln->Value(arg);
	r_8 lag1 = Ln_ap1->Value(arg);
	r_8 lag2 = Ln_ap2->Value(arg);
	r_8 lag3 = Ln_ap3->Value(arg);
267 268

	
269 270 271 272
	r_8 limber3 =  limb3Norm * facts[n]
	  * ( limb3C0*lag0 + limb3C1 * lag1 + limb3C2 * lag2 + limb3C3 * lag3 );
 
	
273
        JlnpApprox = limber3; 
274 275


276 277 278 279 280
	int is_limberApprox = 1;
	// 3) if the two Limber aprox are different => compute by CC integration
	if(fabs(limber0-limber3)>tol4Limber){
	is_limberApprox = 0;

281
#endif
282 283 284 285 286 287
	  
	  IntFunc1D f(l,klptau,Ln);

	  //upper bound of integration is the Max between the 2nd BesselRoot rescaled and the Laguerre Function "end point" 
	  r_8 intUppBound = std::max(integUpperBound[n], ql2scaled);

288

289 290 291
	  theQuad.SetFuncBounded(f,intLowBound,intUppBound);
	  integ_val = theQuad.GlobalAdapStrat(tol4Integral);
	  JlnpApprox = integ_val.first * (facts[n]*sqrt2Ovpi);
292 293 294 295

#if LIMBER>0
	}//end limber test
#endif
296 297 298

	Jlnp_[ploff7 + n*Jstride] = JlnpApprox;

299

300 301
	ofs <<l<<" "
	    <<n<<" "
302 303 304
 	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] 
	    << endl;

305 306 307 308
	
	Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;

	delete Ln;
309
#if LIMBER>0
310 311 312
	delete Ln_ap1;
	delete Ln_ap2;
	delete Ln_ap3;
313
#endif
314 315 316 317 318 319 320 321 322 323

      }//n
    }//l
  }//p


  ofs.close();

  tstack_pop("J Integ compute ");

324 325
  cout << "(JEC) ComputeJInteg END " << endl;

326

327
}//ComputeJlpnInteg
328 329


330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460

/*!
  Computation of the Jlnp = Jln(k_{lp}) = Sqrt(2/pi) tau^3/2 Int_0^infty jl(k_{lp} tau x) K_n(x) x^2 dx
  with
  K_n(x) = Exp[-x/2] Ln(x) Sqrt(n!/(n+alpha)!) ->cf. LaguerreFuncQuad
  k_{lp}*tau = q_{lp}/r_{Nmax-1}
  Apply a Discrete Spherical transform Lanusse et al A&A 540, A92 (2012) Eq.40 with l'=l
  Jlnp ~ Kapa_l^(-3) tau^(3/2) Sqrt(2 Pi) Sum_{p'=0}^{PPmax-1} Kn(r_{lp'}) jl(q_{lp'}q_{lp}/q_{l,PPmax-1})/[jl+1(q_{lp})]^2
 with r_{lp'} = q_{lp'}/q_{l,PPmax-1}r_{Nmax-1}
 notice PPmax as nothing to do in principle with Pmax the Fourier-Bessel max value of the 3rd index
 Kapa_l = q_{l,PPmax-1}/r_{Nmax-1}

JEC: 12/11/15 Add the contribution of Integral_{r_{Nmax-1}}^{r_{Nmax-1} + 3(r_{Nmax-1}-r_{Nmax-2})} by  Clenshaw-Curtis quadrature
 */


// void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){

//   r_8 tau = lagTrans_->GetTau();
//   r_8 tau3sqrt = pow(tau,(r_8)(3./2.));

// //   cout << "(JEC) tau3sqrt= " << tau3sqrt << endl;
  
//   //JEC 12/11/15 Start
//   r_8 rNm1 = Rmax_ / tau; // R = r[N-1] voir acces direct  = (lagTrans_->GetNodes()).back()
//   r_8 rNm2 = (lagTrans_->GetNodes())[Nmax_-2];
//   r_8 deltaR = 3.0 * (rNm1-rNm2);
//   r_8 integLowerBound = rNm1;
//   r_8 integUpperBound = rNm1+deltaR; //for the integral 
//   r_8 tol = 1e-10;
//   typedef  ClenshawCurtisQuadrature<double> Integrator_t;
//   string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-20.txt";
//   Integrator_t theQuad(20,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
  
//   Quadrature<r_8,Integrator_t>::values_t integ_val;
//   //JEC 12/11/15 end
  
//   int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
//   Jlnp_.resize(Nmax_*Jstride);  

//   //JEC 16/11/15 START: Todo all that might be in common somewhere (hint: laguerreTransform.h)
//   int alpha = lagTrans_->GetAlpha();
//   r_8 alphaFact = 1; // alpha!
//   for(int i=1;i<=alpha; i++) alphaFact *= i;
//   alphaFact  = sqrt((r_8)alphaFact); //sqrt(alpha!)
//   r_8 invalphaFact = 1.0/alphaFact; // 1/sqrt(alpha)!                   
  
//   vector<r_8> facts(Nmax_); //sqrt(n!/(n+alpha)!) en iteratif
//   facts[0] = invalphaFact;
//   for(int n=1;n<Nmax_;n++) {
//     facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) );
//   }
//   //END
  
    
//   int PPmaxm1 = PPmax_ -1;

//   tstack_push("J Integ compute ");

//   stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<"-PP"<<PPmax_<<".txt";
//   std::ofstream ofs;
//   string fname(jlnpDir+"/"+ss.str());
//   ofs.open (fname.c_str(), std::ofstream::out);


// //   cout << "(JEC) Jlpn recompute..." << endl;
//   for(int p=0; p<Pmax_; p++){
//     for(int l=0; l<Lmax_; l++){
      
//       int ploff7 = p + l*Pmax_;
//       r_8 qlp = (*jnu_)(l,p);                //ql,p
      
//       r_8 qlpx = (*jnu_)(l,PPmaxm1);         //ql,PPmax-1
      
//       r_8 Kapa_l = qlpx/rNm1;                  // Kapa_l = (ql,PPmax-1)/R
//       r_8 Kapa_l3 =  Kapa_l *  Kapa_l *  Kapa_l;
      
//       r_8 coeff = sqrt(2.0 *M_PI)/Kapa_l3;  // Sqrt(2 Pi)/Kapa_l^3
      
//       for (int n=0; n<Nmax_; n++){
// 	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
	
// 	r_8 sum = 0.;
// 	//contribution of the integral: [0,rNm1] estimated
// 	//using Discrete Fourier-Bessel Transform

// 	for (int pp=0; pp<PPmax_; pp++){
// 	  r_8 qlpp = (*jnu_)(l,pp);           //ql,pp
// 	  r_8 jlp1 = Bessel::jn(l+1,qlpp);       //j(l+1,ql,pp)
	  
// 	  r_8 rlpp = qlpp/qlpx * rNm1;          // rl,pp = (ql,pp/ql,PPmax-1) * rNm1
// 	  r_8 jlqq = Bessel::jn(l,qlpp * qlp/qlpx);
	  
// 	  sum +=  jlqq/(jlp1*jlp1) * Ln->Value(rlpp);

// 	}//pp-loop
	

// 	//	Jlnp_[ploff7 + n*Jstride] = (coeff * sum) * (tau3sqrt*facts[n]);
// 	//JEC 23/11/15 save the Rmax-independant part
// 	Jlnp_[ploff7 + n*Jstride] = coeff * facts[n] * sum;
	
// 	//contribution [rNm1, rNm1+deltaR] estimated with Clenshaw-Curtis quadrature (JEC 12/11/15) Start
// 	IntFunc1D f(l,qlp/rNm1,Ln);
// 	theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
// 	integ_val = theQuad.GlobalAdapStrat(tol,2,1);
// 	//	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (tau3sqrt*facts[n]*sqrt(2.0/M_PI));
// 	//JEC 23/11/15 save the Rmax-independant part
// 	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (facts[n]*sqrt(2.0/M_PI));
	
		
// 	//JEC 23/11/15 save the Rmax-independant part
// 	ofs <<l<<" "
// 	    <<n<<" "
// 	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
		
// 	Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;

// 	delete Ln;
	
//       }//n
//     }//l
//   }//p
  
//   ofs.close();

//   tstack_pop("J Integ compute ");
  
// }//ComputeJlpnInteg


461
/*!
462 463
  FBlmp = FB_{lm}(k_{lp}) = Sum_n FL_{lmn} J_{ln}(k_{lp}) = Sum_n FL_{lmn} J_{lnp}
  the J_{lnp} includes the sqrt(2/pi) contrary to McEwen paper Exact Wavelet on the Ball
464 465 466
  This is the Fourrier-Bessel coefficients related to McEwen & Leistedh definition
  In case one wants to use Lanusse & Starck definition then
  FBlmp_Lanusse = k_{lp} * FBlmp_McEwen
467
*/
468 469 470
void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
				    vector< complex<r_8> >& FBlmp){

471 472
  
  int Jstride = Lmax_*Pmax_;
473

474
#if DEBUG>=1
475 476
  std::ofstream ofs;
  ofs.open ("FBlmp.txt", std::ofstream::out);
477
#endif
478

479 480 481
  //try here a pedestrian computation
  for(int p=0; p<Pmax_; p++){
    for(int l=0; l<Lmax_; l++){
482

483
      int ploff7 = p + l*Pmax_;
484

485 486
      for (int m=0; m<=l; m++) {
	int almoff7 = l+m*Lmax_-m*(m+1)/2;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
487
	int idBess = p*Nalm_ + almoff7;
488
	complex<r_8> sum = 0.0;
489

490
	for (int n=0; n<Nmax_; n++){
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
491
	  int idLag =  n*Nalm_ + almoff7;
492
	  sum += FLlmn[idLag]*Jlnp_[ploff7 + n*Jstride];	  
493
	}//n
494
	FBlmp[idBess] = sum;
495

496
#if DEBUG>=1
497 498 499
       	ofs <<l<<" "
	    <<m<<" "
	    <<p<< " " << setprecision(12) << 	FBlmp[idBess] << endl;
500
#endif
501

502 503 504
      }//m
    }//l
  }//p
505

506
#if DEBUG>=1
507
  ofs.close();
508
#endif
509 510
  

511 512 513
}//Transform


Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
514 515
void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp, 
				vector< complex<r_8> >& FBalmk,
516
				vector<r_8>& fijk
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
517 518 519 520
				) {

  int NtotCoeffFB = Nalm_*Pmax_; 
  if(FBlmp.size() != (size_t)NtotCoeffFB) throw LagSHTError("Laguerre2Bessel::Synthesis size error");
521 522

  cout << "Dump  Laguerre2Bessel::Synthesis START "<<endl;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
523 524 525 526 527 528
    

  int NtotCoeffFL = Nalm_*Nmax_; 
  FBalmk.resize(NtotCoeffFL); 
  
  /*
529 530
    This is the McEwen & Leistedh definition of the Fourier-Bessel Transform
    FBalmk = alm^{FB}(R_k) = Sum_{p=0}^{Pmax-1} Kapa_{lp} FBlmp  j_l(k_{lp} R_k)
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
531 532 533
    with
    kapa_{lp} = Sqrt(2pi) Rmax^{-3}/j_{l+1}^2(q_{lp})
    k_{lp} = q_{lp}/Rmax
534
    R_k = tau r_k
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
535
    r_k: k-th nodes of Laguerre of order Nmax as for the Spherical-Laguerre Transform.
536
    tau: Rmax/r_{N-1} from Laguerre-Transform
537
    => k_{lp} R_k = q_{lp} r_k/r_{N-1}
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
538 539

    Rdeminder: Fourier-Laguerre case where
540
    alm^{FL}(R_k) = Sum_{n=0}^{Nmax-1} FLlmn LaguerreFunc_n(r_k)
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
541 542 543 544 545 546 547
    
    in the Fourier-Bessel case the element of the sum, except the FBlmp, are depending
    on "l" also throw kapa_{lp}, k_{lp} and j_l(k_{lp}.) 
   
   */
  
  vector<r_8> vrk = lagTrans_->GetNodes();
548 549
  r_8 rNm1 = vrk[Nmax_-1];

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
550 551 552
  //pedestrian computation
  r_8 fact = sqrt(2.*M_PI)/(Rmax_*Rmax_*Rmax_);

553
  vector< complex<r_8> > sumElmt(Pmax_);
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
554
  for(int k=0;k<Nmax_;k++){
555 556 557

    r_8 rk_red = vrk[k]/rNm1;

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
558
    for(int l=0;l<Lmax_;l++){
559
      for(int m=0;m<=l;m++){ 
560

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
561 562
	int lmOff7 = l+m*Lmax_-m*(m+1)/2;
	int id0 = k*Nalm_ + lmOff7;
563

564 565 566 567 568
	for(int p=0;p<Pmax_;p++){ //Fill
 	  r_8 qlp = (*jnu_)(l,p);
 	  r_8 jlp1 = Bessel::jn(l+1,qlp);	  
 	  r_8 kapalp = fact/(jlp1*jlp1);
 	  int id1 = p*Nalm_ + lmOff7;
569

570 571
  	  sumElmt[p] = kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
	}
572
	//kahan summation	FBalmk[id0] = LagSHT::accumulate(sumElmt.begin(),sumElmt.end());
573
	FBalmk[id0] = std::accumulate(sumElmt.begin(),sumElmt.end(),complex<r_8>(0.,0.));
574

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
575 576 577 578
      }//loop on m
    }//loop on l
  }//loop on k

579

580 581
  //from Almk to Fijk with libsharp

582
  int Ntot3DPix = Npix_*Nmax_;
583
  fijk.resize(Ntot3DPix); 
584

585 586 587 588 589
   int off7n, off7k;
   //perform the synthesis per shell
   for(int n =0; n<Nmax_; n++){
     off7n = n*Nalm_;
     off7k = n*Npix_; 
590

591
     sphere_->alm2map(reinterpret_cast<const r_8*>(&(FBalmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
592

593
   }//loop on shell
594 595
  

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
596 597 598
}//Synthesis


599 600 601
void Laguerre2Bessel::Print(ostream& os, int what){
  if(what>0){

602
    //Debug the Jlnp coefficients
603
    int stride = Lmax_*Pmax_;
604 605 606
    for (int pc = 0; pc < Pmax_; pc++) {
      for (int lc = 0; lc < Lmax_; lc++) {
	int ploff7 = pc + lc*Pmax_;
607
	for (int nc = 0; nc < Nmax_; nc++) {
608
	  os << lc<<" "<<nc<<" "<<pc<<" " << setprecision(12)<< Jlnp_[ploff7 + stride*nc] << endl;
609 610 611 612 613 614 615
	}//loop nc
      }//loop pc
    }//loop lc
  }//what
}//Print

}//namespace