laguerre2bessel.cc 14.8 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

16 17
namespace LagSHT {

18 19 20 21
  /* 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)
     Version 1: on fait une integration a la Clenshaw-Curtis avec startegie Localeo ou Globale
     recursive
22
     Version 2: Jlnp calcule via une Discrete Fourier-Bessel Transform and a  Clenshaw-Curtis quadrature and global integration strategy. 
23 24
     Version 3: calcul via une quadrature de  Clenshaw-Curtis quadrature and global integration en controlant les bornes
                d'integration
25
  */
26 27 28 29 30 31 32 33
Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
				 LaguerreTransform* lagTrans,
				 int Nmax, int Pmax, r_8 Rmax,
				 string clenshawDir,
				 string jlnpDir,
				 bool recomputeJlnp
				 ):
  sphere_(sphere),lagTrans_(lagTrans),
34 35 36 37 38
  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
  //JEC 1/2/16  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmaax), PPmax_(Pmax) {

  //  PPmax_ = 4096;

39 40 41 42 43 44 45 46 47 48 49
  
  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");
  }
  
50 51
  //  jnu_ = new Bessel(Lmax_,std::max(std::max(Pmax_,PPmax_),Nmax_));  //the max is here just in case PPmax =/= Pmax
  jnu_ = new Bessel(Lmax_,std::max(Pmax_,Nmax_));  //the max is here just in case Nmax =/= Pmax
52 53 54 55
  
  if(!recomputeJlnp){
    
    //read the Rmax-independant values of the jlnp and rescale
56
    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
57
    
58 59
    string fname(jlnpDir+"/"+ss.str());
    std::ifstream ifs(fname.c_str(), std::ifstream::in);
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
    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
// 	  if( (l!=lcur) || (n!=ncur) || (p!=pcur))                         //check
// 	    throw LagSHTError("ERROR: jlnp file read failed.:"+ss.str());
	  Jlnp_[ploff7 + n*Jstride] =  tau3sqrt * val;                     //rescaling
	}//n
      }//l
    }//p
    
  } else {
    ComputeJInteg(clenshawDir,jlnpDir);
  }//jlnp omputation/load
  
86 87
}//Ctor

88

89 90 91 92 93 94 95 96 97 98
/*!
  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}
99 100 101 102 103



 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

104
void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
105 106 107

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

//   cout << "(JEC) tau3sqrt= " << tau3sqrt << endl;
110
  
111 112 113
  //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];
114
  r_8 deltaR = (rNm1-rNm2);
115
  r_8 integLowerBound = rNm1;
116
  r_8 integUpperBound = rNm1+3.0*deltaR; //for the integral 
117 118
  r_8 tol = 1e-10;
  typedef  ClenshawCurtisQuadrature<double> Integrator_t;
119 120
  string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-20.txt";
  Integrator_t theQuad(20,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
121 122 123
  
  Quadrature<r_8,Integrator_t>::values_t integ_val;
  //JEC 12/11/15 end
124 125
  
  int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
126
  Jlnp_.resize(Nmax_*Jstride);  
127 128

  //JEC 16/11/15 START: Todo all that might be in common somewhere (hint: laguerreTransform.h)
129 130 131 132 133 134
  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)!                   
  
135
  vector<r_8> facts(Nmax_); //sqrt(n!/(n+alpha)!) en iteratif
136 137 138 139
  facts[0] = invalphaFact;
  for(int n=1;n<Nmax_;n++) {
    facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) );
  }
140
  //END
141
  
142
    
143
  int PPmaxm1 = PPmax_ -1;
144

145
  tstack_push("J Integ compute ");
146

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

152 153

//   cout << "(JEC) Jlpn recompute..." << endl;
154 155 156 157 158
  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
159 160 161


      r_8 klp = qlp/rNm1; //JEC 25/1/16
162 163 164
      
      r_8 qlpx = (*jnu_)(l,PPmaxm1);         //ql,PPmax-1
      
165
      r_8 Kapa_l = qlpx/rNm1;                  // Kapa_l = (ql,PPmax-1)/R
166 167
      r_8 Kapa_l3 =  Kapa_l *  Kapa_l *  Kapa_l;
      
168
      r_8 coeff = sqrt(2.0 *M_PI)/Kapa_l3;  // Sqrt(2 Pi)/Kapa_l^3
169 170 171
      
      for (int n=0; n<Nmax_; n++){
	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
172 173 174


	
175 176
	
	r_8 sum = 0.;
177 178
	//contribution of the integral: [0,rNm1] estimated
	//using Discrete Fourier-Bessel Transform
179

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

189
	}//pp-loop
190 191
	

192 193 194
	//	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;
195 196 197 198 199
	
	//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);
200 201 202
	//	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));
203
	
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220

	//Limber approx JEC 25/1/16 the sqrt(2/Pi) is canceled by the definition of SphericalBesselJ function used in this approx
	r_8 limber = pow(klp,(r_8)(-3.0)) * pow((r_8)(l+0.5),(r_8)1.5) * Ln->Value(((r_8)(l)+0.5)/klp)*facts[n];

	//New Approx JEC 27/1/16
	theQuad.SetFuncBounded(f,0.8*r_8(l),rNm1+4.0*deltaR);
	integ_val = theQuad.GlobalAdapStrat(tol,2,1);
	

	cout <<l<<" "
	     <<n<<" "
	     <<p<< " "
	     << setprecision(30) << Jlnp_[ploff7 + n*Jstride]  
	     << " " << limber
	     << " " << integ_val.first *  (facts[n]*sqrt(2.0/M_PI))
	     << endl;
		
221
		
222
	//JEC 23/11/15 save the Rmax-independant part
223 224
	ofs <<l<<" "
	    <<n<<" "
225 226 227 228
	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
		
	Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;

229 230 231 232 233 234
	delete Ln;
	
      }//n
    }//l
  }//p
  
235
  ofs.close();
236

237 238
  tstack_pop("J Integ compute ");
  
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
}//ComputeJlpnInteg
*/

/*!
  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}
  
  Use a Clenshaw-Curtis quadrature with 2*40-1 pts between LowBound and UppBound
  - LowBound is (80% * l_value)
  - 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
      - q(l,1)/(k_{lp}*tau) this is an estimate of first oscillation of Bessel function
*/
void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){


  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()
  r_8 tol = 1e-10;
  
  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;
  int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
  Jlnp_.resize(Nmax_*Jstride);  

  r_8 alpha = lagTrans_->GetAlpha();
  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)) );

 
  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
  }


  stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
  std::ofstream ofs;
  string fname(jlnpDir+"/"+ss.str());
  ofs.open (fname.c_str(), std::ofstream::out);
  
  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 klptau = qlp/rNm1;
      r_8 intLowBound = 0.8*l;

      for (int n=0; n<Nmax_; n++){


	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
	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], (*jnu_)(l,1)/klptau);
	

	theQuad.SetFuncBounded(f,intLowBound,intUppBound);
	integ_val = theQuad.GlobalAdapStrat(tol);
	Jlnp_[ploff7 + n*Jstride] = integ_val.first * (facts[n]*sqrt(2.0/M_PI));

 	cout <<l<<" "
	     <<n<<" "
	     <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
	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 ");


342
}//ComputeJlpnInteg
343 344


345
/*!
346 347
  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
348 349 350
  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
351
*/
352 353 354
void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
				    vector< complex<r_8> >& FBlmp){

355 356
  
  int Jstride = Lmax_*Pmax_;
357

358
#if DEBUG>=1
359 360
  std::ofstream ofs;
  ofs.open ("FBlmp.txt", std::ofstream::out);
361
#endif
362

363 364 365
  //try here a pedestrian computation
  for(int p=0; p<Pmax_; p++){
    for(int l=0; l<Lmax_; l++){
366

367
      int ploff7 = p + l*Pmax_;
368

369 370
      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
371
	int idBess = p*Nalm_ + almoff7;
372
	complex<r_8> sum = 0.0;
373

374
	for (int n=0; n<Nmax_; n++){
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
375
	  int idLag =  n*Nalm_ + almoff7;
376
	  sum += FLlmn[idLag]*Jlnp_[ploff7 + n*Jstride];	  
377
	}//n
378
	FBlmp[idBess] = sum;
379

380
#if DEBUG>=1
381 382 383
       	ofs <<l<<" "
	    <<m<<" "
	    <<p<< " " << setprecision(12) << 	FBlmp[idBess] << endl;
384
#endif
385

386 387 388
      }//m
    }//l
  }//p
389

390
#if DEBUG>=1
391
  ofs.close();
392
#endif
393 394
  

395 396 397
}//Transform


Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
398 399
void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp, 
				vector< complex<r_8> >& FBalmk,
400
				vector<r_8>& fijk
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
401 402 403 404
				) {

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

  cout << "Dump  Laguerre2Bessel::Synthesis START "<<endl;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
407 408 409 410 411 412
    

  int NtotCoeffFL = Nalm_*Nmax_; 
  FBalmk.resize(NtotCoeffFL); 
  
  /*
413 414
    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
415 416 417
    with
    kapa_{lp} = Sqrt(2pi) Rmax^{-3}/j_{l+1}^2(q_{lp})
    k_{lp} = q_{lp}/Rmax
418
    R_k = tau r_k
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
419
    r_k: k-th nodes of Laguerre of order Nmax as for the Spherical-Laguerre Transform.
420
    tau: Rmax/r_{N-1} from Laguerre-Transform
421
    => k_{lp} R_k = q_{lp} r_k/r_{N-1}
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
422 423

    Rdeminder: Fourier-Laguerre case where
424
    alm^{FL}(R_k) = Sum_{n=0}^{Nmax-1} FLlmn LaguerreFunc_n(r_k)
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
425 426 427 428 429 430 431
    
    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();
432 433
  r_8 rNm1 = vrk[Nmax_-1];

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
434 435 436
  //pedestrian computation
  r_8 fact = sqrt(2.*M_PI)/(Rmax_*Rmax_*Rmax_);

437
  vector< complex<r_8> > sumElmt(Pmax_);
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
438
  for(int k=0;k<Nmax_;k++){
439 440 441

    r_8 rk_red = vrk[k]/rNm1;

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
442
    for(int l=0;l<Lmax_;l++){
443
      for(int m=0;m<=l;m++){ 
444

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

448 449 450 451 452
	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;
453

454 455
  	  sumElmt[p] = kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
	}
456
	//kahan summation	FBalmk[id0] = LagSHT::accumulate(sumElmt.begin(),sumElmt.end());
457
	FBalmk[id0] = std::accumulate(sumElmt.begin(),sumElmt.end(),complex<r_8>(0.,0.));
458

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
459 460 461 462
      }//loop on m
    }//loop on l
  }//loop on k

463

464 465
  //from Almk to Fijk with libsharp

466
  int Ntot3DPix = Npix_*Nmax_;
467
  fijk.resize(Ntot3DPix); 
468

469 470 471 472 473
   int off7n, off7k;
   //perform the synthesis per shell
   for(int n =0; n<Nmax_; n++){
     off7n = n*Nalm_;
     off7k = n*Npix_; 
474

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

477
   }//loop on shell
478 479
  

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
480 481 482
}//Synthesis


483 484 485
void Laguerre2Bessel::Print(ostream& os, int what){
  if(what>0){

486
    //Debug the Jlnp coefficients
487
    int stride = Lmax_*Pmax_;
488 489 490
    for (int pc = 0; pc < Pmax_; pc++) {
      for (int lc = 0; lc < Lmax_; lc++) {
	int ploff7 = pc + lc*Pmax_;
491
	for (int nc = 0; nc < Nmax_; nc++) {
492
	  os << lc<<" "<<nc<<" "<<pc<<" " << setprecision(12)<< Jlnp_[ploff7 + stride*nc] << endl;
493 494 495 496 497 498 499
	}//loop nc
      }//loop pc
    }//loop lc
  }//what
}//Print

}//namespace