Changeset 4030 in Sophya for trunk/Cosmo/RadioBeam/pknoise.cc

Timestamp:

Oct 26, 2011, 7:42:21 PM (14 years ago)

Author:

ansari

Message:

Corrections papiers avec les commentaires du referee (1) - Reza 26/10/2011

File:

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

trunk/Cosmo/RadioBeam/pknoise.cc

@@ -50 +50 @@
        << "    -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x1.5) \n"   
        << "    -z redshift : redshift  (default=1.0) \n"
+       << "    -pnf : Apply frequency dependent noise factor=ComovDa^2*(1+z)^2/Hz (default= NON) \n"
        << "    -scut SCutValue (default= -100.) \n"
        << "       if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
@@ -91 +92 @@
   int box3dsz[3]={512,512,256};
   double cellsz[3]={3.,3.,3.};
+  bool fgnoisefreq=false;
   int NBINPK=256;
   int prtlev=0;
@@ -122 +124 @@
     else if (strcmp(arg[ka],"-nbin")==0) {
       NBINPK=atoi(arg[ka+1]);    ka+=2;
+    }
+    else if (strcmp(arg[ka],"-pnf")==0) {
+      fgnoisefreq=true;   ka++;
     }
     else if (strcmp(arg[ka],"-prt")==0) {
@@ -191 +196 @@
               << " DoL=" << DIAMETRE/lambda << " ) " << endl;
     Histo2D h2drep = arep_p->GetResponse();
+    Histo2D h2drepuv = arep_p->GetResponseUV();
     double repmax= h2drep.VMax();
     if (fgautoscut) {
@@ -210 +216 @@
     double angscale=1.;
     Vector angscales(box3dsz[2]);
+    double pnoise0=comov_dA_z*comov_dA_z*(1.+z_Redshift)*(1.+z_Redshift)/H_z;
+    Vector pnoisevec(box3dsz[2]);
+    pnoisevec=1.;
     angscales=angscale;
     if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[2]=128; }
@@ -222 +231 @@
     cout << " pknoise[2.c]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")=" 
          << conv.getFrequency() << " zmin=" << 1420.4/freq0-1. << " zmax=" << 1420.4/(freq0+box3dsz[2]*dfreq)-1. << endl;
-    cout << " pknoise[2.d]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << "," 
+    if(fgnoisefreq) {
+      for(int kz=0; kz<box3dsz[2]; kz++) {
+        double zreds=1420.4/(freq0+(double)kz*dfreq)-1.;
+        double comda=redshift2da(zreds);
+        pnoisevec(kz)=comda*comda*(1.+zreds)*(1.+zreds)/redshift2hz(zreds)/pnoise0;
+      }
+      cout << " pknoise[2.d]: PnoiseFactor:" << pnoisevec(0) << " ... " << pnoisevec(pnoisevec.Size()-1) << endl;
+    }
+    cout << " pknoise[2.e]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << "," 
          << box3dsz[2] << ")" << endl;
     Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
-    cout << " pknoise[2.d]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc 
+    cout << " pknoise[2.f]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc 
          << ") cell size (Mpc) : " << cellsz[0] << "x" << cellsz[1] << "x" << cellsz[2] << endl;
     m3d.SetCellSize(dkxmpc, dkympc, dkzmpc);
@@ -235 +252 @@
       pkn.SetFreqRange(freq0, dfreq);
       pkn.SetAngScaleConversion(angscales);
+      pkn.SetPNoiseFactor(pnoisevec);
       pkn.SetPrtLevel(prtlev,prtmod);
       HProf hpn = pkn.Compute(NBINPK);
@@ -255 +273 @@
     po << PPFNameTag("h1drep") << h1drep;
     po << PPFNameTag("h2drep") << h2drep;
-  
+    po << PPFNameTag("h2drepuv") << h2drepuv;
+
     rc = 0;
   }  // End of try bloc