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

#include "machdefs.h"
#include <iostream.h>
#include <math.h>

#include "sambainit.h"
#include "spherehealpix.h"
#include "fiospherehealpix.h"
#include "fitsioserver.h"
#include "nbrandom.h"
#include "bruit.h"
#include "timing.h"

// Test program to generate different type of random skies
template <class T>
void  MeanSig(PixelMap<T> const & map, double& gmoy, double& gsig);

// ------------- Main program --------------
int main(int narg, char* arg[]) 
{
  double teta,phi;
  double gmoy=0., gsig=0.;

  SophyaInit();
  InitTim();   // Initializing the CPU timer 
  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;   
    exit(0);
    }

  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;
  m = 0.;
  s = 1.;
  a = b = 1.;
  K = 1.;
  if ((typ < 0) || (typ > 2)) 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 { 
    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;
    }

  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 {
    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, 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 ");
   }


   
  FitsIoServer fios;
  fios.save(sph, arg[4]);
  cout << "SphereHEALPix<float> written to FITS file " << (string)(arg[4]) << endl;
  PrtTim("End of WriteFITS ");

  cout <<  " ============ End of tgsky program ======== " << endl;
  return 0;
}

/* Nouvelle-Methode */
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);
}