source: Sophya/trunk/SophyaProg/Tests/tspm.cc@ 3950

#include <iostream>
#include "sopnamsp.h"
#include "skymapinit.h"
#include "skymap.h"

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

template <class T>
void  MeanSig(PixelMap<T> const & map, double& gmoy, double& gsig)
{
  gmoy=0.;
  gsig = 0.;
  double valok;
  for(int k=0; k<map.NbPixels(); k++) {
    valok = map(k);
    gmoy += valok;  gsig += valok*valok;
  }
  gmoy /= (double)map.NbPixels();
  gsig = gsig/(double)map.NbPixels() - gmoy*gmoy;
  if (gsig >= 0.) gsig = sqrt(gsig);
}

int main(int narg, char* arg[]) 
{
  double teta,phi;
  double gmoy, gsig;

  SophyaInit();
  InitTim();   // Initializing the CPU timer 
  if ((narg > 1) && (strcmp(arg[1],"-h") == 0) )  {
    cout << " tspm [HEALPix_M=32] [M_TetaPhi=64] : Spherical Map Test " << endl;
    exit(0);
    }

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

  SphereHEALPix<double> sph(m);

  cout << "Filling spherical map NPixels= " << sph.NbPixels() << endl;
  int j;
  for (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 
  MeanSig(sph, gmoy, gsig);
  cout << "SphereHEALPix<double> Mean= " << gmoy << "  Sigma = " << gsig << endl;
  PrtTim("End of Mean-Sig ");

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

  // Reading from the file
  {
  FIO_SphereHEALPix<double> fiog(string("sphg.ppf"));
  SphereHEALPix<double>  sph2 = fiog;
  PrtTim("End of Write/Read ");

  cout << " Spherical map from file sphg.ppf NPixels= " << sph2.NbPixels() << endl;

  int ndiff = 0;
  for(int k=0; k<sph2.NbPixels(); k++) {
//    if ((sph2(k)-sph(k)) > 1.e-49)  { 
    if ( sph2(k) != sph(k) )   { 
      ndiff++;
      if (ndiff < 20)  cout << "!!!Diff: K= " << k << " SPH2= " << sph2(k) << " SPH= " << sph(k) << endl;
    }
  }
  MeanSig(sph, gmoy, gsig);
  cout << "SphMapFromFile Mean= " << gmoy << "  Sigma = " << gsig << endl;
  cout << " NDiff = " << ndiff << " (should be zero = 0) " << endl;
  PrtTim("End of Mean-Sig ");
  }


  int  mt=64;
  if (narg > 2)  mt = atoi(arg[2]);
  cout <<  "\n ===== ThetaPhi Spherical Map Test MT= " << mt << endl;

  SphereThetaPhi<float> spht(m);

  cout << "Filling spherical map NPixels= " << spht.NbPixels() << endl;
  for (j=0;j<spht.NbPixels();j++) 
    {
      spht.PixThetaPhi(j,teta,phi);
      spht(j)= 0.5* cos(5.*teta)*sin(10.*phi);
    }
  PrtTim("End of Fill ");

  // Computing mean and sigma on the sphere 
  MeanSig(spht, gmoy, gsig);
  cout << "SphereThetaPhi<float> Mean= " << gmoy << "  Sigma = " << gsig << endl;
  PrtTim("End of Mean-Sig ");

//  Writing to a PPF file
  {
  POutPersist s("spht.ppf"); 
  FIO_SphereThetaPhi<float> fiog(spht) ;
  fiog.Write(s);
  cout << "SphereThetaPhi<float> written to spht.ppf "  << endl;
  }

  // Reading from the file
  {
  FIO_SphereThetaPhi<float> fiog("spht.ppf");
  SphereThetaPhi<float>  sph2 = fiog;
  PrtTim("End of Write/Read ");

  cout << " Spherical map from file sph.ppf NPixels= " << sph2.NbPixels() << endl;

  int ndiff = 0;
  for(int k=0; k<sph2.NbPixels(); k++) {
//    if ((sph2(k)-sph(k)) > 1.e-49)  { 
    if ( sph2(k) != spht(k) )   { 
      ndiff++;
      if (ndiff < 20)  cout << "!!!Diff: K= " << k << " SPH2= " << sph2(k) << " SPH= " << sph(k) << endl;
    }
  }
  MeanSig(sph, gmoy, gsig);
  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;
}