Changeset 3796 in Sophya for trunk

Timestamp:

Jun 29, 2010, 7:26:01 PM (15 years ago)

Author:

ansari

Message:

Suite modif utilisation de reponse calcule (tabulee) de l'interferometre pour soustraction avant-plans, Reza 29/06/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• applobe.cc (modified) (4 diffs)
• calcpk2.cc (modified) (4 diffs)
• cubedef.h (modified) (1 diff)
• mdish.cc (modified) (2 diffs)
• mdish.h (modified) (2 diffs)
• repicon.cc (modified) (5 diffs)
• subtractradsrc.cmd (modified) (2 diffs)

trunk/Cosmo/RadioBeam/applobe.cc

@@ -3 +3 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-  Usage: applobe In3DPPFName Out3DPPFName [ResmapleFactor=0.5,0.333...] [TargetBeamDoL]
+  Usage: applobe Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName 
+                 [TargetBeamDoL] [ResmapleFactor=0.5,0.333...] 
 ---------------------------------------------------------------  */
 
@@ -50 +51 @@
 
   if (narg < 3) {
-    cout << "Usage: applobe In3DPPFName Out3DPPFName [ResmapleFactor=0.5,0.333...] [TargetBeamDoL]\n" << endl;
+    cout << "Usage: applobe Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName " 
+         << "             [TargetBeamDoL] [ResmapleFactor=0.5,0.333...] \n" << endl;
     return 1;
   }
 
   // decodage arguments 
-  string outname = arg[2];
-  string inname = arg[1];
-  bool fgresample=false;
-  double fresamp=1.;
-  if (narg>3) {
-    fresamp=atof(arg[3]);
-    if ((fabs(fresamp)>1.e-2)&&(fabs(fresamp-1.)>1.e-2))  fgresample=true; 
+  bool fgresptbl=true;
+  double DIAMETRE=100.;
+  string resptblname;
+  if (isdigit(*arg[1])) {
+    fgresptbl=false;
+    DIAMETRE=atof(arg[1]);
   }
+  else resptblname=arg[1];
+
+  string inname = arg[2];
+  string outname = arg[3];
   bool fgcorrbeam=false;
   double tbeamDoL=135;
@@ -68 +73 @@
     tbeamDoL=atof(arg[4]);
     if (tbeamDoL>1)  fgcorrbeam=true;
+  }
+  bool fgresample=false;
+  double fresamp=1.;
+  if (narg>5) {
+    fresamp=atof(arg[5]);
+    if ((fabs(fresamp)>1.e-2)&&(fabs(fresamp-1.)>1.e-2))  fgresample=true; 
   }
 
@@ -101 +112 @@
     conv.setFrequency(Freq0MHz);
     double lambda = conv.getLambda();
-    Four2DResponse fresp(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda, lambda);    
-    BeamEffect beam(fresp);
+    
+    Four2DResponse fresp(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
+    Four2DResponse* fresp_p=&fresp;
+    Four2DRespTable resptbl;
+    if (fgresptbl) {
+      cout << "applobe[2.b]: initializing Four2DRespTable from file" << resptblname << endl;
+      resptbl.readFromPPF(resptblname);
+      fresp_p=&resptbl;
+    }
+    else cout << " applobe[2.b]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
+              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
+    BeamEffect beam(*fresp_p);
 
     if (fgcorrbeam) {
       double DoL = tbeamDoL;
       double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
-      Four2DResponse tbeam(2, DoL, DoL );
+      int typcb = 2;
+      if (fgresptbl) typcb=22;
+      Four2DResponse tbeam(typcb, DoL, DoL );
       cout << "applobe[3]: calling Correct2RefLobe() with target beam D/Lambda=" << DoL  
-           << " -> arcmin " << tbeamarcmin << endl;
+           << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
       beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq);
     }

trunk/Cosmo/RadioBeam/calcpk2.cc

@@ -7 +7 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk 
-               [PixNoiseLevel] [TargetBeamArcmin] [NSigSrcThr]
+Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPkFile 
+               [PixNoiseLevel] [Diameter/Four2DRespTableFile] [TargetBeamArcmin] [NSigSrcThr]
 ---------------------------------------------------------------  */
 
@@ -44 +44 @@
 {
   if (narg<6) {
-    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk \n" 
-         << "        [PixNoiseLevel] [TargetBeamArcmin] [NSigSrcThr]" << endl;
+    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPkFile \n" 
+         << "        [PixNoiseLevel] [Diameter/Four2DRespTableFile] [TargetBeamArcmin] [NSigSrcThr]" << endl;
     return 1;
   }
@@ -65 +65 @@
     
     bool fgcorrbeam=true;
+
+    bool fgresptbl=false;
+    double DIAMETRE=100.;
+    string resptblname;
+    if (narg>7) {
+      if (isdigit(*arg[7])) {
+        fgresptbl=false;
+        DIAMETRE=atof(arg[7]);
+      }
+      else { 
+        resptblname=arg[7];
+        fgresptbl=true;
+      }
+    }
     double tbeamDoL=135;
-    if (narg>7) {
-      tbeamDoL=atof(arg[7]);
+    if (narg>8) {
+      tbeamDoL=atof(arg[8]);
       if (tbeamDoL<1.)  fgcorrbeam=false;
     }
     bool fgclnsrc=true;
     double nsigsrc=5.;
-    if (narg>8) {
-      nsigsrc=atof(arg[8]);
+    if (narg>9) {
+      nsigsrc=atof(arg[9]);
       if (nsigsrc<1.e-6)  fgclnsrc=false;
     }
@@ -119 +133 @@
     conv.setFrequency(Freq0MHz);
     double lambda = conv.getLambda();
-    Four2DResponse arep(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda, lambda);
+    Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
+    Four2DResponse* arep_p=&arep;
+    Four2DRespTable resptbl;
+    if (fgresptbl) {
+      cout << "calcpk2[3.a]: initializing Four2DRespTable from file" << resptblname << endl;
+      resptbl.readFromPPF(resptblname);
+      arep_p=&resptbl;
+    }
+    else cout << " calcpk2[3.a]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
+              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
+
     double DoL = tbeamDoL;
     double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
-    Four2DResponse tbeam(2, DoL, DoL );
-    ForegroundCleaner  cleaner(arep, tbeam, skycube); 
+    int typcb = 2;
+    //    if (fgresptbl) typcb=22;
+    Four2DResponse tbeam(typcb, DoL, DoL );
+
+    ForegroundCleaner  cleaner(*arep_p, tbeam, skycube); 
     if (fgcorrbeam) {
-      cout << "calcpk2[3.a] : calling cleaner.BeamCorrections() for target beam D/Lambda=" << DoL 
-           << "  -> arcmin " << tbeamarcmin << endl;
+      cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam D/Lambda=" << DoL 
+           << "  -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
       cleaner.BeamCorrections();
     }
-    cout << " calcpk2[3.b] : calling cleaner.CleanNegatives() ... " << endl;
+    cout << " calcpk2[3.c] : calling cleaner.CleanNegatives() ... " << endl;
     cleaner.CleanNegatives();
     if (fgclnsrc) {
-      cout << "calcpk2[3.c] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
+      cout << "calcpk2[3.d] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
       cleaner.CleanPointSources(nsigsrc);
     }

trunk/Cosmo/RadioBeam/cubedef.h

@@ -39 +39 @@
 static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index  
 
-// Parametres pour la reponse en Fourier de l'instrument :
-static double InterfArrayDiametre = 50.;  // Taille du reseau en metres 
-
 /*
 static sa_size_t NTheta=256;

trunk/Cosmo/RadioBeam/mdish.cc

@@ -33 +33 @@
       return ( (wk<1.)?(1.-wk):0.);
       break;
+    case 22:   // Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda + trou au centre
+      wk = sqrt(kx*kx+ky*ky)/dx_/2./M_PI;
+      if (wk<0.025) return 39.*wk;
+      else if (wk<1.) return (1.-wk);
+      else return 0.;
+      break;
     case 3:   // Reponse rectangle Dx x Dy  Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda) 
-      wkx = kx/2./M_PI/dx_; 
-      wky = ky/2./M_PI/dy_; 
+      wkx = fabs(kx)/2./M_PI/dx_; 
+      wky = fabs(ky)/2./M_PI/dy_; 
       return ( ((wkx<1.)&&(wky<1.))?((1.-wkx)*(1-wky)):0.);
       break;
@@ -78 +84 @@
   hrep_.FindBin(kx, ky, i, j);
   return hrep_(i, j);
+}
+
+double Four2DRespTable::renormalize(double max)
+{
+  double cmx = hrep_.VMax();
+  hrep_ *= (max/cmx);
+  return cmx;
 }
 

trunk/Cosmo/RadioBeam/mdish.h

@@ -20 +20 @@
 #endif
 
-// -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (angles theta,phi)
+// -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (frequences angulaires theta,phi)
 // typ=1 : Reponse gaussienne parabole diametre D exp[ - 0.5  (lambda  k_g / D )^2 ]
 // typ=2 : Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda   
 // typ=3 : Reponse rectangle Dx x Dy  Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda) 
+// typ=22 : Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda  avec un trou au centre  
+
 class Four2DResponse  {
 public:
@@ -65 +67 @@
   // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
   Four2DRespTable(Histo2D const & hrep, double d, double lambda=1.);
+  // Apres renormalisaton Value(kx,ky) <= max
+  double renormalize(double max=1.);
   // Return the 2D response for wave vector (kx,ky)
   virtual double Value(double kx, double ky); 

trunk/Cosmo/RadioBeam/repicon.cc

@@ -36 +36 @@
 {
   if (((narg>1)&&(strcmp(arg[1],"-h")==0))||(narg<3))  {
-    cout << " Usage: repicon configId OutPPFName [Lambda] " 
+    cout << " Usage: repicon configId OutPPFName [Lambda=0.357] [RenormalizeMax]" 
          << " configs: f4x4 , f8x8 , f20x20 Filled array of nxn dishes \n" 
          << "          confA , confB, confC : semi-filled array of dishes \n"
@@ -54 +54 @@
   if (outfile==".")  outfile = "repicon.ppf";
   double LAMBDA=0.357 ;  // 21 cm at z=0.7 
-  if (narg>3) LAMBDA = atof(arg[2]);
-
+  if (narg>3) LAMBDA = atof(arg[3]);
+  bool fgrenorm=false;
+  double rmax=1.;
+  if (narg>4) {
+    rmax=atof(arg[4]);
+    fgrenorm=true;
+  }
   //-- end command line arguments
   
@@ -61 +66 @@
   try {  // exception handling try bloc at top level
 
-    cout << " 2.c/ Creating MultiDish Filled Array " << endl;
     double Ddish = 5.;
     double Ddish2 = 7.5;
@@ -122 +126 @@
     double LMAX = D;
     bool fgnoauto = true;
-    int NRX=100;
-    int NRY=100;
+    int NRX=200;
+    int NRY=200;
 
     MultiDish mdish(LAMBDA, LMAX, vdishes, fgnoauto);
@@ -134 +138 @@
     PrtTim("Apres mdish.GetResponse()");
 
+    Four2DRespTable mdresp(hrep, Dol, LAMBDA);
+    if (fgrenorm) {
+      cout << " repicon[2.b] call to mdresp.renormalize(" << rmax << ")"; 
+      double omax=mdresp.renormalize(rmax);
+      cout << " Old Max=" << omax << endl;
+    }
     cout << " repicon[3] : saving Four2DRespTable for config " << config << " to PPF file " << outfile << endl;
-    Four2DRespTable mdresp(hrep, Dol, LAMBDA);
     mdresp.writeToPPF(outfile);
 

trunk/Cosmo/RadioBeam/subtractradsrc.cmd

@@ -50 +50 @@
 csh> spiapp -term -exec sumcubes.pic 
 
-## Step 3/ Apply lobe effect on foreground cube and LSS cube
-csh> ./Objs/applobe fgndcube.ppf fgndcube_lobe.ppf
-csh> ./Objs/applobe lsscube.ppf lsscube_lobe.ppf
+## Step 3/ Apply lobe (50 meter diameter array) effect on foreground cube and LSS cube
+csh> ./Objs/applobe 50. fgndcube.ppf fgndcube_lobe.ppf
+csh> ./Objs/applobe 50. lsscube.ppf lsscube_lobe.ppf
 ## Step 3.b/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 130  
-csh> ./Objs/applobe lsscube_lobe.ppf lsscube_corlobe.ppf 1 130
-csh> ./Objs/applobe fgndcube_lobe.ppf fgndcube_corlobe.ppf 1 130
+csh> ./Objs/applobe 50. lsscube_lobe.ppf lsscube_corlobe.ppf 130
+csh> ./Objs/applobe 50. fgndcube_lobe.ppf fgndcube_corlobe.ppf 130
 
 ### Step 4/ Compute power spectra 
@@ -77 +77 @@
 
 # 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam 
-csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 0. 0.
+csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 0. 0.
 # 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130 
-csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 130. 3.
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 130. 3.
 
 ### Step 5 / Check the results using spiapp