Changeset 3785 in Sophya for trunk/Cosmo

Timestamp:

Jun 16, 2010, 11:53:54 PM (15 years ago)

Author:

ansari

Message:

Ajout programme srcat2cube.cc de calcul de carte des radio source, Reza 16/06/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• README (modified) (2 diffs)
• cubedef.h (added)
• makefile (modified) (3 diffs)
• srcat2cube.cc (added)
• syncube.cc (modified) (8 diffs)
• tjyk.cc (modified) (2 diffs)

trunk/Cosmo/RadioBeam/README

@@ -19 +19 @@
   - calcpk2.cc : LSS deltaT/T P(k) computation from LSS+foreground cubes after cleaning
   - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
+  - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
   - tjyk.cc : test de la classe H21Conversions 
 
@@ -25 +26 @@
   - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations 
   - radutil.h radutil.cc : utilitaire de conversion a 21 cm
+  - cubedef.h : size and limits for 3D temperature cubes generated by  syncube.cc , srcat2cube.cc
  
 E/ P(k) in temperature plot + "cosmic variance" + noise ...

trunk/Cosmo/RadioBeam/makefile

@@ -6 +6 @@
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : pknoise calcpk calcpk2 syncube tjyk
+all : pknoise calcpk calcpk2 syncube srcat2cube tjyk
 
 clean : 
@@ -23 +23 @@
 syncube : Objs/syncube
         echo 'makefile : syncube made'
+
+srcat2cube : Objs/srcat2cube
+        echo 'makefile : srcat2cube made'
 
 tjyk : Objs/tjyk
@@ -55 +58 @@
         $(CXXCOMPILE) -o Objs/syncube.o syncube.cc
 
-### programme syncube
+### programme srcat2cube
+Objs/srcat2cube : Objs/srcat2cube.o 
+        $(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o  Objs/radutil.o $(SOPHYAEXTSLBLIST)
+
+Objs/srcat2cube.o : srcat2cube.cc radutil.h 
+        $(CXXCOMPILE) -o Objs/srcat2cube.o srcat2cube.cc
+
+### programme tjyk
 Objs/tjyk : Objs/tjyk.o Objs/radutil.o 
         $(CXXLINK) -o Objs/tjyk Objs/tjyk.o Objs/radutil.o $(SOPHYAEXTSLBLIST)

trunk/Cosmo/RadioBeam/syncube.cc

@@ -44 +44 @@
 #include "fiosinit.h"     
 
-
 #include "timing.h"
 #include "ctimer.h"
+
+#include "cubedef.h"
 
 //----------------------------------------------------------------------------
@@ -65 +66 @@
     return 1;
   }
-
-  //---  Parametres des lois de puissance en frequence 
-  double AmpPL1 = 1.;   // amp max PowerLaw 2 
-  double PLidx1 = -2.5;  // index de la loi de puissance synchrotron 
-  double sigPLidx1 = 0.1;  // Sigma de la variation (gaussienne) de index1 
-  // Amplitude max de la 2eme composante en loi de puissance (tirage plat 0 ... AmpPL2)
-  double AmpPL2 = 0.1;   // amp max PowerLaw 2 
-  double PLidx2 = -3.2; 
-  double sigPLidx2 = 0.15;
 
   // decodage arguments 
@@ -101 +93 @@
     cout << "syncube[1]: Input resolution doubled : " << inmap << endl;
 
-    sa_size_t NTet=256;
-    sa_size_t NPhi=256;
-    LocalMap<r_4> outmap(NPhi,NTet,60.,60.,110.,150.);
+   LocalMap<r_4> outmap(NPhi,NTheta,60.,60.,110.,150.);
     for(int_4 kk=0; kk<outmap.NbPixels(); kk++) {
       double theta, phi;  // Theta, Phi en radians
@@ -110 +100 @@
     }    
 
-    TArray<r_4> omap(NPhi,NTet);
-    double tet0 = Angle(60.,Angle::Degree).ToRadian();
-    double phi0 = Angle(120.,Angle::Degree).ToRadian();
-    double dtet = Angle(60.,Angle::Degree).ToRadian()/(double)NTet;
-    double dphi = Angle(60.,Angle::Degree).ToRadian()/(double)NPhi;
+    TArray<r_4> omap(NPhi,NTheta);
+    double tet0 = Angle(Theta0Degre,Angle::Degree).ToRadian();
+    double phi0 = Angle(Phi0Degre,Angle::Degree).ToRadian();
+    double dtet = Angle(ThetaSizeDegre,Angle::Degree).ToRadian()/(double)NTheta;
+    double dphi = Angle(PhiSizeDegre,Angle::Degree).ToRadian()/(double)NPhi;
     for (sa_size_t j=0; j<omap.SizeY(); j++)  {
       double theta = j*dtet+tet0;
@@ -122 +112 @@
       }
     }
+
+    double mean, sigma;
+    MeanSigma(omap, mean, sigma);
+    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
+    omap.Show(); 
+
     //    Sph2Sph(inmap,outmap);
     
-    cout << "syncube[2]: Output rectangular map computed " << outmap << endl;
     if (narg > 3) {
       string ppfname = arg[3];
@@ -134 +129 @@
     }
 
-    sa_size_t NFreq = 128;
-    TArray<r_4> ocube(NPhi,NTet,NFreq);
+    TArray<r_4> ocube(NPhi,NTheta,NFreq);
 
     double infreq = 400.; // Frequence carte input en MHz
-    double freq0 = 840.;  // Freq0 du cube de sortie 
-    double dfreq = 1.;  
+    double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
+    double dfreq = FreqSizeMHz/(double)NFreq;  
 
     ThSDR48RandGen rg;
@@ -162 +156 @@
     // On sauve le cube de sortie
     {
-      cout << " syncube[2]: Saving output cube to -> " << outname << endl;
+      cout << " syncube[4]: Saving output cube to -> " << outname << endl;
       POutPersist poc(outname);
       poc << ocube;
@@ -170 +164 @@
   }
   catch (PThrowable& exc) {
-    cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
+    cerr << " syncube.cc catched Exception " << exc.Msg() << endl;
     rc = 77;
   }  

trunk/Cosmo/RadioBeam/tjyk.cc

@@ -2 +2 @@
 #include <iostream>
 
+// Appel : tjyk [z=redshift] [resol_arcmin]
+ 
 int main(int narg, char* arg[])
 {
@@ -7 +9 @@
   H21Conversions conv;
   
-  conv.setRedshift(0.6);
-  conv.setOmegaPixArcmin2(10.*10.);
-  cout << " --------------------------------------------------------------- " << endl;
+  double z = 0.;
+  double resol = 10.;
+  if (narg>1) z = atof(arg[1]);
+  if (narg>2) resol = atof(arg[2]);
+
+  conv.setRedshift(z);
+  conv.setOmegaPixArcmin2(resol*resol);
+  cout << " ------------------------------ H21Conversions ------------------------- " << endl;
+  cout << "  Redshift=z= " << z << " Resol(arcmin)= " << resol << endl;
   cout << " H21Conversions: z=" << conv.getRedshift() << " Freq=" << conv.getFrequency()
        << " MHz , Lambda=" << conv.getLambda() << " m , OmegaPix=" << conv.getOmegaPix() << " srad" << endl;