source: Sophya/trunk/Cosmo/RadioBeam/srcat2cube.cc@ 3826

/*  ------------------------ Projet BAORadio -------------------- 
  Programme de fabrication d'un cube 3D (angles,fre) 
  a partir du catalogue de source radio (NVSS) 
    R. Ansari , C. Magneville - Juin 2010 

  Usage: srccat2cube CatalogFitsName Out3DPPFName [Out2DMapName]
---------------------------------------------------------------  */

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

#include "array.h"
#include "histats.h"

#include "swfitsdtable.h"
#include "fitshdtable.h"

#include "randr48.h"

#include "xastropack.h"     // Pour faire les conversions de coordonnees celestes

#include "radutil.h"

// Pour l'initialisation des modules 
#include "tarrinit.h"
#include "histinit.h"
#include "fiosinit.h"     

#include "timing.h"
#include "ctimer.h"

#include "cubedef.h"

//----------------------------------------------------------------------------
//----------------------------------------------------------------------------   
int main(int narg, char* arg[])
{
  // Sophya modules initialization
  TArrayInitiator  _inia;
  HiStatsInitiator  _inih;
  FitsIOServerInitiator  _inif;
  //------- AU LIEU DE ------>  SophyaInit();  

  InitTim();   // Initializing the CPU timer
  Timer tm("srcat2cube");

  if (narg < 3) {
    cout << "Usage: srccat2cube NVSS_CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
    return 1;
  }


  // decodage arguments 
  string outname = arg[2];
  string inname = arg[1];
  int rc = 91;

  cout << " ====== srccat2cube :   Input NVSS catalog name= " << inname << " OutName=" << outname;
  bool fginmap=true;
  try {
    DataTable nvss;
    cout << "srccat2cube[1]: reading NVSS catalog from " << inname << endl;
    {
      FitsInOutFile fis(inname, FitsInOutFile::Fits_RO);
      fis >> nvss;
    }
    cout << nvss;
    sa_size_t idxa = nvss.IndexNom("C_RAJ2000");
    sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
    sa_size_t idxf = nvss.IndexNom("S1_4");
    sa_size_t idxmajax = nvss.IndexNom("MajAxis");
    sa_size_t idxminax = nvss.IndexNom("MinAxis");

    cout << " ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf 
         << " MajAxis: " << idxmajax << " MajAxis: " << idxminax << endl;

    TArray<r_4> omap(NPhi,NTheta);
    double tet0 = Theta0Degre; 
    double phi0 = Phi0Degre; 
    double tetmax = tet0+ThetaSizeDegre;
    double phimax = phi0+PhiSizeDegre;

    cout << "srccat2cube[2]: projecting sources to map ..." << endl;

    sa_size_t srccnt=0;
    sa_size_t extendedsrccnt=0;

    double meanflx=0.;
    double flxmin=9.e99;
    double flxmax=-9.e99;

    double dtet = ThetaSizeDegre/(double)NTheta;
    double dphi = PhiSizeDegre/(double)NPhi;
    double mpixsizarcmin = 0.5*(dtet+dphi)*60.;

    for (sa_size_t n=0; n<nvss.NRows(); n++)  {
      r_8* pline=nvss.GetLineD(n);
      double alpha=pline[idxa];  // alpha en degre
      double delta=pline[idxd];  // delta en degre 
      double flx=pline[idxf]*1.e-3;  // flux en Jy
      double srcszarcmin=0.5*(pline[idxmajax]+pline[idxminax])/60.;  // taille (extension de la source en arcmin 
      if (srcszarcmin<1.) srcszarcmin=1.;
      double tet = 90.-delta;
      double phi = alpha;
      sa_size_t i = (phi-phi0)/dphi;
      sa_size_t j = (tet-tet0)/dtet;
      if ((i<0)||(i>=omap.SizeX()))  continue;
      if ((j<0)||(j>=omap.SizeY()))  continue;
      double srat = (4.*srcszarcmin*srcszarcmin)/(mpixsizarcmin*mpixsizarcmin);
      if (srcszarcmin<(0.5*mpixsizarcmin)) {   // Toute l'energie dans un seul pixel
        omap(i,j) += flx*srat;
      }
      else {  // on repartit l'energie de la source dans plusieurs pixels
        extendedsrccnt++;
        for(int bi=-1;bi<=1;bi++) {
          for(sa_size_t bj=-1; bj<=1; bj++) {
            sa_size_t ii = (phi-phi0+bi*srcszarcmin/60.)/dphi;
            sa_size_t jj = (tet-tet0+bj*srcszarcmin/60.)/dtet;
            if ((ii<0)||(ii>=omap.SizeX()))  continue;
            if ((jj<0)||(jj>=omap.SizeY()))  continue;      
            if ((bi==0)&&(bj==0))  omap(ii,jj) += flx*srat*0.3;
            else omap(ii,jj) += flx*srat*0.7/8.;
          }
        }
      }
      srccnt++;   meanflx+=flx;
      if (flx<flxmin) flxmin=flx;
      if (flx>flxmax) flxmax=flx;
    }

    cout << "srccat2cube[3]: Output rectangular map computed " << endl;
    meanflx /= (double)srccnt;
    cout << " SrcCount in map: " << srccnt << " extended=" << extendedsrccnt 
         << " -> meanFlx=" << meanflx << " min=" << flxmin 
         << " max=" << flxmax << " Jy" << endl;

    double mean, sigma;
    r_4 mintemp, maxtemp;
    omap.MinMax(mintemp, maxtemp);
    MeanSigma(omap, mean, sigma);
    cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
    omap.Show(); 

    H21Conversions conv;
    conv.setRedshift(0.);
    conv.setOmegaPixDeg2(dphi*dtet);
    cout << "srccat2cube[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad"  
         << " toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
    omap *= (r_4)conv.toKelvin(1.);
    MeanSigma(omap, mean, sigma);
    cout << " After conversion  : Mean=" << mean << " Sigma=" << sigma << "  Kelvin " << endl;
    
    if (narg > 3) {
      string ppfname = arg[3];
      cout << " srccat2cube[4]: Saving inmap/outmap tp PPF file-> " << ppfname << endl;
      POutPersist po(ppfname);
      po << PPFNameTag("omap") << omap;
    }

    TArray<r_4> ocube(NPhi,NTheta,NFreq);

    double infreq = 1420.; //  frequence de reference du flux des sources 
    double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
    double dfreq = FreqSizeMHz/(double)NFreq;  

    ThSDR48RandGen rg;
    for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
      for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
        double freqexpo = rg.Gaussian(sigPLidxSrc,PLidxSrc);
        for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
          double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
          ocube(i,j,k) = AmpPL1*omap(i,j)*rapfreq;
        }
      }
    }
    
    // On sauve le cube de sortie
    {
      cout << " srccat2cube[5]: Saving output cube to -> " << outname << endl;
      POutPersist poc(outname);
      poc << ocube;
    }

    rc = 0;
  }
  catch (PThrowable& exc) {
    cerr << " srccat2cube.cc catched Exception " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n srccat2cube.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
  }
  catch (...) {
    cerr << " srccat2cube.cc catched unknown (...) exception  " << endl; 
    rc = 78; 
  } 

  cout << ">>>> srccat2cube[9] ------- FIN ----------- Rc=" << rc << endl;
  return rc;
}