source: Sophya/trunk/SophyaProg/PrgMap/map2cl.cc@ 3044

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

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

#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 PrgMap
  \file map2cl.cc
  \brief \b map2cl: Computes power spectra (C(l)) on spherical maps.

  \verbatim

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

  map2cl : Spherical harmonics analysis - HEALPix map -> Power spectrum C_l
  Usage: map2cl [-float/-r4] [-lmax lval] [-thetacut dtdeg]
         [-fitsin] [-fitsout] InFileName OutFileName
   -float (-r4): single precision C_l and map (default = double)
   -lmax lval: Maximum value for l index (default=255)
   -thetacut dtdeg : Symmetric delta-theta cut (in degree) along equator
                    (default=0 -> no cut)
   -iter_order lval : 1,2,3,4... order of an iterative analysis,
              3rd order is usually optimal (default=0 -> standard analysis)
   -keepalm : Save Alm's as a 2-D array (Matrix) 
   -fitsout: Select the FITS format for the output vector/alm's (default PPF format)
   -fitsin : Select the FITS format for the input map (default PPF format)
   InFileName : Input file name (HEALPix map)
   OutFileName : Output file name (the C_l vector / and Alm's)

  \endverbatim
*/

/* Nouvelle-Fonction */
void Usage(bool fgerr)
{
  cout << endl;
  if (fgerr) {
    cout << " map2cl : Argument Error ! map2cl -h for usage " << endl; 
    exit(1);
  }
  else {
    cout << " map2cl : Spherical harmonics analysis - HEALPix map -> Power spectrum C_l \n" 
         << " Usage: map2cl [-float/-r4] [-lmax lval] [-thetacut dtdeg] [-iter_order lval]\n" 
         << "        [-keepalm AlmFileName] [-fitsin] [-fitsout] InFileName OutFileName \n" 
         << "   -float (-r4): single precision C_l and map (default = double) \n"
         << "   -lmax lval: Maximum value for l index (default=255)\n" 
         << "   -thetacut dtdeg : Symmetric delta-theta cut (in degree) along equator \n" 
         << "                    (default=0 -> no cut)\n"
         << "   -iter_order lval : 1,2,3,4... order of an iterative analysis\n" 
         << "       3rd order is usually optimal (default=0 -> standard analysis)\n"
         << "   -keepalm : Save Alm's as a 2-D complex array (Matrix<complex>) - only with PPF format  \n"
         << "   -fitsout: Select the FITS format for the output files (Map/Alm) (default PPF format) \n"
         << "   -fitsin : Select the FITS format for the input map (default PPF format) \n"
         << "   InFileName : Input file name (HEALPix map) \n" 
         << "   OutFileName : Output file name (the C_l vector) \n" << endl; 
    exit(0);
  }  
}

/* Nouvelle-Fonction */
template <class T> 
class _Map2Cl {
public :
static void ComputeCl(string & infile, string & outfile, int lmax, double tcut, 
                      int iterationOrder, bool fgfitsin, bool fgfitsout, bool keepalm)
{
  double deg2rad =  M_PI/180.;
  double minute2rad =  M_PI/(180.*60.);
 
  SphereHEALPix<T> sph;
  if (fgfitsin) {
    cout << "--- Reading Input FITS file " << infile << endl; 
    FitsInFile fii(infile);
    fii >> sph;
  }
  else {
    cout << "--- Reading Input PPF file " << infile << endl; 
    PInPersist ppi(infile);
    ppi >> sph;
  } 

  cout << " Input map : NbPixels= " <<  sph.NbPixels() << " NSide= " 
       << sph.SizeIndex() << " Resolution= "
       << sqrt(sph.PixSolAngle(0))/minute2rad << " Arcminutes " << endl;
  
  double ctcut = (tcut < 1.e-19) ? 0. : cos((90.-tcut)*deg2rad);
  cout << "--- Calling  DecomposeToCl() (lmax= " << lmax 
       << " cos_theta_cut= " << ctcut << " iter_order= " << iterationOrder << ") theta_cut=" << tcut << " deg"  << endl;
  // Decomposition de la carte en C_l
  SphericalTransformServer<T> sphtr;
  TVector<T> clvec;
  TMatrix< complex<T> > almmtx;
  if (keepalm) {
    Alm<T> alm(lmax);
    sphtr.DecomposeToAlm(sph, alm, lmax, ctcut, iterationOrder);
    almmtx.SetSize(alm.Lmax()+1, alm.Lmax()+1);
    for(int il=0; il<alm.Lmax()+1; il++)
      for(int im=0; im<il; im++)
        almmtx(il, im) = alm(il, im);
    clvec = alm.powerSpectrum();
  }
  else {
    clvec = sphtr.DecomposeToCl(sph, lmax, ctcut, iterationOrder);
  }
  T min, max;
  double mean, sigma;
  clvec.MinMax(min, max);
  MeanSigma(clvec, mean, sigma);
  cout << "--- Statistics on the computed C_l vector: Size=" << clvec.Size() << endl; 
  cout << " C_l.Min= " << min << " C_l.Max= " << max 
       << " C_l.Mean= " << mean << " C_l.Sigma= " << sigma << endl;
  
  if (fgfitsout) {
    FitsOutFile fio(outfile, FitsFile::clear);
    cout << "--- Writing C_l vector to Output FITS file " << outfile << endl; 
    fio << clvec;
  }
  else {
    POutPersist ppo(outfile);
    if (keepalm) {
      cout << "--- Writing C_l vector to Output PPF file " << outfile << endl; 
      ppo.PutObject(clvec, "Cl");
      cout << "--- Writing Alm Matrix to Output PPF file " << outfile << endl; 
      ppo.PutObject(almmtx, "Alm");      
    }
    else {
      cout << "--- Writing C_l vector to Output PPF file " << outfile << endl; 
      ppo << clvec;
    }
  }
}

};

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

  int ko=1;
  for (int k=1; k<narg; k++)   {
    if (strcmp(arg[k], "-lmax") == 0)  {
      if (k == narg-1) Usage(true);  // -lmax est suivi d'un argument 
      lmax = atoi(arg[k+1]);  k++;       // k++ pour sauter au suivant
    }
    else if (strcmp(arg[k], "-thetacut") == 0)  {
      if (k == narg-1) Usage(true);  // -thetacut est suivi d'un argument 
      tcut = atof(arg[k+1]);  k++;       // k++ pour sauter au suivant
    }
    else if (strcmp(arg[k], "-iter_order") == 0)  {
      if (k == narg-1) Usage(true);  // -iter_order est suivi d'un argument 
      iterationOrder = 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], "-keepalm") == 0) {
      fgkeepalm = 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 (fgr4) {
      cout << " SphereHEALPix<r_4>  --> Power spectrum C_l<r_4> (float)" << endl;
      _Map2Cl<r_4>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder, 
                              fgfitsin, fgfitsout, fgkeepalm);
    }
    else {
      cout << " SphereHEALPix<r_8>  --> Power spectrum C_l<r_8> (double)" << endl;
      _Map2Cl<r_8>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder,
                              fgfitsin, fgfitsout, fgkeepalm);
    }
  }
  catch (PThrowable & exc) {   // Exceptions de SOPHYA
    cerr << " map2cl: Catched Exception " << (string)typeid(exc).name()
         << " - Msg= " << exc.Msg() << endl;
  }
  catch (...) {    // Autres Exceptions
    cerr << " map2cl: some other exception was caught ! " << endl;
  }
  
  PrtTim("End of map2cl ");
  return(0);
}