#include <string.h> //strcmp
#include <utility>      // std::pair

//LagSHT
#include "LagSHT/lagsht_numbers.h"
#include "LagSHT/lagsht_exceptions.h"
#include "LagSHT/lagsht_utils.h" //getMemorySize
#include "LagSHT/walltimer.h"    
#include "LagSHT/lagSphericTransform.h"



using namespace LagSHT;


template <typename T>
class Vector2d {
public:
  Vector2d(size_t d1=0, size_t d2=0, T const & t=T()) :
    d1(d1), d2(d2), data(d1*d2, t){ totsize = d1*d2; }
  
  T & operator()(size_t i, size_t j) {
    return data[i + j*d1];
  }
  
  T const & operator()(size_t i, size_t j, size_t k) const {
    return data[i + j*d1];
  }
  
  size_t size() const { return totsize; }

private:
  size_t d1,d2;
  size_t totsize;
  vector<T> data;
};


template <typename T>
class Vector3d {
public:
  Vector3d(size_t d1=0, size_t d2=0, size_t d3=0, T const & t=T()) :
    d1(d1), d2(d2), d3(d3), data(d1*d2*d3, t){ totsize = d1*d2*d3; }
  
  T & operator()(size_t i, size_t j, size_t k) {
    return data[i + j*d1 + k*d1*d2];
  }
  
  T const & operator()(size_t i, size_t j, size_t k) const {
    return data[i + j*d1 + k*d1*d2];
  }
  
  size_t size() const { return totsize; }

private:
  size_t d1,d2,d3;
  size_t totsize;
  vector<T> data;
};


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


//---------- Utility functions and class for this test program

namespace LagSHT {


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

  static double gaussrand(int* state,double mean=0., double sigma=1.0) {
    r_8 rd1 = drand(0.,1.,state);
    while (rd1 == 0.) rd1=drand(0.,1.,state);
    r_8 rd2 = drand(0.,1.,state);
    return (sqrt(-2.*log(rd1))*cos(2.*M_PI*rd2))*sigma + mean;
  }

  class TestTransform : public GaussGeometry {
  public:
    //! Constructor
    TestTransform(int Lmax, int Mmax, int Nrings, int Nphi): 
      GaussGeometry(Lmax, Mmax, Nrings, Nphi) {
    }//Ctor
    
    //! Destructor
    virtual ~TestTransform() {}

    //! Forward map -> coeff
    void Forward(const vector<r_8>& map, vector<complex<r_8> >& coeff){
//       cout << "Forward map: " << map.size() << endl;
//       cout << "coeff: " << coeff.size() << endl;
      map2alm(&(map.data()[0]),reinterpret_cast<r_8*>(&coeff.data()[0]),false);
    }

    //! Backward coeff -> map
    void Backward(const vector< complex<r_8> >& coeff, vector<r_8>& map) {
//       cout << "Backward map: " << map.size() << endl;
//       cout << "coeff: " << coeff.size() << endl;
      alm2map(reinterpret_cast<const r_8*>(&coeff.data()[0]),&(map.data()[0]),false);
    }
    
    //
    void ComputeCl(const vector<complex<r_8> >& alm, vector<r_8>& cl){
      if(cl.size() != (size_t)L_) {
	cout << "WARNING: resising Cl" << endl;
	cl.resize(L_,0.);
      }
      for(int l=0;l<L_;l++) {
       for (int m=0;m<=l;m++){
 	 //sharp LambdaLM
	 size_t i = l+m*L_-m*(m+1)/2;
	 //	 cout << "(l,m): " << l << "," << m << " a-> " << alm[i] << endl;
	 r_8 alm2 = alm[i].real()*alm[i].real() +  alm[i].imag()*alm[i].imag();	 
	 cl[l] += (m==0) ? alm2 : 2.0*alm2;
      }
       cl[l] /= (r_8)(2*l+1);
      }
    }//ComputeCl

    void GenAlmFromCl(const vector<r_8>& cl, vector< complex<r_8> >& alm){
      int state=1234567 + 8912 ; //random seed
     
      for(int l=0;l<L_;l++) {
	r_8 rms = sqrt(cl[l]);
       for (int m=0;m<=l;m++){
 	 //sharp LambdaLM
	 size_t i = l+m*L_-m*(m+1)/2;
	 if(m==0){
	   alm[i] = complex<r_8>(gaussrand(&state,0.,rms));
	 } else {
	   complex<r_8>tmp = complex<r_8>(gaussrand(&state,0.,1.),
					  gaussrand(&state,0.,1.));
	   alm[i] = ((r_8)(rms*M_SQRT1_2))*tmp;
	 }
       }
      }      
    }//GenAlmFromCl


  };


}//end namespace



//------------------------------------
//test3: 
// generate Clin(k) k here index of shells "a la" CMB Cl gaussian profile
// compute Alm(k) from Clin(k)
// compute a map(k) a from Inverse Spherical Harmonic Transform and fill fijk 3D-vector
// proceed to LagSHT Analysis transform to get flmn and AlmOut(k) coeff.
// compute Clout(k) from AlmOut(k) which should be identical to Clin(k)
// compute Clout(n) a la CMB by reducing the "m" index
//------------------------------------
void test3() {

  //initial values
  int Lmax = param.Lmax;
  int Nmax = param.N;
  r_8 Rmax = param.R;
  string geom = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;


  tstack_push("Generate fijk 3D-ball");
  cout << "Generate Cl(k) spectra and fill a 3D-ball" << endl;
  //Class to perform 2D Ylm transform with the convention in lagSHT
  TestTransform blst(Lmax,-1,-1,-1);

  int npix = blst.Npix(); 
  int nalms = blst.Nalm();


  vector<r_8> fijk(npix*Nmax,0); // the 3D pixelization

  Vector2d<r_8> Clkin(Lmax,Nmax,0); // Generate Cl(k) with k index of subshell
  for(int k=0;k<Nmax;k++){ //shell index

    //    cout << "Shell[" << k << "]" << endl;
    
    // Generate true Cl  Profile as in 2D case but evolving in k-direction
    vector<r_8> Clin(Lmax); // Cl on the subshell
    for(int l=0;l<Lmax;l++){
      r_8 mean = Lmax/3.*(1. + k/(Nmax-1.)); // a k dependance Lmax/3 -> 2Lmax/3
      r_8 sigma = Lmax/20.; //fix
      r_8 arg = (l-mean)/sigma; 
      arg *= arg;
      Clkin(l,k) = exp(-arg/2.);
      Clkin(l,k) /= (r_8)(2*l+1); //nomalisation
      Clin[l] = Clkin(l,k);
    }//l-loop

    //From the Clin compute the Alm and generate a 2D map
    vector< complex<r_8> > alm(nalms);
    vector<r_8> shmap(npix,0);
    blst.GenAlmFromCl(Clin,alm);
    blst.Backward(alm,shmap);
    
    //transfert the 2D Map into the 3D pixelization
    //    cout << "Transfert Start at fijk[" << k*npix << "]" << endl;
    copy(shmap.begin(),shmap.end(),fijk.begin()+k*npix);

  }//shell-loop

  tstack_pop("Generate fijk 3D-ball");

  tstack_push("Lag Spherical Transf.");
  cout << "Perform the Laguerre Spherical Transf. to get both the flmn coef. and the alm on each shell" << endl;
  //Perform 3D fijk -> flmn & flmk
  LaguerreSphericalTransform 
    sphlagtrans(geom, //type of geometry
		Lmax, //l=0,...,Lmaw-1
		-1,   //Mmax default = Lmax set if <0
		Nmax, //Nmax : n:0...Nmax-1
		Rmax, //Rmax : R in [0, Rmas]
		ntheta, //number of collatitude rings (or nsides ) default computed if <0
		nphi //number of phi per rings default computed if <0
		); 

  int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
  int Nshell = sphlagtrans.GetOrder();
  int Ntot = Nshell * Nalm; //total number of 3D-pixels & number of Coefficients of the Spherical Laguerre transform
  //  int Npix = sphlagtrans.GetSphericalGeometry()->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> > flmn(Ntot);
  vector< complex<r_8> > almk(Ntot);
  sphlagtrans.Analysis(fijk,flmn, almk);    // <---------- the job

  tstack_pop("Lag Spherical Transf.");

  tstack_push("Clk and Cln OUT");
  cout << "Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn" << endl;

  //Compute the Clkout with k the index of subshell and should be comapred to Clkin from the almk
  //And the Cl(n) n the 3rd index of the flmn
  //Compute back the Cl(l,n)
  //and the Cl(l,k) from the alm(k)
  Vector2d<r_8> Clnout(Lmax,Nmax,0);
  Vector2d<r_8> Clkout(Lmax,Nmax,0);
  for(int n=0;n<Nmax;n++){
    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 
	r_8 flmn2 = flmn[id].real()*flmn[id].real() + flmn[id].imag()*flmn[id].imag();
	r_8 almk2 = almk[id].real()*almk[id].real() + almk[id].imag()*almk[id].imag();
	Clnout(l,n) += (m==0) ? flmn2 : 2.0*flmn2;
	Clkout(l,n) += (m==0) ? almk2 : 2.0*almk2;
      }
      Clnout(l,n) /= (r_8)(2*l+1); //nomalisation "a la Cl 2D" 
      Clkout(l,n)  /= (r_8)(2*l+1);
    }
  }

  tstack_pop("Clk and Cln OUT");

  tstack_push("error analysis");

  cout << "Error analysis on Cl(k) input and output of the Laguerre Spherical Transf." << endl;
  r_8 err_abs(0.);
  r_8 err_rel(0.);
  pair<int,int> imax = make_pair(-1,-1);
  for(int k=0;k<Nmax;k++){ //shell index
    for(int l=0;l<Lmax;l++){
      r_8 diff = Clkin(l,k)-Clkout(l,k); diff *=diff;
      if (diff > err_abs) {
	err_abs = diff;
	imax = make_pair(k,l);
      }
      r_8 relatdiff = diff/fabs((r_8)Clkin(l,k));
      if((relatdiff)>err_rel) err_rel = relatdiff;
    }
  }
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  cout << " Err. Max. " << err_abs << " [" << imax.first << ", " 
       << imax.second << "], Err. Rel. " << err_rel << endl;
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  tstack_pop("error analysis");


//   {//save

// #ifdef SOPHYALIB

//     char buff[80];
//     sprintf(buff,"test3-L%d-N%d.ppf",Lmax,Nmax);
//     POutPersist po(buff);
//     po << PPFNameTag("clkin") << Clkin;
//     po << PPFNameTag("clkout") << Clkout;
//     po << PPFNameTag("clnout") << Clnout;
    
//   }//save  
}//TEST 3




//-------------------------------------
// test5: 
// o generated a 3D Cube with Nc^3 pixels filled with random gaussian number N(0,sigma) : fixed sigma
// o fill 3D Ball embedded into the Cube 
// o perform a LagSHT Analysis to get the flmn and almk coefficients
// o compute the Cl(n) and the Cl(k) by reducing the "m" index. Notice that the Cl(k) are CMB like Cl on each
//   subshells
void test5() {

  using namespace LagSHT;
  
  //input values
  int Lmax = param.Lmax;
  int Nmax = param.N;
  r_8 Rmax = param.R;        // <----------------- will be changed below
  string geom = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;



  //the 3D Cube of Nc^3 pixels 
  int state=1234567 + 8912;

  int Nc = 256*4;

  
  tstack_push("Generate a 3D-Cube");
  cout << "Generate a 3D-Cube " << Nc << "^3" << endl; 
  
  Vector3d<r_8> cube(Nc,Nc,Nc);
  r_8 rms = 1.0; // constant RMS for the moment
  for(size_t ix=0; ix<(size_t)Nc; ix++) {
    for(size_t iy=0; iy<(size_t)Nc; iy++) {
      for(size_t iz=0; iz<(size_t)Nc; iz++) {
	cube(ix,iy,iz) = (r_8)gaussrand(&state,0.,rms);
      }
    }
  }

  tstack_pop("Generate a 3D-Cube");

  

  tstack_push("Fill the  3D-Ball from the Cube");
  cout << "Fill the 3D-Ball from the Cube" << endl;

  r_8 Lo = 256 * 8.0 ; //8Mpc per cell
  r_8 Lx = Lo; //Total dimension of the Cube  
  r_8 Ly = Lo; 
  r_8 Lz = Lo; 

  //Gives physical dimensions of the Cube
  r_8 Dx = Lx/Nc; //Mpc a cell X dimension
  r_8 Dy = Ly/Nc; //Mpc a cell Y dimension
  r_8 Dz = Lz/Nc; //Mpc a cell Z dimension
  

  
  //define Rmax to include a 3D Ball in the 3D Cube
  Rmax = min(min(Lx,Ly),Lz)/2.*0.8; //80% of the minimal half length of the cube
  


  cout << " Fill a 3D-Ball (" << Dx << "," << Dy << "," << Dz << ")Mpc^3: "
       << "Rmax is set to : " << Rmax << endl;


  LaguerreSphericalTransform 
    sphlagtrans(geom, //type of geometry
		Lmax, //l=0,...,Lmaw-1
		-1,   //Mmax default = Lmax set if <0
		Nmax, //Nmax : n:0...Nmax-1
		Rmax, //Rmax : R in [0, Rmas]
		ntheta, //number of collatitude rings (or nsides ) default computed if <0
		nphi //number of phi per rings default computed if <0
		); 

  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
  if(sphere == NULL) throw LagSHTError("test5 sphere == NULL");

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

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


  LaguerreTransform* lagtrans = sphlagtrans.GetLagTransform(); //the pure Laguerre Part
  if(lagtrans == NULL) throw LagSHTError("test5 lagtrans == NULL");

  vector<r_8> fijk(NpTot,0); // the 3D-ball

  //scan the pixels of the different shells of the 3D Ball
  for(int ks=0; ks<Nshell; ks++){ //shell k
    r_8 rShell = lagtrans->Index2R(ks);
    //    cout << "Shell["<<ks<<"] radius = " << rShell << endl;
    if(rShell > Rmax) throw LagSHTError("test5 a shell radius too large...");
    for(int ipix=0; ipix<Npix; ipix++){ //pixels of the shell
      r_8 theta; //here by convention theta = 0 means North Pole
      r_8 phi; 
      sphere->PixThetaPhi(ipix,theta,phi);
      //3D xyz position of the pixel
      r_8 xpix = rShell * sin(theta) * cos(phi); 
      r_8 ypix = rShell * sin(theta) * sin(phi); 
      r_8 zpix = rShell * cos(theta);
	    
      //find the corresponding Cube cell ijk indices
      int i = int((xpix + Lx/2.0 - Dx/2.0)/Dx);
      if(i<0||i>Nc)throw LagSHTError("test5 i out of range");
      int j = int((ypix + Ly/2.0 - Dy/2.0)/Dy);
      if(j<0||j>Nc)throw LagSHTError("test5 j out of range");
      int k = int((zpix + Lz/2.0 - Dz/2.0)/Dz);
      if(k<0||k>Nc)throw LagSHTError("test5 k out of range");
      
      //affect the 3D-cube ijk cell value to the shell ipix
      int iter = ks*Npix+ipix; //location of the ipix of k-th shell
      if(iter<0 || iter>NpTot)throw LagSHTError("test5 iter out of range");
      
      fijk[iter] = cube(i,j,k);
    }//end loop pixels of a shell
  }//end loop shell

  tstack_pop("Fill the  3D-Ball from the Cube");

  tstack_push("Lag Spherical Transf.");
  cout << "Perform the Laguerre Spherical Transf. to get both the flmn coef. and the alm on each shell" << endl;

  
  //Get the flmn and almk coefficients
  vector< complex<r_8> > flmn(Ntot);
  vector< complex<r_8> > almk(Ntot);
  sphlagtrans.Analysis(fijk,flmn, almk);    // <---------- the job


  tstack_pop("Lag Spherical Transf.");

  tstack_push("Clk and Cln OUT");
  cout << "Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn" << endl;

  //Compute the Clkout with k the index of subshell and should be comapred to Clkin from the almk
  //And the Cl(n) n the 3rd index of the flmn
  //Compute back the Cl(l,n)
  //and the Cl(l,k) from the alm(k)
  Vector2d<r_8> Clnout(Lmax,Nmax,0);
  Vector2d<r_8> Clkout(Lmax,Nmax,0);
  for(int n=0;n<Nmax;n++){
    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 
	r_8 flmn2 = flmn[id].real()*flmn[id].real() + flmn[id].imag()*flmn[id].imag();
	r_8 almk2 = almk[id].real()*almk[id].real() + almk[id].imag()*almk[id].imag();
	Clnout(l,n) += (m==0) ? flmn2 : 2.0*flmn2;
	Clkout(l,n) += (m==0) ? almk2 : 2.0*almk2;
      }
      Clnout(l,n) /= (r_8)(2*l+1); //nomalisation "a la Cl 2D" 
      Clkout(l,n)  /= (r_8)(2*l+1);
    }
  }

  tstack_pop("Clk and Cln OUT");

//   //Compute the Clkout with k the index of subshell and should be comapred to Clkin from the almk
//   //And the Cl(n) n the 3rd index of the flmn
//   //Compute back the Cl(l,n)
//   //and the Cl(l,k) from the alm(k)
//   TMatrix<r_8> Clnout(Lmax,Nmax);
//   TMatrix<r_8> Clkout(Lmax, Nmax);
//   for(int n=0;n<Nmax;n++){
//     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 
// 	r_8 flmn2 = flmn[id].real()*flmn[id].real() + flmn[id].imag()*flmn[id].imag();
// 	r_8 almk2 = almk[id].real()*almk[id].real() + almk[id].imag()*almk[id].imag();
// 	Clnout(l,n) += (m==0) ? flmn2 : 2.0*flmn2;
// 	Clkout(l,n)   += (m==0) ? almk2 : 2.0*almk2;
//       }
//       Clnout(l,n) /= (r_8)(2*l+1); //nomalisation "a la Cl 2D" 
//       Clkout(l,n)  /= (r_8)(2*l+1);
//     }
//   }

//   {
//     char buff[80];
//     sprintf(buff,"test5-Nc%d-L%d-N%d.ppf",Nc,Lmax,Nmax);
//     POutPersist po(buff);
//     po << PPFNameTag("clkout") << Clkout;
//     po << PPFNameTag("clnout") << Clnout;
//   }  
}//test 5

//----------------------------------------------
//               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 = 5.;
  int Lmax = 1024;
  int Pmax = 0;

  int test=3;

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

  int ka=1;
  while (ka<narg) {
    if (strcmp(arg[ka],"-h")==0) {
      cout << "usage: ./bin/test3D -t <test number> [1]\n"
	   << " [-n <Nmax>] [-l <Lmax>]\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.]]"
	   << endl; 
      return 0;
    }
    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 ka++;
  }//eo while

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

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


  int rc=0;
  try {
    
    switch(test) {
    case 3:
      tstack_push("Test 3");
      test3();
      tstack_pop("Test 3");
      tstack_report("Test 3");
      break;
    case 5:
      tstack_push("Test 5");
      test5();
      tstack_pop("Test 5");
      tstack_report("Test 5");
      break;
    default:
      throw LagSHTError("EROOR Test type unkwown");
    }//end of switch


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