#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 "lagSphericTransform.h"
#include "laguerre2bessel.h"

using namespace LagSHT;

#define DEBUG 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;   //JEC 18/1/16
  //  int alpha;
  double alpha;
  string geometry;
  int ntheta;
  int nphi;
  string clenshawDir; //JEC 23/11/15
  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;   //JEC 23/11/15
  string jlnpDir = param.jlnpDir;           //       "
  bool recomputeJlnp = param.recomputeJlnp; //       "

  cout << " ______________ TestJlnpComputation  START ___________ " << endl;
  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );

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


  tstack_push("Ctor lag2bess");
  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
			   clenshawDir,jlnpDir,recomputeJlnp);

  tstack_pop("Ctor lag2bess");

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

}


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

  //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 1
    cout << "Dump FL or FB reconstructed Alm(r_k)" <<endl;

#if DEBUG >=2 
    std::ofstream ofs;
    ofs.open ("Almi.txt", std::ofstream::out);
#endif
    
    for(int n=0;n<Nmax;n++){
      r_8 err_abs(0.);
      r_8 err_rel(0.);

      for(int l=0;l<Lmax;l++){
	for (int m=0;m<=l;m++) {
	  int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 

#if DEBUG >=2 
	  ofs <<l<<" "
	      <<m<<" "
	      <<n<< " " << setprecision(12) << FLalmk[id].real() << " " <<   FLalmk[id].imag()
	      << "  " << FBalmk[id].real() << " " << FBalmk[id].imag()  << endl;
#endif	  
	  complex<r_8> cdiff = (FLalmk[id] -  FBalmk[id])*conj(FLalmk[id] -  FBalmk[id]);
	  r_8 diff = sqrt(cdiff.real());
	  if(diff>err_abs){
	    err_abs = diff;
	  }
	  complex<r_8> cref= FLalmk[id]*conj(FLalmk[id]);
	  r_8 ref = sqrt(cref.real());
	  r_8 relatdiff = diff/ref;
	  if(relatdiff>err_rel) err_rel = relatdiff;
	}
      }
      
      cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<n<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
      cout << "Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
      cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 

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

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

  //JEC 23/11/15
  char* pLAGSHTROOTDIR; pLAGSHTROOTDIR = getenv("LAGSHTROOTDIR");
  string clenshawDir = string(pLAGSHTROOTDIR)+"/data";
  string jlnpDir = string(pLAGSHTROOTDIR)+"/data";
  bool recomputeJlnp = true;


  //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 5:
      {
	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();
	TestJlnpComputation();
 	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