#include <stdlib.h>
#include <iostream>
#include <iomanip>      // std::setprecision
#include <fstream>
#include <sstream>
#include <string.h>
#include <limits>       // std::numeric_limits
#include <string>
#include <typeinfo>
using namespace std;

#include "lagsht_exceptions.h"
#include "lagsht_utils.h" //getMemorySize
#include "walltimer.h"    //timing
#include <omp.h>          //OpenMP for sampling

#include "lagSphericTransform.h"
#include "laguerre2bessel.h"

//Test FFT
#include <vector>
#include <algorithm>
#include <functional>
#include <math.h>
#include <fftw3.h>



using namespace LagSHT;

#define DEBUG 10
//JEC 12/1/16: alpha becomes double instead of int

//-------- Parameters set in the main and  used in the different test functions
struct PARAM {
  int Lmax;
  int N;
  r_8 R;
  int Pmax;
  int spin;
  //  int alpha;
  double alpha;
  string geometry;
  int ntheta;
  int nphi;
  string clenshawDir;
  string jlnpDir;
  bool recomputeJlnp;
} param ; 

//simple random generator (taken from sharp_testsuite
static double drand (double min, double max, int *state){
  *state = (((*state) * 1103515245) + 12345) & 0x7fffffff;
  return min + (max-min)*(*state)/(0x7fffffff+1.0);
}//rand


//----------------------------------------------------
void TestBasic() {

  // to access some general parameters... use param.



  cout << " ___________  BASIC TEST _____________ " << endl;
  

  cout << " sizeof(double): " << sizeof(r_8) << " bytes" 
       << " min " << std::numeric_limits<r_8>::min() << "-ln(min): " << -log(std::numeric_limits<r_8>::min()) 
       << " max " << std::numeric_limits<r_8>::max()  
       << " min exp_10 " << std::numeric_limits<r_8>::min_exponent10
       << " max exp_10 " << std::numeric_limits<r_8>::max_exponent10  
       <<  endl;
  


  cout << " sizeof(long double): " << sizeof(r_16) << " bytes" 
       << " min " << std::numeric_limits<r_16>::min() <<  "-ln(min): " << -logl(std::numeric_limits<r_16>::min()) 
       << " max " << std::numeric_limits<r_16>::max()  
       << " min exp_10 " << std::numeric_limits<r_16>::min_exponent10
       << " max exp_10 " << std::numeric_limits<r_16>::max_exponent10  
       <<  endl;
}//TestBasic

//----------------
void LaguerreQuadrature() {

  int N = param.N;
  //  int alpha = param.alpha;
  double alpha = param.alpha;

  cout << " ___________  LaguerreQuadrature  TEST _____________ " << endl;

  LaguerreFuncQuad lagFunc(N,alpha);
  cout << "N = " << lagFunc.GetOrder() << " alpha = " << lagFunc.GetAlpha() << endl;
  vector<r_8> nodes;
  vector<r_8> weights;
  lagFunc.QuadWeightNodes(nodes,weights);

  {
    char buff1[80], buff2[80];
    ofstream ofs1, ofs2;
    sprintf(buff1,"lagNodes-%d-Func.txt",N);
    sprintf(buff2,"lagWeights-%d-Func.txt",N);
    ofs1.open(buff1,ofstream::out);
    ofs2.open(buff2,ofstream::out);
    for (int i=0;i<N;i++){
      ofs1 << setprecision(20) << nodes[i] << endl;
      ofs2 << setprecision(20) << weights[i] << endl;
    }
  }


}//LaguerreQuadrature

//--------------------------------------
void MultiLaguerreTransformSdtAlone() {
  int N = param.N;
  r_8 R  = param.R;
  r_8 alpha = param.alpha; //JEC 12/1/16 to be tested with alpha =/= 2 (eg. alpha = 0)

  cout << " ___________   MultiLaguerreTransformSdtAlone  TEST _____________ " << endl;

  //    Timer tm("Test 3 (r_16)");
      
  int mapSize=2; //this is the stride
      
  int Ntot = N*mapSize;

#if DEBUG >= 1
  cout << "Start random generation...." << endl;
#endif
  vector< complex<r_8> > fnOrig(Ntot);
  int state=1234567 + 8912 ; //random seed
  for(int i=0;i<Ntot;i++){
    r_8 rv = drand(-1,1,&state);
    r_8 iv = drand(-1,1,&state);
    fnOrig[i] = complex<r_8>(rv,iv);
  }
#if DEBUG >=2 
  std::copy(fnOrig.begin(), fnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
#endif
    
  LaguerreTransform trans(N,R,alpha);  //JEC 12/1/16 add alpha

  vector< complex<r_8> > fi(Ntot);
#if DEBUG >= 1
  cout << "Start MultiSynthesis...." << endl;
#endif
  trans.MultiSynthesis(fnOrig,fi,mapSize);
    
#if DEBUG >=2 
  std::copy(fi.begin(), fi.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
#endif

  vector< complex<r_8> > fn(Ntot);
  trans.MultiAnalysis(fi,fn,mapSize);

  r_8 err_abs(0.);
  r_8 err_rel(0.);
  int imax=-1;
  for(int i=0;i<Ntot;i++){
    if(i<10)cout << "("<<i<<"): " << fnOrig[i] << " <-> " << fi[i] << " <-> "  
		 << fn[i] << endl;
    complex<r_8> cdiff = (fnOrig[i] -  fn[i]) * conj(fnOrig[i] -  fn[i]);
    r_16 diff = sqrt(cdiff.real());
    if(diff>err_abs){ 
      err_abs = diff;
      imax = i;
    }
    complex<r_8> foriCAbs = fnOrig[i]*conj(fnOrig[i]);
    r_8 foriAbs = sqrt(foriCAbs.real());
    r_8 relatdiff = diff/foriAbs;
	
    if((relatdiff)>err_rel) err_rel = relatdiff;
  }
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  cout << "r_16 Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  

}//MultiLaguerreTransformSdtAlone

//---------------------------------------
void MultiSphericalLaguerreTransform() {

  r_8 alpha = param.alpha;  //JEC 13/1/16
  int Nmax = param.N;
  int Lmax = param.Lmax;
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;


  cout << " ___________  MultiSphericalLaguerreTransform   TEST _____________ " << endl;

  tstack_push("data input");

  int state = 1234567 + 8912 ; //random seed

  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi,
					 alpha
					 );


  int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
  int Nshell = sphlagtrans.GetOrder();
  int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
  int Npix =  sphlagtrans.GetSphericalGeometry()->Npix();
  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels

  cout << "Verif: Ntheta, Nphi, Npix, Nptot, Nalm, Nshell, Ntot: "
       << sphlagtrans.GetSphericalGeometry()->NTheta() << " "
       << sphlagtrans.GetSphericalGeometry()->NPhi() << " "
       <<  Npix << " "
       <<  NpTot << " "
       <<  Nalm << " "
       <<  Nshell << " "
       <<  Ntot << endl;

  unsigned int maxmemsize =  getMemorySize()/1e6;
  unsigned int estimateMem = 8*((uint_8)(NpTot+2*Ntot))/1e6;
  cout << "Estimate usage memory: " << estimateMem << " MBytes" << endl;
  if ( estimateMem > (0.9*maxmemsize) ) {
    cout << ">>>>> Warning: estimate Mem " << estimateMem <<" MB > 90\% "<< maxmemsize << " MB" << endl;
  }//test memory

  vector< complex<r_8> > flmnOrig(Ntot);
  int id=-1;
  for(int n=0;n<Nmax;n++){ //on each shell
    for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
      for(int l=m;l<Lmax;l++){
	id++;
	r_8 rv = drand(-1,1,&state);
	r_8 iv = (m==0) ? 0.0 : drand(-1,1,&state);
	flmnOrig[id] = complex<r_8>(rv,iv);	    
      }//end l-loop
    }//end m-loop
  }//end loop on shell

#if DEBUG >=2 
  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
#endif
      
  
  tstack_pop("data input");
  tstack_push("processing");
  tstack_push("processing part Synthesis");

#if DEBUG >= 1
  cout << "Main:Synthesis r_8 function..." << endl;
#endif
  vector<r_8> fijk(NpTot);
  sphlagtrans.Synthesis(flmnOrig,fijk);
#if DEBUG >= 2
  for (int i=0; i<NpTot; i++){
    cout << "fijk("<<i<<"): " << fijk[i] << endl;
  }
#endif
  
  tstack_pop("processing part Synthesis");
  tstack_push("processing part Analysis");
  
#if DEBUG >= 1
  cout << "Main:Analysis r_8 function..." << endl;
#endif
  vector< complex<r_8> > flmn(Ntot);
  sphlagtrans.Analysis(fijk,flmn);
    
  tstack_pop("processing part Analysis");
  tstack_pop("processing");  



  cout << "Error analysis..." << endl;

  r_8 err_abs(0.);
  r_8 err_rel(0.);
  int imax = -1;
  for(int i=0;i<Ntot;i++){
    if(i>(Ntot-9)) cout << "(" << i << ") : flnm Orig " <<  flmnOrig[i] << " <-> flmn Rec " <<  flmn[i] << endl;

    complex<r_8> cdiff = (flmnOrig[i] -  flmn[i]) * conj(flmnOrig[i] -  flmn[i]);
    r_16 diff = sqrt(cdiff.real());
    if(diff>err_abs){ 
      err_abs = diff;
      imax = i;
    }
    complex<r_8> foriCAbs = flmnOrig[i]*conj(flmnOrig[i]);
    r_8 foriAbs = sqrt(foriCAbs.real());
    r_8 relatdiff = diff/foriAbs;
	
    if((relatdiff)>err_rel) err_rel = relatdiff;

  }
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  cout << "Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 


}//MultiSphericalLaguerreTransform

//---------------------------------------
void SpinMultiSphericalLaguerreTransform() {

  r_8 alpha = param.alpha;  
  int Nmax = param.N;
  int Lmax = param.Lmax;
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;
  int spin = param.spin; //JEC 18/1/16

  if(spin == 0)  
    throw LagSHTError("WARNING: spin =0 : use test MultiSphericalLaguerreTransform"); 

  cout << " ___________  SpinMultiSphericalLaguerreTransform   TEST _____________ " << endl;

  tstack_push("data input");

  int state = 1234567 + 8912 ; //random seed

  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi,
					 alpha
					 );


  int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
  int Nshell = sphlagtrans.GetOrder();
  int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
  int Npix =  sphlagtrans.GetSphericalGeometry()->Npix();
  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels


  cout << "Verif: Ntheta, Nphi, Npix, Nptot, Nalm, Nshell, Ntot: "
       << sphlagtrans.GetSphericalGeometry()->NTheta() << " "
       << sphlagtrans.GetSphericalGeometry()->NPhi() << " "
       <<  Npix << " "
       <<  NpTot << " "
       <<  Nalm << " "
       <<  Nshell << " "
       <<  Ntot << endl;


  unsigned int maxmemsize =  getMemorySize()/1e6;
  unsigned int estimateMem = 2*8*((uint_8)(NpTot+2*Ntot))/1e6;
  cout << "Estimate usage memory: " << estimateMem << " MBytes" << endl;
  if ( estimateMem > (0.9*maxmemsize) ) {
    cout << ">>>>> Warning: estimate Mem " << estimateMem <<" MB > 90\% "<< maxmemsize << " MB" << endl;
  }//test memory


  vector< complex<r_8> > ElmnOrig(Ntot);
  vector< complex<r_8> > BlmnOrig(Ntot);
  int id
=-1;

  for(int n=0;n<Nmax;n++){ //on each shell
    for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
      for(int l=m;l<Lmax;l++){
	id++;
	r_8 rv = (l<abs(spin)) ? 0.0 : drand(-1,1,&state);
	r_8 iv = (m==0 || l<abs(spin)) ? 0.0 : drand(-1,1,&state);
	ElmnOrig[id] = complex<r_8>(rv,iv);	    
	rv = (l<abs(spin)) ? 0.0 : drand(-1,1,&state);
	iv = (m==0 || l<abs(spin)) ? 0.0 : drand(-1,1,&state);
	BlmnOrig[id] = complex<r_8>(rv,iv);	    
      }//end l-loop
    }//end m-loop
  }//end loop on shell

#if DEBUG >=2 
  std::copy(ElmnOrig.begin(), ElmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
  std::copy(BlmnOrig.begin(), BlmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
#endif
      
  
  tstack_pop("data input");
  tstack_push("processing");
  tstack_push("processing part Synthesis");

#if DEBUG >= 1
  cout << "Main:Synthesis complex function..." << endl;
#endif
  vector<r_8> fijk_re; //NpTot
  vector<r_8> fijk_im;
//   vector< complex<r_8> > Elmk; //calcul intermediaire
//   vector< complex<r_8> > Blmk;

//  sphlagtrans.Synthesis(ElmnOrig, BlmnOrig, spin, fijk_re, fijk_im, Elmk, Blmk);
  sphlagtrans.Synthesis(ElmnOrig, BlmnOrig, spin, fijk_re, fijk_im);
#if DEBUG >= 2
  for (int i=0; i<NpTot; i++){
    cout << "fijk("<<i<<"): " << fijk_re[i] << ", " << fijk_im[i] << endl;
  }
#endif
  
  tstack_pop("processing part Synthesis");
  tstack_push("processing part Analysis");
  
#if DEBUG >= 1
  cout << "Main:Analysis complex function..." << endl;
#endif
  vector< complex<r_8> > Elmn(Ntot);
  vector< complex<r_8> > Blmn(Ntot);

//   vector< complex<r_8> > Elmk_ana;
//   vector< complex<r_8> > Blmk_ana;

//  sphlagtrans.Analysis(fijk_re, fijk_im, spin, Elmn, Blmn, Elmk_ana, Blmk_ana);
  sphlagtrans.Analysis(fijk_re, fijk_im, spin, Elmn, Blmn);
    
  tstack_pop("processing part Analysis");
  tstack_pop("processing");  


  cout << "Error analysis..." << endl;

  r_8 err_abs(0.);
  r_8 err_rel(0.);
  int imax = -1;
  for(int i=0;i<Ntot;i++){
    if(i>(Ntot-9)) {
      cout << "(" << i << ") : Elnm Orig " <<  ElmnOrig[i] << " <-> Elmn Rec " <<  Elmn[i] << endl;
      cout << "............  Blnm Orig " <<  BlmnOrig[i] << " <-> Blmn Rec " <<  Blmn[i] << endl;
    }

    complex<r_8> cdiff = ((ElmnOrig[i] -  Elmn[i]) * conj(ElmnOrig[i] -  Elmn[i])
			  + (BlmnOrig[i] -  Blmn[i]) * conj(BlmnOrig[i] -  Blmn[i]))/2. ;
    r_16 diff = sqrt(cdiff.real());
    if(diff>err_abs){ 
      err_abs = diff;
      imax = i;
    }
    complex<r_8> foriCAbs = (ElmnOrig[i]*conj(ElmnOrig[i]) + BlmnOrig[i]*conj(BlmnOrig[i]))/2. ;
    r_8 foriAbs = sqrt(foriCAbs.real());
    r_8 relatdiff = diff/foriAbs;
	
    if((relatdiff)>err_rel) err_rel = relatdiff;

  }
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  cout << "Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 


}//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;
  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");
  recomputeJlnp = true;
  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
			   clenshawDir,jlnpDir,recomputeJlnp);
  tstack_pop("Ctor lag2bess");

  cout << " ______________ TestJlnpComputation  End  ___________ " << endl;

}



//---------------------------------------------------------------------------------
// void TestJlnpComputationByFFTAdaptative() {
//   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;
//   double alpha = param.alpha;

//   string jlnpDir = param.jlnpDir;
//   string clenshawDir = param.clenshawDir;
 
//   cout << " ______________ TestJlnpComputation FFT Adaptative START ___________ " << endl;
//   LaguerreSphericalTransform sphlagtrans(geometry,
// 					 Lmax,
// 					 -1,
// 					 Nmax,
// 					 Rmax,
// 					 ntheta,
// 					 nphi
// 					 );


//   tstack_push("Ctors....");  
//   LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
//   Bessel jnu(Lmax,std::max(Pmax,Nmax));  //the max is here just in case Nmax =/= Pmax


//   r_8 tau = lagTrans->GetTau();
//   r_8 rNm1 = Rmax / tau; // R = r[N-1] with N the original transform value. Todo (lagTrans_->GetNodes()).back()

//   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;


//   tstack_pop("Ctors....");  
  

//   tstack_push("Bounds init...");  

//   vector<r_8> integUpperBound(Nmax);
//   integUpperBound[0] = min(100.,rNm1); //this is the case of L_0^{\alpha}(x) =1 with no root
//   integUpperBound[1] = min(120.,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
//   }
//   tstack_pop("Bounds init...");  

//   cout << "(JEC) ComputeJInteg  integUpperBound END " << endl;


//   tstack_push("Common FFTW Plan init...");  

//   cout << "(JEC) Plan Init START" << endl;

//   //Try a common order for K & J function this make the FFTW plan for once
//   int iOrdGen = 8; //initial 
//   int nOrdGen=pow(2.,(r_8)iOrdGen)-1; //the -1 is just to use FFTW with 2^iOrdGen 
//   vector<r_8> VecDCT1(nOrdGen+1,0.);  //but for  FFTW_REDFT00 it is requiered n+1... see FFTW doc
//   FFTPlanning planJK(nOrdGen,VecDCT1);

//   cout << "(JEC) Plan Init END" << endl;

//   cout << "(JEC) Plan Inverse Init START" << endl;

//   //Unique Product plan
//   int nOrdProd= 2*nOrdGen+1; //order max of product is 2*nOrdGen but to get a power of 2 for the FFT I add 1.
//   vector<r_8> VecDCT1Inv(nOrdProd+1,0.);
//   FFTPlanning planInv(nOrdProd,VecDCT1Inv);

//   cout << "(JEC) Plan Inverse Init END" << endl;

//   tstack_pop("Common FFTW Plan init...");  
  
//   tstack_push("Common CC weights init...");  

//   cout << "(JEC) CC weights START" << endl;

//   //Clenshow-Curtis single quadratuer weight
//   vector<r_8> wCC(nOrdProd+1,0);
//   FFTPlanning planCC(nOrdProd, wCC);

//   ClenshawCurtisWeightsFast(nOrdProd+1,planCC,wCC);
// #if DEBVEC>0
//   VectorDebug("DCT CCurtis W [0,10]: ",wCC,0,10);
//   VectorDebug("DCT CCurtis W [last-10,last]: ",wCC,nOrdProd-11,nOrdProd+1);
// #endif



//   cout << "(JEC) CC weights END" << endl;

//   tstack_pop("Common CC weights init...");  

//   stringstream ss; ss << "Adap-OMP-FFTClass-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 p=0;p<10;p++){
    
    
//     for(int l=0; l<Lmax; l++){
//       //for(int l=0; l<10; l++){

//       r_8 qlp = jnu(l,p);

//       r_8 klptau = qlp/rNm1; 
//       JBess jFunc(l,klptau);
  
//       r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau;
//       r_8 ql2scaled = jnu(l,1)/klptau;   //  the 2nd BesselRoot rescaled
  


//       for (int n=0; n<Nmax; n++) {
// 	//for (int n=0; n<10; n++) {
  
// 	r_8 intUppBound = std::max(integUpperBound[n], ql2scaled);
	
// #if DEBUG>10

// 	cout << "lowBnd, uppBnd: " << setprecision(30) 
// 	     << intLowBound << ", " << intUppBound 
// 	     << setprecision(6)
// 	     << endl;
// 	cout << "l,n,p,klptau: " 
// 	     << l << ", " 
// 	     << n << ", " 
// 	     << p << ", " 
// 	     << klptau << endl;
// #endif
  
// 	tstack_push("Integral Computation...");
  
// 	tstack_push("Reconfigure FFT planning");
	
// 	//estimate if one needs to reconfigure the FFT plan
// 	int nJ = EstimChebyshevOrdSpherBessel(l,klptau,intLowBound,intUppBound);
// 	int nK =0;
// 	{
// 	  LaguerreFuncQuad lag(n,alpha);
// 	  vector<r_8> nodes;
// 	  vector<r_8> weights;
// 	  lag.QuadWeightNodes(nodes,weights);
// 	  nK = EstimChebyshevOrdLagFunc(nodes,intLowBound,intUppBound);
// 	}
// 	int newiOrdGen=max(nJ,nK);

// 	if (newiOrdGen>iOrdGen) {
// 	  //destroy old and create new plan
// 	  iOrdGen = newiOrdGen;

// 	  cout << "(JEC) new order : " <<  newiOrdGen << " lnp = " 
// 	       << l << " "
// 	       << n << " "
// 	       << p << " "
// 	       << endl;

// 	  planJK.DestroyPlan();
// 	  nOrdGen=pow(2.,(r_8)iOrdGen)-1;
// 	  VecDCT1.resize(nOrdGen+1); //here NO RE-INITIALIZATION to specific value 
// 	  planJK.CreatePlan(nOrdGen,VecDCT1);
	  
// 	  planInv.DestroyPlan();
// 	  nOrdProd= 2*nOrdGen+1;
// 	  VecDCT1Inv.resize(nOrdProd+1);
// 	  planInv.CreatePlan(nOrdProd,VecDCT1Inv);

// 	  planCC.DestroyPlan();
// 	  wCC.resize(nOrdProd+1);
// 	  planCC.CreatePlan(nOrdProd, wCC);

// 	  //recompute Clenshaw-Curtis weights
// 	  ClenshawCurtisWeightsFast(nOrdProd+1,planCC,wCC);
// 	}
	
// 	tstack_pop("Reconfigure FFT planning");


// 	tstack_push("Chebyshev J...");  


// 	//Get Chebyshev coeff of Bessel
// 	//  int iOrdJ = 8; int nOrdJ=pow(2.,(r_8)iOrdJ);
// 	int nOrdJ = nOrdGen; //USE COMMON dimension
// 	vector<r_8> coefJ(nOrdJ+1,0.);
// 	{
// 	  tstack_push("Chebyshev J SAMPLING...");  
// 	  fill(VecDCT1.begin(), VecDCT1.end(), 0.0); //may be unusefull here
// 	  ChebyshevSampling(nOrdJ,VecDCT1,jFunc,intLowBound,intUppBound);

// #if DEBVEC>0
// 	  VectorDebug("Sampling J-func",VecDCT1,0,10);
// #endif
// 	  tstack_pop("Chebyshev J SAMPLING...");  
// 	  tstack_push("Chebyshev J FFT...");  

// 	  ChebyshevCoeffFFT(nOrdJ,planJK, VecDCT1);
// 	  copy(VecDCT1.begin(),VecDCT1.end(),coefJ.begin());

// 	  tstack_pop("Chebyshev J FFT...");  

// #if DEBVEC>0
// 	  VectorDebug("Che J-func",coefJ,0,10);
// #endif
// 	}

// 	tstack_pop("Chebyshev J...");  


// 	tstack_push("Chebyshev K...");  
  
// 	//Get Chebyshev coeff of Laguerre Func
// 	//  int iOrdK = 6; int nOrdK = pow(2.,(r_8)iOrdK);
// 	int nOrdK =  nOrdGen; //USE COMMON dimension
// 	vector<r_8> coefK(nOrdK+1,0.);
// 	{
// 	  tstack_push("Chebyshev K SAMPLING..."); 
// 	  LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
// 	  KLag kLagFunc(Ln, facts[n]);
// 	  fill(VecDCT1.begin(), VecDCT1.end(), 0.0); //may be unusefull here
// 	  ChebyshevSampling(nOrdK,VecDCT1,kLagFunc,intLowBound,intUppBound);

// #if DEBVEC>0
// 	  VectorDebug("Sampling K-func",VecDCT1,0,10);
// #endif
// 	  tstack_pop("Chebyshev K SAMPLING..."); 
// 	  tstack_push("Chebyshev K FFT..."); 
// 	  ChebyshevCoeffFFT(nOrdK, planJK, VecDCT1);
// 	  copy(VecDCT1.begin(),VecDCT1.end(),coefK.begin());

// 	  tstack_pop("Chebyshev K FFT..."); 
// #if DEBVEC>0
// 	  VectorDebug("Che K-func",coefK,0,10);
// #endif

// 	  delete Ln;
// 	}
// 	tstack_pop("Chebyshev K...");  
  

// 	tstack_push("Chebyshev Product...");  

// 	//Product size
// 	//   int nOrdProd = nOrdJ*nOrdK;
// 	//perform Inverse Chebyshev transform in the new basis

// 	tstack_push("Inverse Chebyshev J...");  

// 	vector<r_8>coefJExt(nOrdProd+1,0.); 
// 	fill(VecDCT1Inv.begin(), VecDCT1Inv.end(), 0.0);
// 	copy(coefJ.begin(),coefJ.end(),VecDCT1Inv.begin());  //make the VecDCT1Inv[0:nOrdJ-1] = coefJ; then VecDCT1Inv is comleted by 0
// 	InverseChebyshevCoeffFFT(nOrdProd,planInv,VecDCT1Inv);
// 	copy(VecDCT1Inv.begin(), VecDCT1Inv.end(), coefJExt.begin());

// #if DEBVEC>0
// 	VectorDebug("Inv Che J",coefJExt,0,10);
// #endif

// 	tstack_pop("Inverse Chebyshev J...");  
// 	tstack_push("Inverse Chebyshev K...");  

// 	vector<r_8>coefKExt(nOrdProd+1,0.); 
// 	fill(VecDCT1Inv.begin(), VecDCT1Inv.end(), 0.0);
// 	copy(coefK.begin(),coefK.end(),VecDCT1Inv.begin());
// 	InverseChebyshevCoeffFFT(nOrdProd,planInv,VecDCT1Inv);
// 	copy(VecDCT1Inv.begin(), VecDCT1Inv.end(), coefKExt.begin());

// #if DEBVEC>0
// 	VectorDebug("Inv Che K",coefKExt,0,10);
// #endif

// 	tstack_pop("Inverse Chebyshev K...");  


// 	tstack_push("Inverse Chebyshev J*K...");  

// 	vector<r_8>invCoefProd(nOrdProd+1,0.); //is there an in-place product?
// 	transform(coefJExt.begin(),coefJExt.end(),
// 		  coefKExt.begin(),
// 		  invCoefProd.begin(),
// 		  multiplies<r_8>());

// #if DEBVEC>0
// 	VectorDebug("Inv Che Prod",invCoefProd,0,10);
// #endif
  
// 	tstack_pop("Inverse Chebyshev J*K...");  
// 	tstack_pop("Chebyshev Product...");  
  
// 	//Compute integral
// 	tstack_push("CC quadrature...");
// 	r_8 integral = inner_product(invCoefProd.begin(),invCoefProd.end(),wCC.begin(),0.);
// 	integral *= (intUppBound - intLowBound)/2.0; //the 2 division comme from CC quadrature weights definition in [-1,1];
// 	tstack_pop("CC quadrature...");

// #if CHECK_BY_CCQUAD>0
// 	cout << "Integral via FFT               = " << setprecision(20) << integral <<  setprecision(6) << endl;
// #endif
// 	ofs <<l<<" "
// 	    <<n<<" "
//  	    <<p<< " " << setprecision(30) << integral 
// 	    << endl;

// 	tstack_pop("Integral Computation...");


// #if CHECK_BY_CCQUAD>0
// 	tstack_push("Recursive CC quadrature");
 
// 	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;
// 	r_8 tol4Integral = 1e-10;
  
// 	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
// 	IntFunc1D f(l,klptau,Ln);
// 	theQuad.SetFuncBounded(f,intLowBound,intUppBound);
// 	integ_val = theQuad.GlobalAdapStrat(tol4Integral);
// 	delete Ln;
// 	tstack_pop("Recursive CC quadrature");

// 	cout << "Integral via Recursive CC Quad = " 
// 	     << setprecision(20) 
// 	     << integ_val.first * sqrt(2.0/M_PI) * facts[n]
// 	     << setprecision(6) 
// 	     << endl;
// #endif

//       }//n-loop
//     }//l-loop
//   }//p-loop

//   //file to save computations
//   ofs.close();

//   cout << " ______________ TestJlnpComputation FFT Adaptative START ___________ " << endl;


// }//TestJlnpComputationByFFTAdaptative


//-------------------------------------------------------------
void TestLaguerre2Bessel() {

  //Todo introduce Pmax !

  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 << " ______________ TestLaguerre2Bessel TEST ___________ " << endl;
  
   cout << "tstack 1 Start" <<endl;
   tstack_push("data input");

  int state = 1234567 + 8912 ; //random seed

  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );

  
  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
  LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();

  int Nalm = sphere->Nalm();
  int Nshell = sphlagtrans.GetOrder();
  int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
  int Npix =  sphere->Npix();
  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels

  cout << "Verif: Npix, Nptot, Nalm, Nshell, Ntot: "
       <<  Npix << " "
       <<  NpTot << " "
       <<  Nalm << " "
       <<  Nshell << " "
       <<  Ntot << endl;

  vector< complex<r_8> > flmnOrig(Ntot);
  int id=-1;
  for(int n=0;n<Nmax;n++){ //on each shell
    for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
      for(int l=m;l<Lmax;l++){
	id++;
	
	//verif
	int idLag =  n*Nalm + l+m*Lmax-m*(m+1)/2 ;
	if(id != idLag)
	  throw LagSHTError("ERROR: FLag index error"); 

	r_8 rv = drand(-1,1,&state);
	r_8 iv = (m==0) ? 0.0 : drand(-1,1,&state);
	flmnOrig[id] = complex<r_8>(rv,iv);	    
      }//end l-loop
    }//end m-loop
  }//end loop on shell

#if DEBUG >=2 
  cout << "F-Lag coeff orig ......... START " << endl;
  //  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, "\n"));
  
  for(int l=0; l<Lmax; l++){
    for (int m=0; m<=l; m++) {
      int almoff7 = l+m*Lmax-m*(m+1)/2;
      for (int n=0; n<Nmax; n++){
	int idLag =  n*Nalm + almoff7;
	cout << l << " " 
	     << m << " "
	     << n << " "
	     << setprecision(12) << flmnOrig[idLag] << endl;
      }
    }
  }
  
  cout << "F-Lag coeff orig ......... END " << endl;
#endif
        
  cout << "tstack 1 End" <<endl;
  tstack_pop("data input");


  cout << "tstack 2 Start" <<endl;
  tstack_push("processing");
  cout << "tstack 2a Start" <<endl;
  tstack_push("Init Laguerre 2 Bessel");
  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
			   clenshawDir,jlnpDir,recomputeJlnp);
  cout << "tstack 2a End" <<endl;
  tstack_pop("Init Laguerre 2 Bessel");

  cout << "tstack 2b Start" <<endl;
  tstack_push("Compute Fourier-Bessel coeffs.");

  vector< complex<r_8> > FBlmp(Nalm*Pmax);
  lag2bess.Lag2BessCoeff(flmnOrig,FBlmp);
#if DEBUG >=2 
  cout << "F-Bessel coeff ......... START " << endl;
  std::copy(FBlmp.begin(), FBlmp.end(), std::ostream_iterator< complex<r_8> >(std::cout, "\n"));
  cout << "F-Bessel coeff ......... END " << endl;
#endif
  cout << "tstack 2b End" <<endl;
  tstack_pop("Compute Fourier-Bessel coeffs.");
  

  cout << "tstack 2c Start" <<endl;
  tstack_push("Compute Alm(rk) from FB coeffs.");
  vector< complex<r_8> > FBalmk;
  vector<r_8> fFBijk;
  lag2bess.Synthesis(FBlmp,FBalmk,fFBijk);
  cout << "tstack 2c End" <<endl;
  tstack_pop("Compute Alm(rk) from FB coeffs.");
  

  cout << "tstack 2d Start" <<endl;
  tstack_push("Compute Alm(rk) from FL coeffs.");
  vector< complex<r_8> > FLalmk;
  vector<r_8> fFLijk;
  sphlagtrans.Synthesis(flmnOrig,fFLijk,FLalmk);
#if DEBUG >=2 
  cout << "F-Lag fijk coeff ......... START " << endl;
  std::copy(fFLijk.begin(), fFLijk.end(), std::ostream_iterator< complex<r_8> >(std::cout, "\n"));
  cout << "F-Lag fijk coeff ......... END " << endl;
#endif

  cout << "tstack 2d End" <<endl;
  tstack_pop("Compute Alm(rk) from FL coeffs.");
  
  cout << "tstack 2 End" <<endl;

  tstack_pop("processing");


  //verif F-Laguerre part
  tstack_push("F-L Analysis verif");
  vector< complex<r_8> > flmn(Ntot);
  sphlagtrans.Analysis(fFLijk,flmn);
  {
    r_8 err_abs(0.);
    r_8 err_rel(0.);
    int imax = -1;
    for(int i=0;i<Ntot;i++){
      
      complex<r_8> cdiff = (flmnOrig[i] -  flmn[i]) * conj(flmnOrig[i] -  flmn[i]);
      r_16 diff = sqrt(cdiff.real());
      if(diff>err_abs){ 
	err_abs = diff;
	imax = i;
      }
      complex<r_8> foriCAbs = flmnOrig[i]*conj(flmnOrig[i]);
      r_8 foriAbs = sqrt(foriCAbs.real());
      r_8 relatdiff = diff/foriAbs;
	
      if((relatdiff)>err_rel) err_rel = relatdiff;
    }
    cout << " >>>>>>>>>>>>>>>>>>>>> Fourrier-Laguerre part <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
    cout << "Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
    cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
    
  }
  tstack_pop("F-L Analysis verif");


  
  
  {//Check 1b
    cout << "Dump FL or FB reconstructed Cl(r_k)" <<endl;

#if DEBUG >=2 
    std::ofstream ofs;
    ofs.open ("Cli.txt", std::ofstream::out);
#endif
    
    for(int n=0;n<Nmax;n++){
      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<<" " << n << " " << setprecision(12) << ClFL << " " << ClFB << endl;
#endif	  

      }//l
    }//n
#if DEBUG >=2 
    ofs.close();
#endif
  }//check 1b






  {//check 2
#if DEBUG >=2 
    std::ofstream ofs;
    ofs.open ("fijk.txt", std::ofstream::out);
#endif
    cout << "Dump FL or FB reconstructed fijk on each shell k" <<endl;

    for (int ish=0;ish<Nshell; ish++){
      r_8 err_abs(0.);
      r_8 err_rel(0.);
      for (int ip=0; ip<Npix; ip++) {
	int id = ish*Npix+ip;
	//	cout << "Shell("<<ish<<") Pix["<<ip<<"] FL fijk = " << fFLijk[id] << " FB fijk = " << fFBijk[id] << endl;

#if DEBUG >=2 
	ofs << ish << " " << ip << " " << setprecision(12) <<  fFLijk[id] << " " << fFBijk[id] << endl;
#endif
	r_8 diff = fabs(fFLijk[id]-fFBijk[id]);
	if(diff>err_abs){
	  err_abs = diff;
	}
	r_8 relatdiff = diff/ fFLijk[id];
	if(relatdiff>err_rel) err_rel = relatdiff;
      }//loop on pîx
      
            
      cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<ish<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
      cout << "Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
      cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 

    }//loop on shell


#if DEBUG >=2 
    ofs.close();
#endif
  }//check 2

  

}//TestLaguerre2Bessel
//-------------------------------------------------------------

void TestPixelization() {

  int Nmax = param.N;
  int Lmax = param.Lmax;
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;


  
  LaguerreSphericalTransform lagsht(geometry,
				    Lmax,
				    -1,
				    Nmax,
				    Rmax,
				    ntheta,
				    nphi
				    );
  
  ntheta = lagsht.GetSphericalGeometry()->NTheta(); //identique a Nrings
  nphi   = lagsht.GetSphericalGeometry()->NPhi();

  int state = 1234567 + 8912 ; //random seed
  
  stringstream ss; ss << "pixels-" << geometry << "-L" << Lmax << ".txt";
  std::ofstream ofs (ss.str().c_str(), std::ofstream::out);

  r_8 dTheta = M_PI/ntheta/2;
  r_8 dPhi   = 2*M_PI/nphi/2;
  
  cout << "half theta Width : " << dTheta << ", half phi Width : " << dPhi << endl; 

  
  int nloop = 100000;

  for (int i=0;i<nloop; i++) {
    r_8 theta_test = drand(0.,M_PI,&state);
    r_8 phi_test = drand(0.,2*M_PI,&state);
    int idx = lagsht.GetSphericalGeometry()->PixIndexSph(theta_test, phi_test);

    r_8 theta_out;
    r_8 phi_out;
    lagsht.GetSphericalGeometry()->PixThetaPhi(idx,theta_out, phi_out);
    
    r_8 dphi = phi_test-phi_out;
    if(dphi>M_PI)dphi=2*M_PI-dphi;

    ofs << theta_test << " " << phi_test << " " << theta_out << " " << phi_out << " " << dphi<< endl;
  }
  
}//TestPixelization

//----------------------------------------------
//               Main 
//----------------------------------------------

int main(int narg, char *arg[]) {

   unsigned int maxmemsize =  getMemorySize()/1e6;
   cout << "Max Memory size: " <<  maxmemsize << " MBytes" <<  endl;

  int N = 128; 
  r_8 R = 1.;
  int Lmax = 128;
  int Pmax = 0;
  int spin = 0;

  int test=1;

  string geometry = "Gauss";
  int ntheta = -1;
  int nphi = -1;

  char* pLAGSHTROOTDIR; pLAGSHTROOTDIR = getenv("LAGSHTROOTDIR");
  if(pLAGSHTROOTDIR == 0) 
    throw LagSHTError("ERROR: LAGSHTROOTDIR no defined => execute the setup"); 
    
  string clenshawDir = string(pLAGSHTROOTDIR)+"/data";
  string jlnpDir = string(pLAGSHTROOTDIR)+"/data";
  bool recomputeJlnp = false;


  //JEC 12/1/16 alpha as input
  double alpha = 2.0; //default

  int ka=1;
  while (ka<narg) {
    if (strcmp(arg[ka],"-h")==0) {
      cout << "usage: ./Objs/lagsht_testsuite -t <test number> [1]\n"
	   << " [-a <alpha> [2 by default] ] "
	   << " [-n <Nmax> [128]] [-l <Lmax> [128]]\n" 
	   << " [-g <geometry> Gauss|Fejer1|Healpix [Gauss]] \n"
	   << " [-ntheta <number of theta rings> [determined by the geometry]\n     in case of Healpix gives the nside parameter]\n"
	   << " [-nphi <number of pixel per rings> [determined by the geometry]]\n"
	   << " [-p <Pmax for Fourier-Bessel> [Nmax] ]\n"
	   << " [-r <Rmax> [1.]]\n"
	   << " [-spin <spin value> [0]"
	   << endl; 
      return 0;
    }
    else if (strcmp(arg[ka],"-a")==0) {
      alpha = atof(arg[ka+1]);
      ka+=2;      
    }    
    else if (strcmp(arg[ka],"-spin")==0) {
      spin = atoi(arg[ka+1]);
      ka+=2;      
    }    
    else if (strcmp(arg[ka],"-n")==0) {
      N=atoi(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-l")==0) {
      Lmax=atoi(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-t")==0) {
      test=atoi(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-g")==0) {
      geometry=arg[ka+1];
      ka+=2;      
    }
    else if (strcmp(arg[ka],"-ntheta")==0) {
      ntheta = atoi(arg[ka+1]);
      ka+=2;            
    }
    else if (strcmp(arg[ka],"-nphi")==0) {
      nphi = atoi(arg[ka+1]);
      ka+=2;      
    }
    else if (strcmp(arg[ka],"-r")==0) {
      R = atof(arg[ka+1]);
      ka+=2;      
    }    
    else if (strcmp(arg[ka],"-p")==0) {
      Pmax = atof(arg[ka+1]);
      ka+=2;      
    }    
    else if (strcmp(arg[ka],"-clenPath")==0) {
      clenshawDir = arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-jlnpPath")==0) {
      jlnpDir = arg[ka+1];
      ka+=2;
    }
    else if(strcmp(arg[ka],"-jlnp")==0) {
      recomputeJlnp = true;
      ka+=1;
    }
    else ka++;
  }//eo while


  param.Lmax = Lmax;
  param.N    = N;
  param.Pmax = Pmax;
  param.R    = R;
  param.spin  = spin;
  param.alpha = alpha;
  param.geometry = geometry;
  param.ntheta = ntheta;
  param.nphi = nphi;

  param.clenshawDir = clenshawDir;
  param.jlnpDir     = jlnpDir;
  param.recomputeJlnp = recomputeJlnp;

  cout << "Configuration parameters are set to: "
       << " Test number : " << test << "\n"
       << " Lmax, Nmax, Pmax, spin, Rmax , alpha: " 
       << param.Lmax << ", "
       << param.N << ", "
       << param.Pmax << ", "
       << param.spin << ", "
       << param.R << ", "
       << param.alpha << "\n"
       << "Geometry: " << param.geometry << " Ntheta, Nphi: " 
       << param.ntheta << ", " << param.nphi<<"\n";

  int rc=0;
  try {
    
    switch(test) {
    case 0:
      TestBasic();
      break;
    case 1:
      LaguerreQuadrature();
      break;
    case 2:
      MultiLaguerreTransformSdtAlone();
      break;
    case 3:
      TestPixelization();
      break;
    case 4:
      {
	if(N>5500) throw LagSHTError("WARNING: N>5500 the algorithm may fail due to weight estimates"); 
	
      
// 	unsigned int estimateMem = 8*N*(uint_8)((2*Lmax+1)*(2*Lmax+1))/1e6; //Npix * 8Bytes and then in MBytes
// 	if ( estimateMem > (0.9*maxmemsize) ) {
// 	  char buff[80];
// 	  cout << "estimate Mem usage " << estimateMem << " MB" << endl;
// 	  sprintf(buff,"Main FATAL: estimate Mem %u MB > (90%%) Max Mem %u MB ",estimateMem,maxmemsize);
// 	  throw  LagSHTError(buff);
// 	}//test memory
	
	tstack_push("MultiSphericalLaguerreTransform");
	MultiSphericalLaguerreTransform();
	tstack_pop("MultiSphericalLaguerreTransform");
	tstack_report("MultiSphericalLaguerreTransform");  
      }
      break;

    case 50:
      {
	cout << "rebuild Jlnp: \n Clenshaw Path<"<<param.clenshawDir<<">\n and jlnp Path<"
	     << param.jlnpDir << ">"
	     << endl;
	if(Pmax<N){ param.Pmax = N; cout << "(Warning): Pmax<Nmax => modify the value to "<<  param.Pmax << endl;}
	tstack_push("TestJlnpComputation");
	TestJlnpComputation();
 	tstack_pop("TestJlnpComputation");
 	tstack_report("TestJlnpComputation");
      }
      break;

    case 51:
      {
	if(recomputeJlnp){
	  cout << "rebuild Jlnp: \n Clenshaw Path<"<<param.clenshawDir<<">\n and jlnp Path<"
	       << param.jlnpDir << ">"
	       << endl;
	} else {
	  cout << "Jlnp path <" << param.jlnpDir << ">" << endl;
	}
	if(Pmax<N){ param.Pmax = 16*N; cout << "(Warning): Pmax<Nmax => modify the value to "<<  param.Pmax << endl;}
	tstack_push("TestLaguerre2Bessel");
	TestLaguerre2Bessel();
 	tstack_pop("TestLaguerre2Bessel");
 	tstack_report("TestLaguerre2Bessel");
      }
      break;


    case 6:	
      {

	unsigned int estimateMem = 8*2*N*(uint_8)((2*Lmax+1)*(2*Lmax+1))/1e6; //Npix * 8Bytes and then in MBytes
	if ( estimateMem > (0.9*maxmemsize) ) {
	  char buff[80];
	  cout << "estimate Mem usage " << estimateMem << " MB" << endl;
	  sprintf(buff,"Main FATAL: estimate Mem %u MB > (90%%) Max Mem %u MB ",estimateMem,maxmemsize);
	  throw  LagSHTError(buff);
	}//test memory


	tstack_push("SpinMultiSphericalLaguerreTransform");
	SpinMultiSphericalLaguerreTransform();
	tstack_pop("SpinMultiSphericalLaguerreTransform");
	tstack_report("SpinMultiSphericalLaguerreTransform");  
      }
      break;
      
    default:
      throw LagSHTError("EROOR Test type unkwown");
    }//end of switch


    cout << "---/ Fin bloc try ---- " << endl;}
    
  catch (LagSHTError & e) {
    cerr << " lagsh_testsuite.cc: Catched Exception (LagSHTError)" << (string)typeid(e).name() 
	 << " - Msg= " << e.what() << endl;
    rc = 99;
  }
  catch (std::exception & e) {
    cerr << " lagsh_testsuite.cc: Catched std::xception "  
	 << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {
    cerr << " lagsh_testsuite.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  cout << " ---- Programme lagsh_testsuite.cc -  FIN  (Rc=" << rc << ") --- " << endl;
  return rc;
}//main