#include <iostream.h>
#include "sambainit.h"
#include "spheregorski.h"

#include "tod.h"
#include "timing.h"


int main(int narg, char* arg[]) 
{
  double teta,phi;
  PeidaInit();
  InitTim();   // Initializing the CPU timer 
  if ((narg > 1) && (strcmp(arg[1],"-h") == 0) )  {
    cout << " tspm [Gorski_M=32] : Gorski Spherical Map Test " << endl;
    exit(0);
    }

  int  m=32;
  if (narg >1)  m = atoi(arg[1]);
  cout <<  " ===== Gorski Spherical Map Test M= " << m << endl;

  SphereGorski<double> sph(m);

  cout << "Filling spherical map NPixels= " << sph.NbPixels() << endl;
  for (int j=0;j<sph.NbPixels();j++) 
    {
      sph.PixThetaPhi(j,teta,phi);
      sph(j)= 0.2* cos(3.*teta)*sin(8*phi);
    }
  PrtTim("End of Fill ");

  // Computing mean and sigma on the sphere 
  double gmoy=0. , gsig = 0.; 
  double valok;
  for(int k=0; k<sph.NbPixels(); k++) {
    valok = sph(k);
    gmoy += valok;  gsig += valok*valok; 
  }
  cout << "SphMap Mean= " << gmoy << "  Sigma = " << gsig << endl;
  PrtTim("End of Mean-Sig ");

//  Writing to a PPF file
  {
  POutPersist s("sphg.ppf"); 
  FIO_SphereGorski<double> fiog(sph) ;
  fiog.Write(s);
  cout << "SphMap written to sphg.ppf "  << endl;
  }

  // Reading from the file
  {
  FIO_SphereGorski<double> fiog("sphg.ppf");
  double gmoy=0. , gsig = 0.; 
  double valok;
  SphereGorski<double>  sph2 = fiog;
  PrtTim("End of Write/Read ");

  int ndiff = 0;
  for(int k=0; k<sph2.NbPixels(); k++) {
    valok = sph2(k);
    gmoy += valok;  gsig += valok*valok; 
    if ((sph2(k)-sph(k)) > 1.e-49)  ndiff++;
  }
  cout << "SphMapFromFile Mean= " << gmoy << "  Sigma = " << gsig << endl;
  cout << " NDiff = " << ndiff << " (should be zero = 0) " << endl;
  PrtTim("End of Mean-Sig ");
  }

  cout <<  " ===== Fin de TSPM_Test ======== " << endl;
  return 0;
}