source: Sophya/trunk/SophyaProg/PMixer/tgsky.cc@ 1313

#include "pmixer.h"
#include "mapoperation.h"

/*! \ingroup PMixer
 * \file tgsky.cc
 * \brief \b PROGRAM \b tgsky  <BR>
 * \l Program which generates different types of random skies
 */
// ------------- Main program --------------

int main(int narg, char* arg[]) 
{

  if ((narg < 5) || ((narg>1) && (strcmp(arg[1],"-h") == 0)) )  {
    cout << " tgsky : Generation of random skies  " << endl;
    cout << " Usage: tgsky NSide TypSky Params FitsFileName [PPFName] " << endl;
    cout << "   - TypSky = 0  Params= m,sig  , Gaussian(m,sig) " << endl;
    cout << "   - TypSky = 1  Params= m,sig  , OOFNoise(sig) + m" << endl;
    cout << "   - TypSky = 2  Params= K,a,b,m  K*cos(a*teta)*sin(b*phi)+m" << endl;   
    cout << "   - TypSky = 3  Params= theta,phi DIPOLE" << endl;   
    exit(0);
    }

  try {
  double teta,phi;
  double gmoy=0., gsig=0.;
#if !defined(SunOS) || !defined(__GNUG__)
  // ca se plante sur Sun !
  SophyaInit();
#endif
  InitTim();   // Initializing the CPU timer 
  int  ns = atoi(arg[1]);
  SphereHEALPix<float> sph(ns);

  cout << "Filling SphereHEALPix<float> NSide= " << ns 
       << "  NPixels= " << sph.NbPixels() << endl;

  int typ = atoi(arg[2]);
  double m,s,a,b,K;
  double Theta,Phi;
  Theta = 0.;
  Phi   = 0.;
  m = 0.;
  s = 1.;
  a = b = 1.;
  K = 1.;
  if ((typ < 0) || (typ > 3)) typ = 0;
  if (typ < 2)  { 
    sscanf(arg[3],"%lg,%lg",&m,&s);
    cout << "  TypSky= " << typ << " m= " << m << " sig= " << s << endl;
    if (typ == 1) cout << "  ---> OOFNoise() " << endl;
    else cout << "  ---> NoiseGenerator() " << endl;
  }
  else if(typ==2) { 
    sscanf(arg[3],"%lg,%lg,%lg,%lg",&K,&a,&b,&m);
    cout << "  TypSky= " << typ << " K= " << m << " a= " << a << " b=" 
         << b << " m= " << m << endl;
    cout << "    ---> K * cos(a*teta) * sin(b*teta) + m" << endl;
    }
  else if(typ==3) { 
    sscanf(arg[3],"%lg,%lg",&Theta,&Phi);
    cout << "  TypSky= " << typ << " Theta= " << Theta << " phi= " << Phi ;
    cout << "    ---> DIPOLE " << endl;
    }

  PrtTim("End of Init ");

  if (typ == 2) {
    for (int j=0;j<sph.NbPixels();j++)   {
      sph.PixThetaPhi(j,teta,phi);
      sph(j)= K * cos(a*teta)*sin(b*phi) + m;
    }
  }
  else if(typ==3)
    {
      UnitVector vd(Theta,Phi);
      UnitVector vc(Theta,Phi);
      for(int i=0; i<sph.NbPixels(); i++)
        {
          double Thetar,Phir;
          sph.PixThetaPhi(i,Thetar,Phir);      
          vc.SetThetaPhi(Thetar,  Phir);
          sph(i) += vd.Psc(vc);
        }
    }
  else{
    NoiseGenerator * ng;
    if (typ == 0) ng = new NoiseGenerator(s);
    else ng = new OOFNoise(s);    
    for (int j=0;j<sph.NbPixels();j++)   sph(j) = ng->Noise()+m;
    delete ng;
  }
  PrtTim("End of Fill ");

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

  if (narg > 5) {
   POutPersist s(arg[5]); 
   FIO_SphereHEALPix<float> fiog(&sph) ;
   fiog.Write(s);
   cout << "SphereHEALPix<float> written to POutPersist file"  
        << (string)(arg[5]) << endl;
   PrtTim("End of WritePPF ");
   }


  { 
  FitsOutFile fios(arg[4]);
  fios.firstImageOnPrimaryHeader(false); // Use secondary header 
  DVList dvl;
  dvl["PDMTYPE"] = "COMPMAP";
  dvl.SetComment("PDMTYPE", "Planck Data Model Type"); 
  dvl["SOPHYVER"] =  SophyaVersion();
  dvl.SetComment("SOPHYVER", "Sophya Version"); 
  dvl["TGSKYKW"] =  "Test Keyword from tgsky.cc";
  fios.DVListIntoPrimaryHeader(dvl);  
  sph.Info() = dvl;
  fios << sph; 
  cout << "SphereHEALPix<float> written to FITS file " << (string)(arg[4]) << endl;
  PrtTim("End of WriteFITS ");
  }
  }
  catch (PThrowable & exc) {
    cerr << " tgsky: Catched Exception - Msg= " << exc.Msg() 
         << " Type=" << (string)typeid(exc).name() << endl;
  }
  catch (...) {
    cerr << " tgsky: Catched unknown exception (...) " << endl;
  }
  cout <<  " ============ End of tgsky program ======== " << endl;
  return 0;
}