Changeset 3785 in Sophya for trunk/Cosmo/RadioBeam/syncube.cc

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

File:

• trunk/Cosmo/RadioBeam/syncube.cc (modified) (8 diffs)

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;
   }