Changeset 3787 in Sophya for trunk/Cosmo

Timestamp:

Jun 25, 2010, 12:00:30 PM (15 years ago)

Author:

ansari

Message:

Ajout classes / programmes de calcul d'effet de lobe sur les radio-sources, Reza 25/06/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• README (modified) (2 diffs)
• applobe.cc (added)
• calcpk.cc (modified) (6 diffs)
• calcpk2.cc (modified) (8 diffs)
• cubedef.h (modified) (2 diffs)
• lobe.cc (added)
• lobe.h (added)
• makefile (modified) (4 diffs)
• mdish.cc (modified) (3 diffs)
• mdish.h (modified) (2 diffs)
• specpk.cc (modified) (2 diffs)
• subtractradsrc.cmd (added)
• syncube.cc (modified) (1 diff)
• tjyk.cc (modified) (1 diff)

trunk/Cosmo/RadioBeam/README

@@ -16 +16 @@
   - plpkn.pic : spiapp script for plotting pknoise.cc program output 
 A.2/ power spectrum computation and foreground/radio source 
+  - subtractradsrc.cmd : Explanation and commands to compute various power spectra and foreground subtraction 
   - calcpk.cc : power spectrum P(k) computation from 3D mass or temperature map
   - calcpk2.cc : LSS deltaT/T P(k) computation from LSS+foreground cubes after cleaning
+  - applobe.cc : Antenna/array beam effect on 3D sky cube 
   - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
   - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
@@ -24 +26 @@
 A.3/ List of other common files :
   - specpk.h specpk.cc : 3D power spectrum computation classes  
+  - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
   - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations 
   - radutil.h radutil.cc : utilitaire de conversion a 21 cm

trunk/Cosmo/RadioBeam/calcpk.cc

@@ -19 +19 @@
 
 #include "qhist.h" 
+#include "lobe.h" 
+
 #include "histinit.h"
 #include "fftwserver.h"
@@ -33 +35 @@
 {
   if (narg<3) {
-    cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] " << endl;
+    cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] [PixNoiseLevel] " << endl;
     return 1;
   }
@@ -41 +43 @@
     string inppfname = arg[1];
     string outname = arg[2];
-    r_4 renfact=1.;
+    double renfact=1.;
     bool fgrenfact=false;
     if (narg>3) {
@@ -47 +49 @@
       fgrenfact=true;
     }
+    double pixsignoise = 0.;
+    bool fgaddnoise=false;
+    if (narg>4) {
+      pixsignoise=atof(arg[3]);
+      fgaddnoise=true;
+    }
+
     cout << "calcpk[1] : reading 3D map from file " << inppfname << endl;
     TArray<r_4> inmap;
@@ -53 +62 @@
       pin >> inmap;
     }
-    if (fgrenfact) inmap *= renfact;
+    if (fgrenfact) { 
+      cout << "  Applying RenormFactor inmap = inmap*rfact, rfact=" << renfact << endl;
+      inmap *= renfact;
+    }
     double mean, sigma;
     MeanSigma(inmap, mean, sigma);
@@ -59 +71 @@
     inmap.Show(); 
     cout << "  Mean=" << mean << " Sigma=" << sigma << endl;
+    if (fgaddnoise) { 
+      BeamEffect::AddNoise(inmap, pixsignoise);
+    }
+
     tm.Split(" After read ");
     

trunk/Cosmo/RadioBeam/calcpk2.cc

@@ -7 +7 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk 
+Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk [PixNoiseLevel]
 ---------------------------------------------------------------  */
 
@@ -22 +22 @@
 
 #include "qhist.h" 
+#include "lobe.h" 
+#include "cubedef.h" 
+
 #include "histinit.h"
 #include "fftwserver.h"
@@ -30 +33 @@
 typedef ThSDR48RandGen RandomGenerator ;
 //--- Declaration des fonctions 
-TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx);
+TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
+                            TArray<r_8>& synctemp, TArray<r_8>& specidx);
 
 //-------------------------------------------------------------------------
@@ -39 +43 @@
 {
   if (narg<6) {
-    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk  " << endl;
+    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk [PixNoiseLevel] " << endl;
     return 1;
   }
@@ -50 +54 @@
     r_4 rfacsync = atof(arg[4]);
     string outname = arg[5];
+
+    double pixsignoise = 0.;
+    bool fgaddnoise=false;
+    if (narg>6) {
+      pixsignoise=atof(arg[6]);
+      fgaddnoise=true;
+    }
+
+
     TArray<r_4> inmaplss, inmapsync;
     const char * tits[2]={"LSS", "Sync"};
@@ -68 +81 @@
       cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
     }
+
+    if (fgaddnoise) { 
+      cout << " calcpk2: adding noise to LSS input cube ... " << endl;
+      BeamEffect::AddNoise(inmaplss, pixsignoise);
+    }
+
     tm.Split(" After read ");
     TArray<r_8> synctemp, specidx; 
     cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
-    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, synctemp, specidx);
+    double freq0 = Freq0MHz;
+    double dfreq = FreqSizeMHz/(double)NFreq;
+    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, freq0, dfreq, synctemp, specidx);
     double mean, sigma;
     MeanSigma(inmap, mean, sigma);
@@ -117 +138 @@
 
 /* --Fonction-- */
-TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx)
+TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
+                            TArray<r_8>& synctemp, TArray<r_8>& specidx)
+// Inputs : maplss, mapsyc, freq0, dfreq
+// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
+// Return_Array : foreground subtracted LSS signal 
 {
   bool smo;
@@ -131 +156 @@
   Vector vlnf(maplss.SizeZ());
   int nprt = 0;
-  double freq0=840.;  // Frequence premier index en k (MHz)
-  double dfreq=1.;   // largeur en frequence de chaque plan (Mhz)  
+  // double freq0 : Frequence premier index en k (MHz)
+  // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
   for(sa_size_t i=0; i<maplss.SizeX(); i++) 
     for(sa_size_t j=0; j<maplss.SizeY(); j++)  {

trunk/Cosmo/RadioBeam/cubedef.h

@@ -6 +6 @@
 
 // Definition tailles et bornes du cube angles (X,Y), frequence (Z)
-static sa_size_t NTheta=256;
-static sa_size_t NPhi=256;
-static sa_size_t NFreq = 128;
 
-static double Theta0Degre = 60.;  // -> Delta = +30 deg
-static double Phi0Degre = 120.;   // -> alpha = 8h00
-static double ThetaSizeDegre = 60.;  // Taille de la carte en degre selon alpha (axe X)
-static double PhiSizeDegre = 60.;    // Taille de la carte en degre selon delta (axe Y)
+static sa_size_t NTheta=360;
+static sa_size_t NPhi=360;
+static sa_size_t NFreq = 200;
+
+static double Theta0Degre = 75.;  // -> Delta = +30 deg
+static double Phi0Degre = 135.;   // -> alpha = 9h00
+static double ThetaSizeDegre = 30.;  // Taille de la carte en degre selon alpha (axe X)
+static double PhiSizeDegre = 30.;    // Taille de la carte en degre selon delta (axe Y)
 
 static double Freq0MHz = 840.;       // premiere frequence 
-static double FreqSizeMHz = 128.;    // Taille de la carte en MHz (axe Z)
+static double FreqSizeMHz = 100.;    // Taille de la carte en MHz (axe Z)
 
 
@@ -31 +32 @@
 static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index  
 
+/*
+static sa_size_t NTheta=256;
+static sa_size_t NPhi=256;
+static sa_size_t NFreq = 128;
+
+static double Theta0Degre = 60.;  // -> Delta = +30 deg
+static double Phi0Degre = 120.;   // -> alpha = 8h00
+static double ThetaSizeDegre = 60.;  // Taille de la carte en degre selon alpha (axe X)
+static double PhiSizeDegre = 60.;    // Taille de la carte en degre selon delta (axe Y)
+static double Freq0MHz = 840.;       // premiere frequence 
+static double FreqSizeMHz = 100.;    // Taille de la carte en MHz (axe Z)
+*/
+
 #endif

trunk/Cosmo/RadioBeam/makefile

@@ -6 +6 @@
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : pknoise calcpk calcpk2 syncube srcat2cube tjyk
+all : pknoise calcpk calcpk2 syncube srcat2cube tjyk applobe
 
 clean : 
         rm Objs/* 
+
+PKGOLIST =  Objs/lobe.o Objs/specpk.o Objs/mdish.o Objs/qhist.o Objs/radutil.o 
+PKGHLIST =  lobe.h specpk.h mdish.h qhist.h radutil.h cubedef.h 
 
 ### les executables 
@@ -30 +33 @@
         echo 'makefile : tjyk made'
 
+applobe : Objs/applobe
+        echo 'makefile : applobe made'
+
 ### programme pknoise
-Objs/pknoise : Objs/pknoise.o Objs/specpk.o Objs/mdish.o Objs/qhist.o 
-        $(CXXLINK) -o Objs/pknoise Objs/pknoise.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
+Objs/pknoise : Objs/pknoise.o $(PKGOLIST)
+        $(CXXLINK) -o Objs/pknoise Objs/pknoise.o $(PKGOLIST)  $(SOPHYAEXTSLBLIST)
 
-Objs/pknoise.o : pknoise.cc specpk.h  mdish.h qhist.h 
+Objs/pknoise.o : pknoise.cc $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/pknoise.o pknoise.cc
 
 ### programme calcpk
-Objs/calcpk : Objs/calcpk.o Objs/specpk.o Objs/mdish.o Objs/qhist.o 
-        $(CXXLINK) -o Objs/calcpk Objs/calcpk.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
+Objs/calcpk : Objs/calcpk.o $(PKGOLIST)
+        $(CXXLINK) -o Objs/calcpk Objs/calcpk.o $(PKGOLIST)  $(SOPHYAEXTSLBLIST)
 
-Objs/calcpk.o : calcpk.cc specpk.h  mdish.h qhist.h 
+Objs/calcpk.o : calcpk.cc $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/calcpk.o calcpk.cc
 
 ### programme calcpk2
-Objs/calcpk2 : Objs/calcpk2.o Objs/specpk.o Objs/mdish.o Objs/qhist.o 
-        $(CXXLINK) -o Objs/calcpk2 Objs/calcpk2.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
+Objs/calcpk2 : Objs/calcpk2.o $(PKGOLIST)
+        $(CXXLINK) -o Objs/calcpk2 Objs/calcpk2.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
-Objs/calcpk2.o : calcpk2.cc specpk.h  mdish.h qhist.h 
+Objs/calcpk2.o : calcpk2.cc $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/calcpk2.o calcpk2.cc
 
 ### programme syncube
-Objs/syncube : Objs/syncube.o 
-        $(CXXLINK) -o Objs/syncube Objs/syncube.o $(SOPHYAEXTSLBLIST)
+Objs/syncube : Objs/syncube.o $(PKGOLIST)
+        $(CXXLINK) -o Objs/syncube Objs/syncube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
 Objs/syncube.o : syncube.cc 
@@ -59 +65 @@
 
 ### programme srcat2cube
-Objs/srcat2cube : Objs/srcat2cube.o 
-        $(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o  Objs/radutil.o $(SOPHYAEXTSLBLIST)
+Objs/srcat2cube : Objs/srcat2cube.o $(PKGOLIST)
+        $(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
-Objs/srcat2cube.o : srcat2cube.cc radutil.h 
+Objs/srcat2cube.o : srcat2cube.cc $(PKGHLIST)
         $(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)
+Objs/tjyk : Objs/tjyk.o $(PKGOLIST)
+        $(CXXLINK) -o Objs/tjyk Objs/tjyk.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
-Objs/tjyk.o : tjyk.cc radutil.h 
+Objs/tjyk.o : tjyk.cc $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/tjyk.o tjyk.cc
+
+### programme applobe
+Objs/applobe : Objs/applobe.o $(PKGOLIST)
+        $(CXXLINK) -o Objs/applobe Objs/applobe.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
+
+Objs/applobe.o : applobe.cc $(PKGHLIST)
+        $(CXXCOMPILE) -o Objs/applobe.o applobe.cc
+
 
 ### les classes / fonctions 
@@ -79 +93 @@
         $(CXXCOMPILE) -o Objs/mdish.o mdish.cc
 
+Objs/lobe.o : lobe.cc lobe.h mdish.h qhist.h 
+        $(CXXCOMPILE) -o Objs/lobe.o lobe.cc
+
 Objs/qhist.o : qhist.cc qhist.h 
         $(CXXCOMPILE) -o Objs/qhist.o qhist.cc

trunk/Cosmo/RadioBeam/mdish.cc

@@ -12 +12 @@
   : typ_(typ), dx_((dx>1.e-3)?dx:1.), dy_((dy>1.e-3)?dy:1.)
 {
+  setLambdaRef();
+  setLambda();
 }
 
@@ -17 +19 @@
 double Four2DResponse::Value(double kx, double ky)
 {
+  kx *= lambda_ratio_;
+  ky *= lambda_ratio_;
   double wk,wkx,wky;
   switch (typ_) 
@@ -62 +66 @@
 double Four2DRespTable::Value(double kx, double ky)
 {
+  kx *= lambda_ratio_;
+  ky *= lambda_ratio_;
   int_4 i,j;
   if ( (kx<=hrep_.XMin())||(kx>=hrep_.XMax()) || 

trunk/Cosmo/RadioBeam/mdish.h

@@ -34 +34 @@
   { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; return (*this); }
 
+  inline void setLambdaRef(double lambda=1.)
+  { lambdaref_ = lambda; }
+  inline void setLambda(double lambda=1.)
+  { lambda_ = lambda;  lambda_ratio_ = lambdaref_/lambda_; }
+  
   // Return the 2D response for wave vector (kx,ky)
   virtual double Value(double kx, double ky); 
@@ -45 +50 @@
   int typ_;
   double dx_, dy_; 
+  double lambdaref_, lambda_;
+  double lambda_ratio_;     // lambdaref_/lambda_;
+
 };
 

trunk/Cosmo/RadioBeam/specpk.cc

@@ -151 +151 @@
   // fourAmp represent 3-D fourier transform of a real input array. 
   // The second half of the array along Y and Z contain negative frequencies
-  double kxx, kyy, kzz;
+  double kxx, kyy, kzz, rep, amp;
   // sa_size_t is large integer type  
   for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
@@ -158 +158 @@
       kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
       for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
-        double kxx=(double)kx*dkx_;
-        double rep = resp(kxx, kyy);
+        kxx=(double)kx*dkx_;
+        rep = resp(kxx, kyy);
         if (crmask&&(kz==0))  mask(ky,kx)=((rep<1.e-8)?9.e9:(1./rep));
         if (rep<1.e-8)  fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
         else {
-          double amp = 1./sqrt(rep)/sqrt(2.);
+          amp = 1./sqrt(rep)/sqrt(2.);
           fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
         }

trunk/Cosmo/RadioBeam/syncube.cc

@@ -115 +115 @@
     double mean, sigma;
     MeanSigma(omap, mean, sigma);
-    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
+    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << " Sizes:" << endl;
     omap.Show(); 
 

trunk/Cosmo/RadioBeam/tjyk.cc

@@ -6 +6 @@
 int main(int narg, char* arg[])
 {
+  if (narg<2) {
+    cout << " Usage: tjyk redshift [resol_arcmin] " << endl;
+    return 1;
+  }
 
   H21Conversions conv;