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

Timestamp:

Oct 17, 2011, 10:56:41 AM (14 years ago)

Author:

ansari

Message:

Implementatiom prise en compte dA(z) ds la calcul de bruit Pnoise(k) - Reza 17/10/2011

File:

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

trunk/Cosmo/RadioBeam/pknoise.cc

@@ -6 +6 @@
 
   Usage:  pknoise pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName
-          -parname : -noise , -renmax , -scut , -ngen , -z , -bsize , -cszmpc ,  -prt 
+          -parname : -noise , -renmax, -scut, -ngen, -z, -bsize, -cszmpc, -nbin -prt 
 ---------------------------------------------------------------  */
 
@@ -30 +30 @@
 #include "timing.h"
 #include "ctimer.h"
+#include "slininterp.h"
 
 typedef DR48RandGen RandomGenerator ;
@@ -47 +48 @@
        << "       NGen=-1 -> Call ComputeNoisePk(), dont change Four3DPk cell size \n"
        << "    -bsize sx,sy,sz : 3D real space box size (default=512x512x256) \n"   
-       << "    -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x3) \n"   
-       << "    -z z,dA,H_z : redshift,ComovDist(z),H(z) (default=1.0,3330,120.5) \n"
+       << "    -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x1.5) \n"   
+       << "    -z redshift : redshift  (default=1.0) \n"
        << "    -scut SCutValue (default= -100.) \n"
        << "       if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
-       << "    -prt PrtLev,PrtModulo (default=0,10) " << endl;
+       << "    -prt PrtLev,PrtModulo (default=0,10) \n"
+       << "    -nbin PkNBin : nomber of histogram bins for Pnoise(k) (default= 256) \n"
+       << "  NOTE : this program needs redshift_da.txt andredshift_Hz.txt files" << endl;
   return;
 }
@@ -87 +90 @@
   double H_z=120.5;  // Hubble_param(z) 
   int box3dsz[3]={512,512,256};
-  double cellsz[3]={3.5,3.5,3.5};
-
+  double cellsz[3]={3.,3.,1.5};
+  int NBINPK=256;
   int prtlev=0;
   int prtmod=10;
@@ -109 +112 @@
     }
     else if (strcmp(arg[ka],"-z")==0) {
-      sscanf(arg[ka+1],"%lg,%lg,%lg",&z_Redshift,&comov_dA_z,&H_z);  ka+=2;
+      z_Redshift=atof(arg[ka+1]);  ka+=2;
     }
     else if (strcmp(arg[ka],"-bsize")==0) {
@@ -116 +119 @@
     else if (strcmp(arg[ka],"-cszmpc")==0) {
       sscanf(arg[ka+1],"%lg,%lg,%lg",cellsz,cellsz+1,cellsz+2);   ka+=2;
+    }
+    else if (strcmp(arg[ka],"-nbin")==0) {
+      NBINPK=atoi(arg[ka+1]);    ka+=2;
     }
     else if (strcmp(arg[ka],"-prt")==0) {
@@ -147 +153 @@
     cout << " pknoise[0] : Executing, output file= " << outfile << endl;  
     POutPersist po(outfile);
+    
+    cout << " pknoise[0.b] : Creating interpolating objects from redshift_da.txt and redshift_Hz.txt " << endl;
+    SLinInterp1D redshift2da;
+    string interpfilename="redshift_da.txt";
+    redshift2da.ReadXYFromFile(interpfilename);
+    SLinInterp1D redshift2hz;
+    interpfilename="redshift_hz.txt";
+    redshift2hz.ReadXYFromFile(interpfilename);
+    comov_dA_z=redshift2da(z_Redshift);
+    H_z=redshift2hz(z_Redshift);
+    cout << " Redshift= " << z_Redshift << " -> ComovdA= " << comov_dA_z
+         << " Mpc  H(z)= " << H_z << " km/s/Mpc " << endl;
 
     H21Conversions conv;
@@ -184 +202 @@
     RandomGenerator rg;
 
+    double dfreq=H_z*cellsz[2]/(1+z_Redshift)/lambda/1000.;
+    double freq0=conv.getFrequency()-dfreq*box3dsz[2]/2;
+
     double dkxmpc=2.*DeuxPI; 
     double dkympc=2.*DeuxPI; 
     double dkzmpc=2.*DeuxPI; 
     double angscale=1.;
+    Vector angscales(box3dsz[2]);
+    angscales=angscale;
     if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[3]=128; }
     else {
       angscale=comov_dA_z;
+      for(int kz=0; kz<box3dsz[2]; kz++) angscales(kz)=redshift2da(1420.4/(freq0+(double)kz*dfreq)-1.);
       dkxmpc = DeuxPI/(double)box3dsz[0]/cellsz[0];
       dkympc = DeuxPI/(double)box3dsz[1]/cellsz[1];
       dkzmpc = DeuxPI/(double)box3dsz[2]/cellsz[2];
     }
-    cout << " pknoise[2.b]: ComovDist" << comov_dA_z << " Mpc H(z)=" << H_z << " Mpc/km/s -> angscale=" << angscale << endl; 
-    cout << " pknoise[2.c]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << "," 
+    cout << " pknoise[2.b]: ComovDist=" << comov_dA_z << " Mpc H(z)=" << H_z << " Mpc/km/s -> angscale=" << angscale << endl; 
+    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] << "," 
          << box3dsz[2] << ")" << endl;
     Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
@@ -203 +229 @@
     m3d.SetPrtLevel(prtlev,prtmod);
 
-    cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
     if (NMAX>0) {
+      cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
       PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
-      double dfreq=H_z*cellsz[2]/(1+z_Redshift)/lambda/1000.;
-      double freq0=conv.getFrequency()-dfreq*box3dsz[2]/2;
-      cout << " pknoise[3.b]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")=" 
-           << conv.getFrequency() << " AngScale=" << angscale << endl;
       pkn.SetFreqRange(freq0, dfreq);
-      pkn.SetAngScaleConversion(angscale);
+      pkn.SetAngScaleConversion(angscales);
       pkn.SetPrtLevel(prtlev,prtmod);
-      HProf hpn = pkn.Compute();
-      cout << " pknoise[3.c]: writing hpn noise profile to " << outfile << endl;
-      po << hpn;
+      HProf hpn = pkn.Compute(NBINPK);
     }
     else {
-      Histo fracmodok;
-      DataTable dtnoise;
-      HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,angscale,SCut);
-      HProf h1dnoise=arep_p->GetProjNoiseLevel();
-      HProf h1drep=arep_p->GetProjResponse();
-      cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
-      po << PPFNameTag("dtnoise") << dtnoise;
-      po << PPFNameTag("hpnoise") << hpn;
-      po << PPFNameTag("fracmodok") << fracmodok;
-      po << PPFNameTag("h1dnoise") << h1dnoise;
-      po << PPFNameTag("h1drep") << h1drep;
-      po << PPFNameTag("h2drep") << h2drep;
-    }
+      cout << " pknoise[3]: Computing Noise P(k) m3d.ComputeNoisePk(...) " << endl;
+      HProf hpn = m3d.ComputeNoisePk(*(arep_p),angscale,SCut,NBINPK);
+    }
+        
+    HProf hpnoise=m3d.GetPk();
+    DataTable dtnoise;
+    Histo fracmodok=m3d.FillPkDataTable(dtnoise);
+    HProf h1dnoise=arep_p->GetProjNoiseLevel();
+    HProf h1drep=arep_p->GetProjResponse();
+    cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
+    po << PPFNameTag("dtnoise") << dtnoise;
+    po << PPFNameTag("hpnoise") << hpnoise;
+    po << PPFNameTag("fracmodok") << fracmodok;
+    po << PPFNameTag("h1dnoise") << h1dnoise;
+    po << PPFNameTag("h1drep") << h1drep;
+    po << PPFNameTag("h2drep") << h2drep;
+  
     rc = 0;
   }  // End of try bloc