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

/*  ------------------------ 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 nvssTocube(DataTable& nvss, TArray<r_4>& omap,  TArray<r_4>& cube);
int north20Tocube(DataTable& nor, TArray<r_4>& omap,  TArray<r_4>& cube);

//----------------------------------------------------------------------------
//----------------------------------------------------------------------------   
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 < 4) {
    cout << "Usage: srccat2cube -nvss/-north20 CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
    return 1;
  }


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

  cout << " ====== srccat2cube : Input catalog name= " << inname << " OutName=" << outname;
  bool fginmap=true;
  try {
    DataTable cat;
    cout << "srccat2cube[1]: reading source catalog from " << inname << endl;
    {
      FitsInOutFile fis(inname, FitsInOutFile::Fits_RO);
      fis >> cat;
    }
    cout << cat ;
    TArray<r_4> omap(NPhi,NTheta);
    TArray<r_4> ocube(NPhi,NTheta,NFreq);
    if (copt=="-nvss")  nvssTocube(cat, omap, ocube);
    else north20Tocube(cat, omap, ocube);
    
    {  // On sauve le cube de sortie
      cout << " srccat2cube[7]: Saving output cube to -> " << outname << endl;
      POutPersist poc(outname);
      poc << ocube;
    }
    
    if (narg > 4) {
      string ppfname = arg[4];
      cout << " srccat2cube[8]: saving 2D source map to PPF file-> " << ppfname << endl;
      POutPersist po(ppfname);
      po << omap;
    }
    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;
}



/* -- Fonction -- */
int nvssTocube(DataTable& nvss, TArray<r_4>& omap,  TArray<r_4>& ocube)
{
  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 << " NVSS catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf 
       << " MajAxis: " << idxmajax << " MajAxis: " << idxminax << endl;
  
  double tet0 = Theta0Degre; 
  double phi0 = Phi0Degre; 
  double tetmax = tet0+ThetaSizeDegre;
  double phimax = phi0+PhiSizeDegre;
  
  cout << "srccat2cube/NVSS[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;
    if (srcszarcmin<(0.5*mpixsizarcmin)) {   // Toute l'energie dans un seul pixel
      omap(i,j) += flx;
    }
    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*0.3;
          else omap(ii,jj) += flx*0.7/8.;
        }
      }
    }
    srccnt++;   meanflx+=flx;
    if (flx<flxmin) flxmin=flx;
    if (flx>flxmax) flxmax=flx;
  }
  
  cout << "srccat2cube/NVSS[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 minjy, maxjy;
  omap.MinMax(minjy, maxjy);
  MeanSigma(omap, mean, sigma);
  cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky ;  Sizes:" << endl;
  omap.Show(); 
  
  H21Conversions conv;
  conv.setRedshift(0.);
  conv.setOmegaPixDeg2(dphi*dtet);
  cout << "srccat2cube/NVSS[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad"  
       << " toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
  omap *= (r_4)conv.toKelvin(1.);
  MeanSigma(omap, mean, sigma);
  r_4 minT, maxT;
  omap.MinMax(minT, maxT);
  cout << " NVSS/ After conversion  : Mean=" << mean << " Sigma=" << sigma << " Min=" << minT 
       << " Max=" << maxT << "  Kelvin " << endl;

  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 = PLidxSrcMin+rg.Flat01()*(PLidxSrcMax-PLidxSrcMin);
      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;
      }
    }
  }
  cout << "srccat2cube/NVSS[5] data cube created from sources " << endl;
  ocube.MinMax(minT, maxT);
  MeanSigma(ocube, mean, sigma);
  cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
 
  return srccnt;
}

/* -- Fonction -- */
int north20Tocube(DataTable& nor, TArray<r_4>& omap,  TArray<r_4>& ocube)
{

  sa_size_t idxa = nor.IndexNom("ra");
  sa_size_t idxd = nor.IndexNom("dec");
  sa_size_t idxf = nor.IndexNom("flux");
  sa_size_t idxslo = nor.IndexNom("SpLO");
  sa_size_t idxshi = nor.IndexNom("SpHI");
  
  cout << " North20cm catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf 
       << " SpLO: " << idxslo << " SpHI: " << idxshi << endl;
  
  double tet0 = Theta0Degre; 
  double phi0 = Phi0Degre; 
  double tetmax = tet0+ThetaSizeDegre;
  double phimax = phi0+PhiSizeDegre;
  
  cout << "srccat2cube/North20[2]: projecting sources to map ..." << endl;
  
  sa_size_t srccnt=0;
  sa_size_t lowoksrccnt=0;
  
  double meanflx=0.;
  double flxmin=9.e99;
  double flxmax=-9.e99;
  
  double dtet = ThetaSizeDegre/(double)NTheta;
  double dphi = PhiSizeDegre/(double)NPhi;

  double infreq = 1420.; //  frequence de reference du flux des sources 
  double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
  double dfreq = FreqSizeMHz/(double)NFreq;  
  
  for (sa_size_t n=0; n<nor.NRows(); n++)  {
    r_8* pline=nor.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 tet = 90.-delta;
    double phi = alpha*360./24.;
    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;
    omap(i,j) += flx;
    srccnt++;   meanflx+=flx;
    if (flx<flxmin) flxmin=flx;
    if (flx>flxmax) flxmax=flx;
    double slo=pline[idxslo];
    if (slo<9.) {  // source detected at 80 cm 
      lowoksrccnt++;
    }
    else slo=-1.;
    for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
      double rapfreq = pow((freq0+k*dfreq)/infreq, slo);
      ocube(i,j,k) += flx*rapfreq;
    }
  }
  
  cout << "srccat2cube/North20[3]: Output rectangular map computed " << endl;
  meanflx /= (double)srccnt;
  cout << " SrcCount in map: " << srccnt << " SpLowOK=" << lowoksrccnt 
       << " -> meanFlx=" << meanflx << " min=" << flxmin 
       << " max=" << flxmax << " Jy" << endl;
  
  double mean, sigma;
  r_4 minjy, maxjy;
  omap.MinMax(minjy, maxjy);
  MeanSigma(omap, mean, sigma);
  cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;

  ocube.MinMax(minjy, maxjy);
  MeanSigma(ocube, mean, sigma);
  cout << " Cube : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;

  H21Conversions conv;
  conv.setRedshift(0.);
  conv.setOmegaPixDeg2(dphi*dtet);
  cout << "srccat2cube/North20[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad"  
       << " @1400MHz  toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
  
  // Jansky to Kelvin conversion 
  for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
    conv.setFrequency(freq0+k*dfreq);
    //    cout << " DBG* Freq= " << freq0+k*dfreq << " -> toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
    ocube(Range::all(), Range::all(), Range(k)) *= (r_4)conv.toKelvin(1.);
  }
  cout << "srccat2cube/North20[5] data cube in Kelvin computed " << endl;

  r_4 minT, maxT;
  ocube.MinMax(minT, maxT);
  MeanSigma(ocube, mean, sigma);
  cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
  
  return srccnt;
}