(JEC) 1/2/16 Jlnp compute by Clenshaw-Curtis quadrature only and bounds control

src/lagsht_testsuite.cc

@@ -453,6 +453,46 @@ void SpinMultiSphericalLaguerreTransform() {
 
 }//SpinMultiSphericalLaguerreTransform
 
+//-------------------------------------------------------------
+void TestJlnpComputation() {
+
+  int Nmax = param.N;
+  int Lmax = param.Lmax;
+  int Pmax = param.Pmax;
+  r_8 Rmax = param.R;
+  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; //       "
+
+  cout << " ______________ TestJlnpComputation  START ___________ " << endl;
+  LaguerreSphericalTransform sphlagtrans(geometry,
+					 Lmax,
+					 -1,
+					 Nmax,
+					 Rmax,
+					 ntheta,
+					 nphi
+					 );
+
+  
+  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
+  LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
+
+
+  tstack_push("Ctor lag2bess");
+  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
+			   clenshawDir,jlnpDir,recomputeJlnp);
+
+  tstack_pop("Ctor lag2bess");
+
+  cout << " ______________ TestJlnpComputation  End  ___________ " << endl;
+
+}
+
+
 //-------------------------------------------------------------
 void TestLaguerre2Bessel() {
 
@@ -794,7 +834,7 @@ int main(int narg, char *arg[]) {
   char* pLAGSHTROOTDIR; pLAGSHTROOTDIR = getenv("LAGSHTROOTDIR");
   string clenshawDir = string(pLAGSHTROOTDIR)+"/data";
   string jlnpDir = string(pLAGSHTROOTDIR)+"/data";
-  bool recomputeJlnp = false;
+  bool recomputeJlnp = true;
 
 
   //JEC 12/1/16 alpha as input
@@ -944,7 +984,8 @@ int main(int narg, char *arg[]) {
 	}
 	if(Pmax<N){ param.Pmax = 16*N; cout << "(Warning): Pmax<Nmax => modify the value to "<<  param.Pmax << endl;}
 	tstack_push("TestLaguerre2Bessel");
-	TestLaguerre2Bessel();
+	//	TestLaguerre2Bessel();
+	TestJlnpComputation();
  	tstack_pop("TestLaguerre2Bessel");
  	tstack_report("TestLaguerre2Bessel");
       }

src/laguerre2bessel.cc

@@ -20,6 +20,8 @@ namespace LagSHT {
      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,
@@ -29,7 +31,11 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
 				 bool recomputeJlnp
 				 ):
   sphere_(sphere),lagTrans_(lagTrans),
-  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax), PPmax_(Pmax) {
+  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
@@ -41,12 +47,13 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
     throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
   }
   
-  jnu_ = new Bessel(Lmax_,max(Pmax_,PPmax_));  //the max is here just in case PPmax =/= Pmax
+  //  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
   
   if(!recomputeJlnp){
     
     //read the Rmax-independant values of the jlnp and rescale
-    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<"-PP"<<PPmax_<<".txt";
+    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
     
     string fname(jlnpDir+"/"+ss.str());
     std::ifstream ifs(fname.c_str(), std::ifstream::in);
@@ -94,7 +101,6 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
 
  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();
@@ -105,9 +111,9 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
   //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 deltaR = (rNm1-rNm2);
   r_8 integLowerBound = rNm1;
-  r_8 integUpperBound = rNm1+deltaR; //for the integral 
+  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";
@@ -150,6 +156,9 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
       
       int ploff7 = p + l*Pmax_;
       r_8 qlp = (*jnu_)(l,p);                //ql,p
+
+
+      r_8 klp = qlp/rNm1; //JEC 25/1/16
       
       r_8 qlpx = (*jnu_)(l,PPmaxm1);         //ql,PPmax-1
       
@@ -160,6 +169,9 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
       
       for (int n=0; n<Nmax_; n++){
 	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
+
+
+	
 	
 	r_8 sum = 0.;
 	//contribution of the integral: [0,rNm1] estimated
@@ -189,6 +201,23 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
 	//JEC 23/11/15 save the Rmax-independant part
 	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (facts[n]*sqrt(2.0/M_PI));
 	
+
+	//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;
+		
 		
 	//JEC 23/11/15 save the Rmax-independant part
 	ofs <<l<<" "
@@ -207,6 +236,109 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
 
   tstack_pop("J Integ compute ");
   
+}//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 ");
+
+
 }//ComputeJlpnInteg
 
 

src/laguerre2bessel.h

@@ -153,7 +153,7 @@ protected:
   int Pmax_; //!< Max. order of the Bessel coeff.
   r_8 Rmax_; //!< Radial maximal value  
 
-  int PPmax_; //!< cut off in the alm computation
+  //JEC 1/2/16  int PPmax_; //!< cut off in the alm computation
 
   
   Bessel* jnu_; //!< Spherical Bessel fnt (owner of the ptr)

src/quadinteg.h

@@ -36,12 +36,11 @@
 #include <fstream>
 #include <iomanip>
 #include <string>
+#include <algorithm>
 
 #include <limits> //for precision features
 #include <iterator> //ostream_iterator
 
-#include "walltimer.h"    //timing
-
 
 #define DUMPLEV 0
 
@@ -299,7 +298,6 @@ public:
       integral += absw_[i] * fi;
       error    += errw_[i] * fi;
     }
-    
     return ResultClass(integral*boundsDist,fabs((double)error*boundsDist),xi);
   }//IntEval
   
@@ -393,7 +391,9 @@ public:
 
     while( ((error >= tol*fabs(integral)) && (n<MaxNumRecurssion)) || 
 	   (n<MinNumRecurssion) ) {
-      n++;
+      n++; 
+
+      //      cout << "Recursion: " << n << "int: " << integral << " error/int: " << error/fabs(integral) << endl;
 
       curBounds = theResults.front().first; //{x0,x2}
       T xmiddle = (curBounds.first+curBounds.second)*0.5; //x1
@@ -432,10 +432,24 @@ public:
     }//eod while
      
      //debug
-/*     cout << "GlobalAdapStrat end: {" << n << ", " << integral << ", "  */
-/*     	 << error << "}" << endl; */
-
-    
+/*     cout << "GlobalAdapStrat end: {" << n << ", " << integral << ", " 
+       << error << "}" << endl;
+       cout << "debug quadinteg: " << endl;
+       cout << "{";
+       vector<r_8> pts;
+       for(iRes=theResults.begin();iRes!=iEnd;++iRes){
+       result_t iResCur=*iRes;
+       for(size_t i=0; i<npts_;i++){
+       bounds_t bnd = iResCur.first;
+       pts.push_back(bnd.first + absc_[i]*(bnd.second-bnd.first));
+       }
+       }
+       std::sort(pts.begin(),pts.end());
+       for(int i=0;i<pts.size();i++){
+       cout << "{" << pts[i] <<"," << (*func_)(pts[i])<<"},"; 
+       }
+       cout << "};"<<endl;
+*/
 
     return make_pair(integral,error);
   }//GlobalAdapStrat