#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
#define DEBVEC 0
#define NUM_OMP_THREADS   4
#define CHECK_BY_CCQUAD 0

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

}
//---------------------------------------------------------------------------------
class JBess : public ClassFunc1D {
public:
  JBess(int l, r_8 kt):l_(l),kt_(kt) {norm_ = sqrt(2./M_PI);}
  virtual double operator()(double x) const {
    return norm_*Bessel::jn(l_,x*kt_);
  }
  virtual ~JBess() {}
private:
  int l_;
  r_8 kt_;
  r_8 norm_;
};
//------------------
class KLag : public ClassFunc1D {
public:
  KLag(LaguerreFuncQuad* Ln, r_8 norm):Ln_(Ln),norm_(norm){}
  virtual double operator()(double x) const {
    return norm_*x*x*Ln_->Value(x);
  }
  virtual ~KLag() {}
private:
  LaguerreFuncQuad* Ln_; //no ownership
  r_8 norm_;
};
//------------------
class IntFunc1D : public ClassFunc1D {
public:
  IntFunc1D(int l, r_8 klp,LaguerreFuncQuad* Ln, r_8 norm=1.0):
    l_(l),klp_(klp),Ln_(Ln),norm_(norm) {}
  virtual double operator()(double x) const {
    return norm_*(x*x)*(Bessel::jn(l_,x*klp_))*Ln_->Value(x);
  }
  virtual ~IntFunc1D() {}
private:
  int l_;
  r_8 klp_;
  LaguerreFuncQuad* Ln_; //no ownership
  r_8 norm_;
};//class IntFunc1D
//------------------
void ChebyshevSampling(int n, vector<r_8>& val, const ClassFunc1D& f, r_8 a, r_8 b){
  r_8 bma = 0.5*(b-a); r_8 bpa = 0.5*(b+a);
  r_8 cte = M_PI/((r_8)n);

  int k;
  r_8 x;
#pragma omp parallel for private(x) num_threads(NUM_OMP_THREADS)
  for(k=0;k<=n;k++){
    x = cos(k*cte)*bma+bpa;
    val[k] = f(x);
  }
}
//------------------
class FFTPlanning {
public:
  FFTPlanning(int n, vector<r_8>& vec) {
    createPlan(n,vec);
  }
  ~FFTPlanning() { DestroyPlan(); }
  void DestroyPlan() {
    if(plan_) {
      cout << "Destroy Plan ["<< this << "]" << endl;
      fftw_destroy_plan(plan_); 
      plan_=NULL;
    }
  }
  void Execute() const { fftw_execute(plan_); }
private:
  void createPlan(int n, vector<r_8>& vec) {
    cout << "Create Plan ["<< this << "]" << endl;
    plan_ = fftw_plan_r2r_1d(n + 1,&vec.data()[0], &vec.data()[0], FFTW_REDFT00, FFTW_ESTIMATE);
  }
  fftw_plan plan_;
};
//------------------
//size of val = n+1 
//void ChebyshevCoeffFFT(int n, const fftw_plan plan, vector<r_8>& val){
void ChebyshevCoeffFFT(int n, const FFTPlanning& plan, vector<r_8>& val){


  tstack_push("FFTW plan exec...");
  //  fftw_execute ( plan );
  plan.Execute();
  tstack_pop("FFTW plan exec...");
  /*
    Chebyshev Coeff. the noramization of FFTW is propto  1/(val.size()-1) = 1/n
  */
  tstack_push("FFTW norm...");
  r_8 norm = 1./((r_8)n); 
  transform(val.begin(), val.end(), val.begin(),
		std::bind1st(multiplies<r_8>(),norm));

  val[0] *= 0.5;
  val[n] *= 0.5;
  tstack_pop("FFTW norm...");

}
//------------------
void InverseChebyshevCoeffFFT(int n, const FFTPlanning& plan, vector<r_8>& val){
  val[0] *= 2.0;
  val[n] *= 2.0;
  //  fftw_execute ( plan );
  plan.Execute();
  
  transform(val.begin(), val.end(), val.begin(),
		std::bind1st(multiplies<r_8>(),0.5));

}
//------------------
void VectorDebug(const string& msg, const vector<r_8>&vec, int beg, int end){
  cout << "(JEC) debug " << msg << endl;
  for(int i=beg;i<end;i++) {
    cout << vec[i] << ", ";
  }
  cout << endl;
}
//------------------
//extrait de quadinteg.h dans un premier temps.
//rappel l'original calcule les poids (et nodes) pour une quadrature sur [-1, 1]
void ClenshawCurtisWeights(int  order, vector<r_8>& w) throw(string) {
  double b;
  int i;
  int j;
  double pi = M_PI; //JEC use cmath definition of PI
  double theta;

  //Todo this kind of error should be tracked before
  if ( order < 1 ) {
    cerr << "\n";
    cerr << "CLENSHAW_CURTIS_COMPUTE - Fatal error!\n";
    cerr << "  Illegal value of ORDER = " << order << "\n";
    exit (EXIT_FAILURE);
  }
  
 
  if ( order == 1 ) {
    w[0] = 2.0;
  } else {
    
    for ( i = 0; i < order; i++ ) {
      theta = ( double ) ( i ) * pi / ( double ) ( order - 1 );  
      w[i] = 1.0;
      
      for ( j = 1; j <= ( order - 1 ) / 2; j++ ) {
	if ( 2 * j == ( order - 1 ) ) {
	  b = 1.0;
	} else {
	  b = 2.0;
	}
	
	w[i] = w[i] - b *  cos ( 2.0 * ( double ) ( j ) * theta )
	  / ( double ) ( 4 * j * j - 1 );
      }
    }
    
    w[0] = w[0] / ( double ) ( order - 1 );
    for ( i = 1; i < order - 1; i++ ) {
      w[i] = 2.0 * w[i] / ( double ) ( order - 1 );
    }
    w[order-1] = w[order-1] / ( double ) ( order - 1 );
  }
}

//weights defined for [-1 1]
//FFTW Y_k = 2 Sum''_j=0^(N-1) X_j Cos[Pi k j/(N-1)]
//CC weights
//W_k = 4/(N-1) a_k  Sum''_{j=0, j even}^(N-1) 1/(1-j^2) Cos[Pi k j/(N-1)]
void ClenshawCurtisWeightsFast(int n, const FFTPlanning&  plan, vector<r_8>& w) {

  //dim of w is n
  fill(w.begin(), w.end(), (r_8)0.0);

  for(int k=0;k<n; k +=2){
    w[k] = 1./(1.-(r_8)(k*k));
  }
  
  //  fftw_execute(plan);
  plan.Execute();
  

  r_8 norm = 2.0/(r_8)(n-1); //2 * FFTW DCT-1 
  transform(w.begin(), w.end(), w.begin(),
		std::bind1st(multiplies<r_8>(),norm));
  w[0] /= (r_8)2; w[n-1] /= (r_8)2;
  
}

/*
  estimates the order of the Chebyshev expansion approx using the number
  of roots and min/max of the cosine function asymptotic apprx of Spherical Bessel
  in the range [lowBnd, uppBnd]
 */
int EstimChebyshevOrdSpherBessel(int el, r_8 kt, r_8 lowBnd, r_8 uppBnd){
  int n=0;
  int qq = - ( el/2 +1 );
  r_8 val = M_PI/(2.0* kt) * (2*qq + el + 2 );
  while(val<lowBnd) {
    qq++;
    val = M_PI/(2.0* kt) * (2*qq + el + 2);
  }
  while(val<=uppBnd){
    qq++;
    n++;
    val = M_PI/(2.0* kt) * (2*qq + el + 2);
  }

  int nNodesAndMinMax = 2*n; //factor 2 to take into account min/max
  
  return floor(log2((r_8)nNodesAndMinMax))+1;
}
/*
  estimates the order of the Chebyshev expansion approx using the number
  of roots of Laguerre polynomial in the range [lowBnd, uppBnd]
 */
int EstimChebyshevOrdLagFunc(const vector<r_8>& nodes, r_8 lowBnd, r_8 uppBnd) {

  //If C++11 not available find an alternative with Class fucntion
  int n= count_if(nodes.begin(), nodes.end(), [&](r_8 x){ return (x>=lowBnd)&&(x<=uppBnd);} );

  if(n==0)n=1; //no root for Laguerre[0,x]

  int nNodesAndMinMax = 2*n; //factor 2 to take into account min/max
  return floor(log2((r_8)nNodesAndMinMax))+1;

}


//---------------------------------------------------------------------------------
void TestJnlpOrder() {
  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;
  cout << "_______________ TestJnlpOrder Start _____________ " << endl;
  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );

  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()

  vector<r_8> integUpperBound(Nmax);
  integUpperBound[0] = min(100.,rNm1); //this is the case of L_0^{\alpha}(x) =1 with no root (L
  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
  }
  
  stringstream ss; ss << "FFT-order"<<Lmax<<"-N"<<Nmax<<"-P"<<Pmax<<".txt";
  std::ofstream ofs;
  string fname("./"+ss.str());
  ofs.open (fname.c_str(), std::ofstream::out);


  for(int p=0;p<Pmax;p++){
    for(int l=0; l<Lmax; 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++) {
  
	r_8 intUppBound = std::max(integUpperBound[n], ql2scaled);
	
	int nJ = EstimChebyshevOrdSpherBessel(l,klptau,intLowBound,intUppBound);

	LaguerreFuncQuad lag(n,alpha);
	vector<r_8> nodes;
	vector<r_8> weights;
	lag.QuadWeightNodes(nodes,weights);

	int nK = EstimChebyshevOrdLagFunc(nodes,intLowBound,intUppBound);
	
	ofs << l << " " 
	     << n << " "
	     << p << " "
	     << nJ << " " 
	     << nK << endl;

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

  ofs.close();

  cout << "_______________ TestJnlpOrder End _____________ " << endl;

}//TestJnlpOrder
//---------------------------------------------------------------------------------
void TestJlnpComputationByFFT() {
  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 START ___________ " << endl;
  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );


  tstack_push("Ctors....");  
  //  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
  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 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()

  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 = 10; 
  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 << "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 klptau =1.; //JUST FOR TEST wrt M10
  

      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("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 END ___________ " << endl;


}


//---------------------------------------------------------------------------------
// 
//Version 0 du calcul des Jlnp par FFTW mais non optimisee pour les plans.
// 
//
// void ChebyshevCoeffFFT(int n, vector<r_8>& val){
//   /* Planning In-place
//    */
//   tstack_push("FFTW plan init...");
//   fftw_plan plan;
//   plan = fftw_plan_r2r_1d(n + 1,&val.data()[0], &val.data()[0], FFTW_REDFT00, FFTW_ESTIMATE);
//   tstack_pop("FFTW plan init...");

//   tstack_push("FFTW plan exec...");
//   fftw_execute ( plan );
//   tstack_pop("FFTW plan exec...");
//   /*
//     Chebyshev Coeff.
//   */
//   tstack_push("FFTW norm...");
//   r_8 norm = 1./((r_8)n);
//   transform(val.begin(), val.end(), val.begin(),
// 		std::bind1st(multiplies<r_8>(),norm));
//   //  for(int k=0;k<=n;k++) val[k] *= norm;
//   val[0] *= 0.5;
//   val[n] *= 0.5;
//   tstack_pop("FFTW norm...");

//   fftw_destroy_plan(plan);
// }
// //------------------
// void InverseChebyshevCoeffFFT(int n, vector<r_8>& val){
//   val[0] *= 2.0;
//   val[n] *= 2.0;
//   //  for(int k=0;k<=n;k++) val[k] *= n; //not necessary
//   fftw_plan plan;
//   plan = fftw_plan_r2r_1d(n + 1,&val.data()[0], &val.data()[0], FFTW_REDFT00, FFTW_ESTIMATE);
//   fftw_execute ( plan );
//   //  for(int k=0;k<=n;k++) val[k] /= (2.0*n); //not necessary
//   for(int k=0;k<=n;k++) val[k] *= 0.5;
//   fftw_destroy_plan(plan);
// }
// //------------------
// void VectorDebug(const string& msg, const vector<r_8>&vec, int beg, int end){
//   cout << "(JEC) debug " << msg << endl;
//   for(int i=beg;i<end;i++) {
//     cout << vec[i] << ", ";
//   }
//   cout << endl;
// }
// void TestJlnpComputationByFFT() {
//   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 clenshawDir = param.clenshawDir;
 
//   cout << " ______________ TestJlnpComputation FFT START ___________ " << endl;
//   LaguerreSphericalTransform sphlagtrans(geometry,
// 					 Lmax,
// 					 -1,
// 					 Nmax,
// 					 Rmax,
// 					 ntheta,
// 					 nphi
// 					 );


//   tstack_push("Ctors....");  
//   //  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
//   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 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()

//   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] = 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
//   }
//   tstack_pop("Bounds init...");  

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

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

//   int p=256;
//   int l=64;
//   int n=32;
//   r_8 qlp = jnu(l,p);

//   //----------->   r_8 klptau = qlp/rNm1; 
  
//   r_8 klptau =1.; //JUST FOR TEST wrt M10
  
//   cout << "n,l,p,klptau: " 
//        << n << ", " 
//        << l << ", " 
//        << p << ", " 
//        << klptau << endl;

//   r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau;
//   r_8 ql2scaled = jnu(l,1)/klptau;   //  the 2nd BesselRoot rescaled
  
//   r_8 intUppBound = std::max(integUpperBound[n], ql2scaled);

//   cout << "lowBnd, uppBnd: " << setprecision(30) 
//        << intLowBound << ", " << intUppBound 
//        << setprecision(6)
//        << endl;

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


//   //Get Chebyshev coeff of Bessel
//   int iOrdJ = 8; int nOrdJ=pow(2.,(r_8)iOrdJ);
//   vector<r_8> coefJ(nOrdJ+1,0.);
//   {
//     JBess jFunc(l,klptau);
//     tstack_push("Chebyshev J SAMPLING...");  
//     ChebyshevSampling(nOrdJ,coefJ,jFunc,intLowBound,intUppBound);

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

//     ChebyshevCoeffFFT(nOrdJ,coefJ);

//     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);
//   vector<r_8> coefK(nOrdK+1,0.);
//   {  
//     LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
//     KLag kLagFunc(Ln, facts[n]);
//     ChebyshevSampling(nOrdK,coefK,kLagFunc,intLowBound,intUppBound);

// #if DEBVEC>0
//     VectorDebug("Sampling K-func",coefK,0,10);
// #endif

//     ChebyshevCoeffFFT(nOrdK,coefK);

// #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;
//   //extend the Chebyshev coeff vector with 0. Is there a simpler way?
//   vector<r_8>coefJExt(nOrdProd+1,0.); copy(coefJ.begin(),coefJ.end(),coefJExt.begin());
//   vector<r_8>coefKExt(nOrdProd+1,0.); copy(coefK.begin(),coefK.end(),coefKExt.begin());
//   //perform Inverse Chebyshev transform in the new basis
//   InverseChebyshevCoeffFFT(nOrdProd,coefJExt);

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

//   InverseChebyshevCoeffFFT(nOrdProd,coefKExt);

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

//   vector<r_8>invCoefProd(nOrdProd); //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("Chebyshev Product...");  
  
//   tstack_push("CCurtis weights...");
//   int iOrdCC = iOrdJ+iOrdK; int nOrdCC = pow(2.0,(r_8)iOrdCC)+1;
//   vector<r_8> wCC(nOrdCC,0);
//   ClenshawCurtisWeights(nOrdCC,wCC);

// #if DEBVEC>0
//   VectorDebug("CCurtis W [0,10]: ",wCC,0,10);
//   VectorDebug("CCurtis W [last-10,last]: ",wCC,nOrdCC-11,nOrdCC);
// #endif
//   tstack_pop("CCurtis weights...");
  
//   //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...");
//   cout << "Integral via FFT               = " << setprecision(20) << integral <<  setprecision(6) << endl;

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


//   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;
  
  
//   cout << " ______________ TestJlnpComputation FFT END ___________ " << endl;


// }

//-------------------------------------------------------------
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 52:
      {
	if(Pmax<N){ param.Pmax = N; cout << "(Warning): Pmax<Nmax => modify the value to "<<  param.Pmax << endl;}
	tstack_push("TestJlnpComputation   FFT");
	TestJlnpComputationByFFT();
	tstack_pop("TestJlnpComputation   FFT");
	tstack_report("TestJlnpComputation   FFT");
      }
      break;


    case 53:
      {
	TestJnlpOrder();
      }
      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