Changeset 3829 in Sophya for trunk/Cosmo/RadioBeam/srcat2cube.cc

Timestamp:

Aug 3, 2010, 6:38:00 PM (15 years ago)

Author:

ansari

Message:

Adaptation programme srcat2cube.cc pour utilisation catalogue North20cm avec indice spectral, Reza 03/08/2010

File:

• trunk/Cosmo/RadioBeam/srcat2cube.cc (modified) (4 diffs)

trunk/Cosmo/RadioBeam/srcat2cube.cc

@@ -34 +34 @@
 
 #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);
 
 //----------------------------------------------------------------------------
@@ -48 +52 @@
   Timer tm("srcat2cube");
 
-  if (narg < 3) {
-    cout << "Usage: srccat2cube NVSS_CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
+  if (narg < 4) {
+    cout << "Usage: srccat2cube -nvss/-north20 CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
     return 1;
   }
@@ -55 +59 @@
 
   // decodage arguments 
-  string outname = arg[2];
-  string inname = arg[1];
+  
+  string copt=arg[1];
+  string outname=arg[3];
+  string inname=arg[2];
   int rc = 91;
 
-  cout << " ====== srccat2cube :   Input NVSS catalog name= " << inname << " OutName=" << outname;
+  cout << " ====== srccat2cube : Input catalog name= " << inname << " OutName=" << outname;
   bool fginmap=true;
   try {
-    DataTable nvss;
-    cout << "srccat2cube[1]: reading NVSS catalog from " << inname << endl;
+    DataTable cat;
+    cout << "srccat2cube[1]: reading source 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;
-
+      fis >> cat;
+    }
+    cout << cat ;
     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;
+    TArray<r_4> ocube(NPhi,NTheta,NFreq);
+    if (copt=="-nvss")  nvssTocube(cat, omap, ocube);
+    else north20Tocube(cat, omap, ocube);
     
-    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;
+    {  // 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;
   }
@@ -206 +113 @@
 
 
+
+/* -- 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 = 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;
+      }
+    }
+  }
+  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=5.;
+    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;
+}