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

Timestamp:

Oct 10, 2011, 11:31:03 PM (14 years ago)

Author:

ansari

Message:

Modif programme pknoise.cc et classes Four3DPk, PkNoiseCalculator pour calcul spectre de bruit en faisant varier la reponse instrumentale avec la frequence, Reza 10/10/2011

File:

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

trunk/Cosmo/RadioBeam/pknoise.cc

@@ -6 +6 @@
 
   Usage:  pknoise pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName
-          -parname : -noise , -renmax , -scut , -ngen , -z , -prt 
+          -parname : -noise , -renmax , -scut , -ngen , -z , -bsize , -cszmpc ,  -prt 
 ---------------------------------------------------------------  */
 
@@ -45 +45 @@
        << "    -ngen NGen (default=0) number of noise fourier amp generations \n"
        << "       NGen=0 -> Call ComputeNoisePk(), else generate Fourier Amplitudes (random) \n"
-       << "    -z redshift (default=0.7) \n"
+       << "       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"
        << "    -scut SCutValue (default= -100.) \n"
        << "       if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
@@ -80 +83 @@
   bool fgautoscut=true;
   double FacSCut=-SCut;
-  double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m  
+  double z_Redshift=1.0 ;  // z=0.7 : 21 cm at z=0.7 -> 0.357 m  
+  double comov_dA_z=3330.;  // Comoving radial distance = (1+z)dA
+  double H_z=120.5;  // Hubble_param(z) 
+  int box3dsz[3]={512,512,256};
+  double cellsz[3]={3.5,3.5,3.5};
+
   int prtlev=0;
   int prtmod=10;
@@ -101 +109 @@
     }
     else if (strcmp(arg[ka],"-z")==0) {
-      z_Redshift=atof(arg[ka+1]);  ka+=2;
+      sscanf(arg[ka+1],"%lg,%lg,%lg",&z_Redshift,&comov_dA_z,&H_z);  ka+=2;
+    }
+    else if (strcmp(arg[ka],"-bsize")==0) {
+      sscanf(arg[ka+1],"%d,%d,%d",box3dsz,box3dsz+1,box3dsz+2);   ka+=2;
+    }
+    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],"-prt")==0) {
@@ -169 +183 @@
     cout << " pknoise[2]: Instanciating object type Four3DPk  " << endl;
     RandomGenerator rg;
-    Four3DPk m3d(rg);
-    m3d.SetCellSize(2.*DeuxPI, 2.*DeuxPI, 2.*DeuxPI); 
+
+    double dkxmpc=2.*DeuxPI; 
+    double dkympc=2.*DeuxPI; 
+    double dkzmpc=2.*DeuxPI; 
+    double angscale=1.;
+    if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[3]=128; }
+    else {
+      angscale=comov_dA_z;
+      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] << "," 
+         << box3dsz[2] << ")" << endl;
+    Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
+    cout << " pknoise[2.d]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc << endl;
+    m3d.SetCellSize(dkxmpc, dkympc, dkzmpc);
     m3d.SetPrtLevel(prtlev,prtmod);
 
@@ -176 +206 @@
     if (NMAX>0) {
       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.SetPrtLevel(prtlev,prtmod);
       HProf hpn = pkn.Compute();
-      cout << " pknoise[3.b]: writing hpn noise profile to " << outfile << endl;
+      cout << " pknoise[3.c]: writing hpn noise profile to " << outfile << endl;
       po << hpn;
     }
@@ -184 +220 @@
       Histo fracmodok;
       DataTable dtnoise;
-      HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,SCut);
+      HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,angscale,SCut);
       HProf h1dnoise=arep_p->GetProjNoiseLevel();
       HProf h1drep=arep_p->GetProjResponse();