source: Sophya/trunk/SophyaProg/PrgMap/cl2map.cc@ 1441

#include "machdefs.h"    // Definitions specifiques SOPHYA

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

#include "nbmath.h"
#include "timing.h"

#include "array.h"
#include "skymap.h"
#include "samba.h"
#include "sambainit.h"
#include "fitsspherehealpix.h"
#include "fitstarray.h"

/*!
  \ingroup PrgUtil
  \file cl2map.cc
  \brief \b cl2map: generates spherical maps from power spectra (C(l))

  \verbatim

  csh> cl2map -h
  SOPHYA Version  1.1 Revision 0 (V_Fev2001) -- Mar  9 2001 15:45:31 cxx

  cl2map : Spherical harmonics synthesis - Power spectrum C_l -> HEALPix map
  Usage: cl2map [-float/-r4] [-msph pixp] [-beam fwhm]
         [-fitsin] [-fitsout] [-autoinirand] InFileName OutFileName
   -float (-r4): single precision C_l and map (default = double)
   -msph pix: Spherical map pixelization parameter
              nside=2^n for HEALPix - default=128
                lmax = 2^n  default=128
   -beam fwhm : Beam moothing (in l-space) FWHM in arcmin, default=0
   -fitsout: Select the FITS format for the output map (default PPF format)
   -fitsin : Select the FITS format for the input vector (default PPF format)
   -autoinirand : Automtatic random number generator initialization
   InFileName : Input file name containing the C_l vector
   OutFileName : Output file name (HEALPix map)

  \endverbatim
*/


/* Nouvelle-Fonction */
void Usage(bool fgerr)
{
  cout << endl;
  if (fgerr) {
    cout << " cl2map : Argument Error ! cl2map -h for usage " << endl; 
    exit(1);
  }
  else {
    cout << " cl2map : Spherical harmonics synthesis - Power spectrum C_l -> HEALPix map \n" 
         << " Usage: cl2map [-float/-r4] [-msph pixp] [-beam fwhm] \n" 
         << "        [-fitsin] [-fitsout] [-autoinirand] InFileName OutFileName \n" 
         << "   -float (-r4): single precision C_l and map (default = double) \n"
         << "   -msph pix: Spherical map pixelization parameter \n" 
         << "              nside=2^n for HEALPix - default=128 \n"
         << "                lmax = 2^n  default=128 \n"
         << "   -beam fwhm : Beam moothing (in l-space) FWHM in arcmin, default=0 \n"
         << "   -fitsout: Select the FITS format for the output map (default PPF format) \n"
         << "   -fitsin : Select the FITS format for the input vector (default PPF format) \n"
         << "   -autoinirand : Automtatic random number generator initialization \n" 
         << "   InFileName : Input file name containing the C_l vector \n" 
         << "   OutFileName : Output file name (HEALPix map) \n" << endl;
    exit(0);
  }  
}

/* Nouvelle-Classe */
template <class T> 
class _Cl2Map {
public :
static void BuildMap(string & infile, string & outfile, int msph, double beam, 
                     bool fgfitsin, bool fgfitsout)
{
  double deg2rad =  M_PI/180.;
  double minute2rad =  M_PI/(180.*60.);

  TArray<T> cla;
  SphereHEALPix<T> sph;
  if (fgfitsin) {
    cout << "--- Reading Input FITS file " << infile << endl; 
    FitsInFile fii(infile);
    fii >> cla;
  }
  else {
    cout << "--- Reading Input PPF file " << infile << endl; 
    PInPersist ppi(infile);
    ppi >> cla;
  }
  TVector<T> clvec(cla);
  cout << " Input vector : " ; 
  clvec.Show(); 
  T min, max;
  double mean, sigma;
  clvec.MinMax(min, max);
  MeanSigma(clvec, mean, sigma);
  cout << " C_l.Min= " << min << " C_l.Max= " << max 
       << " C_l.Mean= " << mean << " C_l.Sigma= " << sigma << endl;

  cout << "--- Calling  GenerateFromCl() (beam= " << beam 
       << "' msph= " << msph << ") beam=" << beam*minute2rad << " rad" << endl;
  // on cree la carte 
  SphericalTransformServer<T> sphtr;
  sphtr.GenerateFromCl(sph, msph, clvec, beam*minute2rad);

  cout << " Output map : NbPixels= " <<  sph.NbPixels() << " NSide= " 
       << sph.SizeIndex() << " Resolution= "
       << sqrt(sph.PixSolAngle(0))/minute2rad << " Arcminutes " << endl;

  if (fgfitsout) {
    cout << "--- Writing map to Output FITS file " << outfile << endl; 
    FitsOutFile fio(outfile);
    fio << sph;
  }
  else {
    cout << "--- Writing map to Output PPF file " << outfile << endl; 
    POutPersist ppo(outfile);
    ppo << sph;
  }
}

};

/* Main-Program */
int main(int narg, char *arg[])
{
  if ((narg > 1) && (strcmp(arg[1],"-h") == 0) ) Usage(false);
  
  int msph = 128;
  double beam = 0.;
  string infile;
  string outfile;
  bool fgfitsin = false;
  bool fgfitsout = false;
  bool fginirand = false;
  bool fgr4 = false;
  cout << " cl2map : Decoding command line options ... " << endl;

  int ko=1;
  for (int k=1; k<narg; k++)   {
    if (strcmp(arg[k], "-msph") == 0)  {
        if (k == narg-1) Usage(true);  // -msph est suivi d'un argument 
      msph = atoi(arg[k+1]);  k++;       // k++ pour sauter au suivant
    }
    else if (strcmp(arg[k], "-beam") == 0)  {
        if (k == narg-1) Usage(true);  // -beam est suivi d'un argument 
      beam = atof(arg[k+1]);  k++;       // k++ pour sauter au suivant
    }
    else if (strcmp(arg[k], "-fitsin") == 0) {
      fgfitsin = true;  
    }
    else if (strcmp(arg[k], "-fitsout") == 0) {
      fgfitsout = true;  
    }
    else if (strcmp(arg[k], "-autoinirand") == 0) {
      fginirand = true;  
    }
    else if ((strcmp(arg[k], "-float") == 0) || (strcmp(arg[k], "-r4") == 0) ) {
      fgr4 = true;  
    }
   
    else { ko = k; break; }  // Debut des noms
  }

  if ((narg-ko) < 2)  Usage(true); 
  infile = arg[ko];
  outfile = arg[ko+1];
  
  // Bloc try de gestion des exception 
  try {
    InitTim();
    SophyaInit();
    if (fginirand) {
      cout << " Setting random number generator seed (using time) " << endl; 
      Auto_Ini_Ranf(1);
    }
    if (fgr4) {
      cout << " Power spectrum C_l (r_4) --> SphereHEALPix<r_4> (float) " << endl;
      _Cl2Map<r_4>::BuildMap(infile, outfile, msph, beam, fgfitsin, fgfitsout);
    }
    else {
      cout << " Power spectrum C_l (r_8) --> SphereHEALPix<r_8> (double) " << endl;
      _Cl2Map<r_8>::BuildMap(infile, outfile, msph, beam, fgfitsin, fgfitsout);
    }
  }
  catch (PThrowable & exc) {   // Exceptions de SOPHYA
    cerr << " cl2map: Catched Exception " << (string)typeid(exc).name()
         << " - Msg= " << exc.Msg() << endl;
  }
  catch (...) {    // Autres Exceptions
    cerr << " cl2map: some other exception was caught ! " << endl;
  }
  
  PrtTim("End of cl2map ");
  return(0);
}