(JEC) 23/2/16 Laguerre->Bessel using CC quadrature only

src/lagsht_testsuite.cc

@@ -19,7 +19,7 @@ using namespace std;
 
 using namespace LagSHT;
 
-#define DEBUG 0
+#define DEBUG 10
 
 
 //JEC 12/1/16: alpha becomes double instead of int
@@ -463,9 +463,9 @@ void TestJlnpComputation() {
   string geometry = param.geometry;
   int ntheta = param.ntheta;
   int nphi = param.nphi;
-  string clenshawDir = param.clenshawDir;   //JEC 23/11/15
-  string jlnpDir = param.jlnpDir;           //       "
-  bool recomputeJlnp = param.recomputeJlnp; //       "
+  string clenshawDir = param.clenshawDir;
+  string jlnpDir = param.jlnpDir;
+  bool recomputeJlnp = param.recomputeJlnp;
 
   cout << " ______________ TestJlnpComputation  START ___________ " << endl;
   LaguerreSphericalTransform sphlagtrans(geometry,
@@ -483,11 +483,14 @@ void TestJlnpComputation() {
 
 
   tstack_push("Ctor lag2bess");
+  recomputeJlnp = true;
   Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
 			   clenshawDir,jlnpDir,recomputeJlnp);
-
   tstack_pop("Ctor lag2bess");
 
+
+
+
   cout << " ______________ TestJlnpComputation  End  ___________ " << endl;
 
 }
@@ -669,49 +672,48 @@ void TestLaguerre2Bessel() {
 
   
   
-  {//Check 1
-    cout << "Dump FL or FB reconstructed Alm(r_k)" <<endl;
+  {//Check 1b
+    cout << "Dump FL or FB reconstructed Cl(r_k)" <<endl;
 
 #if DEBUG >=2 
     std::ofstream ofs;
-    ofs.open ("Almi.txt", std::ofstream::out);
+    ofs.open ("Cli.txt", std::ofstream::out);
 #endif
     
     for(int n=0;n<Nmax;n++){
-      r_8 err_abs(0.);
-      r_8 err_rel(0.);
-
       for(int l=0;l<Lmax;l++){
+
+	r_8 ClFL = 0;
+	r_8 ClFB = 0;
+
 	for (int m=0;m<=l;m++) {
 	  int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
 
+	  r_8 almkFL2 = FLalmk[id].real()*FLalmk[id].real() + FLalmk[id].imag()*FLalmk[id].imag();
+	  r_8 almkFB2 = FBalmk[id].real()*FBalmk[id].real() + FBalmk[id].imag()*FBalmk[id].imag();
+
+	  ClFL += (m==0) ? almkFL2 : 2*almkFL2;
+	  ClFB += (m==0) ? almkFB2 : 2*almkFB2;
+
+	}//m
+
+	ClFL /= (r_8)(2*l+1);
+	ClFB /= (r_8)(2*l+1);
 #if DEBUG >=2 
-	  ofs <<l<<" "
-	      <<m<<" "
-	      <<n<< " " << setprecision(12) << FLalmk[id].real() << " " <<   FLalmk[id].imag()
-	      << "  " << FBalmk[id].real() << " " << FBalmk[id].imag()  << endl;
+	  ofs <<l<<" " << n << " " << setprecision(12) << ClFL << " " << ClFB << endl;
 #endif	  
-	  complex<r_8> cdiff = (FLalmk[id] -  FBalmk[id])*conj(FLalmk[id] -  FBalmk[id]);
-	  r_8 diff = sqrt(cdiff.real());
-	  if(diff>err_abs){
-	    err_abs = diff;
-	  }
-	  complex<r_8> cref= FLalmk[id]*conj(FLalmk[id]);
-	  r_8 ref = sqrt(cref.real());
-	  r_8 relatdiff = diff/ref;
-	  if(relatdiff>err_rel) err_rel = relatdiff;
-	}
-      }
 
-      cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<n<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
-      cout << "Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
-      cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
-
-    }
+      }//l
+    }//n
 #if DEBUG >=2 
     ofs.close();
 #endif
-  }//check 1
+  }//check 1b
+
+
+
+
+
 
   {//check 2
 #if DEBUG >=2 

src/laguerre2bessel.cc

@@ -17,11 +17,12 @@ namespace LagSHT {
 
   /* 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 Locale ou Globale
+     Version 1: on fait une integration a la Clenshaw-Curtis avec startegie Localeo ou Globale
      recursive
      Version 2: Jlnp calcule via une Discrete Fourier-Bessel Transform and a  Clenshaw-Curtis quadrature and global integration strategy. 
      Version 3: calcul via une quadrature de  Clenshaw-Curtis quadrature and global integration en controlant les bornes
                 d'integration
+
   */
 Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
 				 LaguerreTransform* lagTrans,
@@ -34,9 +35,6 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
   Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
   //JEC 1/2/16  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmaax), PPmax_(Pmax) {
   
-  //  PPmax_ = 4096;
-
-  
   if(sphere_) {
     //store localy usefull parameters for the algorithms Analysis/Synthesis
     Npix_ = sphere_->Npix();
@@ -48,12 +46,21 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
   }
   
   //  jnu_ = new Bessel(Lmax_,std::max(std::max(Pmax_,PPmax_),Nmax_));  //the max is here just in case PPmax =/= Pmax
+  
+  cout << "(JEC) lag2bes Ctor  Compute jnu_ START" << endl;
+
   jnu_ = new Bessel(Lmax_,std::max(Pmax_,Nmax_));  //the max is here just in case Nmax =/= Pmax
   
+  cout << "(JEC) lag2bes Ctor  Compute jnu_ END" << endl;
+
+  
   if(!recomputeJlnp){
     
+    cout << "READ mode......." << endl; 
+
     //read the Rmax-independant values of the jlnp and rescale
-    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
+    //    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
+    stringstream ss; ss << "limber0-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
     
     string fname(jlnpDir+"/"+ss.str());
     std::ifstream ifs(fname.c_str(), std::ifstream::in);
@@ -72,14 +79,17 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
 	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());
+ 	  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 {
+
+    cout << "Compute mode......." << endl; 
+
     ComputeJInteg(clenshawDir,jlnpDir);
   }//jlnp omputation/load
   
@@ -87,20 +97,36 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
 
 
 
+
 /*!
   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)
+  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
+      - 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
+
 */
+
 void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
 
+  cout << "(JEC) ComputeJInteg START " << endl;
+
+
+  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
+
 
   if(Nmax_<2)throw LagSHTError("Laguerre2Bessel::ComputeJInteg Nmax<2");
 
@@ -110,21 +136,36 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
   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;
+  r_8 tol4Integral = 1e-10;
+  r_8 tol4Limber   = 1e-6;
 
   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();
+  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;
+
   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)) );
 
+  cout << "(JEC) ComputeJInteg Facts END " << endl;
+
+ 
+  cout << "(JEC) ComputeJInteg  integUpperBound START " << endl;
 
   vector<r_8> integUpperBound(Nmax_);
   integUpperBound[0] = rNm1; //this is the case of L_0^{\alpha}(x) =1 with no root
@@ -137,83 +178,270 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
     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
   }
 
+  cout << "(JEC) ComputeJInteg  integUpperBound END " << endl;
 
-  stringstream ss; ss << "NEW-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
+
+  stringstream ss; ss << "NEWNEW-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
   std::ofstream ofs;
   string fname(jlnpDir+"/"+ss.str());
   ofs.open (fname.c_str(), std::ofstream::out);
 
-//   //JEC TMP
-//   stringstream ssTMP; ssTMP << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
-  
-//   string fnameTMP(jlnpDir+"/"+ss.str());
-//   std::ifstream ifs(fnameTMP.c_str(), std::ifstream::in);
-//   if(!ifs.is_open())
-//     throw LagSHTError("ERROR: jlnp file open failed.:"+ssTMP.str());
-//   //JEC TMP
-
+  cout << "(JEC) ComputeJInteg write on <" <<fname << ">"<< endl;
 
-
-
-//   int pcur, lcur, ncur; //TMP
-//   r_8 val;              //TMP
+  r_8 sqrt2Ovpi= sqrt(2.0/M_PI); 
   
   for(int p=0; p<Pmax_; p++){
+    cout << "p = " << p;
+
     for(int l=0; l<Lmax_; l++){
+      cout << " l = " << l;
 
       int ploff7 = p + l*Pmax_;
+
       r_8 qlp = (*jnu_)(l,p);                //ql,p
-      r_8 klptau = qlp/rNm1;
+      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;
+      
+      
+
       r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau; //JEC 4/2/16 for low l just start at 0   
 
+    
+      r_8 ql2scaled = (*jnu_)(l,1)/klptau;   //  the 2nd BesselRoot rescaled
+
+      r_8 arg = (r8l+0.5)/klptau;            // (L+1/2)/(klp*tau) for Limber development
+      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));
+
       for (int n=0; n<Nmax_; n++){
 	
-// 	ifs >> lcur >> ncur >> pcur >> val;         //TMP read
-// 	if( (l == lcur) || (n == ncur) || (p == pcur)){    //TMP check
-// 	  ofs <<l<<" "
-// 	      <<n<<" "
-// 	      <<p<< " " << setprecision(30) << val << endl;
-// 	  Jlnp_[ploff7 + n*Jstride] =  tau3sqrt * val;
-// 	  continue;
-// 	}
+	cout << " n = " << n << endl;
 
-	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
-	IntFunc1D f(l,klptau,Ln);
+	//1. compute the Limber (Order 0) approximation
+	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
 
+	r_8 limber0 = limb0Norm * Ln->Value(arg)*facts[n];
 
-	//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);
+	
+	//2.  compute a higher order Limber (Order 3) approximation
+	//h2: term of 2nd order; h3: term of 3rd order
+
+	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);
+
+	
+	r_8 limber3 =  limb3Norm * facts[n]
+	  * ( limb3C0*lag0 + limb3C1 * lag1 + limb3C2 * lag2 + limb3C3 * lag3 );
  
 	
+	r_8 JlnpApprox = limber3; 
+
+
+	int is_limberApprox = 1;
+	// 3) if the two Limber aprox are different => compute by CC integration
+	if(fabs(limber0-limber3)>tol4Limber){
+	is_limberApprox = 0;
+
+	  
+	  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);
+
 	  theQuad.SetFuncBounded(f,intLowBound,intUppBound);
-	integ_val = theQuad.GlobalAdapStrat(tol);
-	Jlnp_[ploff7 + n*Jstride] = integ_val.first * (facts[n]*sqrt(2.0/M_PI));
+	  integ_val = theQuad.GlobalAdapStrat(tol4Integral);
+	  JlnpApprox = integ_val.first * (facts[n]*sqrt2Ovpi);
+	}
+
+	Jlnp_[ploff7 + n*Jstride] = JlnpApprox;
 
-//  	cout <<l<<" "
-// 	     <<n<<" "
-// 	     <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
 	ofs <<l<<" "
 	    <<n<<" "
- 	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
+ 	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] 
+	    << endl;
+
 	
 	Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;
 
 	delete Ln;
+	delete Ln_ap1;
+	delete Ln_ap2;
+	delete Ln_ap3;
 
       }//n
     }//l
   }//p
 
 
-//   ifs.close(); //TMP
   ofs.close();
 
   tstack_pop("J Integ compute ");
 
+  cout << "(JEC) ComputeJInteg END " << endl;
+
 
 }//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}
+  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
+
+
 /*!
   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
@@ -369,157 +597,3 @@ void Laguerre2Bessel::Print(ostream& os, int what){
 }//Print
 
 }//namespace
-
-
-/*!
-  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 = (rNm1-rNm2);
-  r_8 integLowerBound = rNm1;
-  r_8 integUpperBound = rNm1+3.0*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);
-