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

Timestamp:

Aug 2, 2010, 7:29:48 PM (15 years ago)

Author:

ansari

Message:

modifs et ajout de programme pour traitement cartes GSM (Global Sky Model), Reza 02/08/2010

File:

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

trunk/Cosmo/RadioBeam/srcat2cube.cc

@@ -72 +72 @@
     sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
     sa_size_t idxf = nvss.IndexNom("S1_4");
-    cout << " ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf << endl;
+    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);
@@ -83 +87 @@
 
     sa_size_t srccnt=0;
+    sa_size_t extendedsrccnt=0;
+
     double meanflx=0.;
     double flxmin=9.e99;
@@ -89 +95 @@
     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);
@@ -94 +102 @@
       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;
@@ -100 +110 @@
       if ((i<0)||(i>=omap.SizeX()))  continue;
       if ((j<0)||(j>=omap.SizeY()))  continue;
-      omap(i,j) += flx;
+      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;
@@ -108 +134 @@
     cout << "srccat2cube[3]: Output rectangular map computed " << endl;
     meanflx /= (double)srccnt;
-    cout << " SrcCount in map: " << srccnt << " -> meanFlx=" << meanflx << " min=" << flxmin 
+    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;