Changeset 4026 in Sophya for trunk/Cosmo

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

Location:

trunk/Cosmo/RadioBeam

Files:

• README (modified) (1 diff)
• pknoise.cc (modified) (7 diffs)
• specpk.cc (modified) (8 diffs)
• specpk.h (modified) (4 diffs)

trunk/Cosmo/RadioBeam/README

@@ -24 +24 @@
   - gsm2cube.cc : Computes 3D temperature cube from a set of GSM (Global Sky Model) healpix PPF maps  
   - tjyk.cc : test de la classe H21Conversions 
+  - ckpowl.cc : Programme de verification comportement en loi de puissance 
 
 A.3/ List of other common files :

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();

trunk/Cosmo/RadioBeam/specpk.cc

@@ -6 +6 @@
 
 #include "specpk.h"
+#include "radutil.h"
 #include "randr48.h"      
 #include "ctimer.h"      
@@ -162 +163 @@
 }
 
+
 // Generate mass field Fourier Coefficient
-void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, bool crmask)
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale, bool crmask)
+// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
+// typically = ComovRadialDistance 
 {
   TMatrix<r_4> mask(fourAmp.SizeY(), fourAmp.SizeX());
@@ -176 +180 @@
       for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
         kxx=(double)kx*dkx_;
-        rep = resp(kxx, kyy);
+        rep = resp(kxx*angscale, kyy*angscale);
         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.);
@@ -193 +197 @@
   if (prtlev_>0)
     cout << " Four3DPk::ComputeNoiseFourierAmp() done ..." << endl;
+}
+
+// Generate mass field Fourier Coefficient
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, double angscale)
+// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
+// typically = ComovRadialDistance 
+{
+  H21Conversions conv;
+  // 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, rep, amp;
+  // sa_size_t is large integer type  
+  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
+    conv.setFrequency(f0+kz*df);
+    resp.setLambda(conv.getLambda());
+    for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+      kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
+      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+        kxx=(double)kx*dkx_;
+        rep = resp(kxx*angscale, kyy*angscale);
+        if (rep<1.e-8)  fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
+        else {
+          amp = 1./sqrt(rep)/sqrt(2.);
+          fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
+        }
+      }
+    }
+  }
+
+  if (prtlev_>1)
+    cout << " Four3DPk::ComputeNoiseFourierAmp(...) Computing FFT along frequency ..." << endl;
+  TVector< complex<TF> >  veczin(fourAmp.SizeZ()), veczout(fourAmp.SizeZ());
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+
+  for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+    for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+      //      veczin=fourAmp(Range(kx), Range(ky), Range::all());
+      for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++)  veczin(kz)=fourAmp(kx,ky,kz);
+      ffts.FFTBackward(veczin,veczout);
+      veczout /= (TF)sqrt((double)fourAmp.SizeZ());
+      //      fourAmp(Range(kx), Range(ky), Range::all())=veczout;
+      for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++)  fourAmp(kx,ky,kz)=veczout(kz);
+    }
+  }
+
+  if (prtlev_>2)  fourAmp.Show();
+  if (prtlev_>0)
+    cout << " Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df) done ..." << endl;
 }
 
@@ -270 +323 @@
 // cells with amp^2=re^2+im^2>s2cut are ignored
 // Output : noise power spectrum (profile histogram)
+// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
+// typically = ComovRadialDistance 
 HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt, 
-                               double s2cut, int nbin, double kmin, double kmax)
+                               double angscale, double s2cut, int nbin, double kmin, double kmax)
 {
   // The second half of the array along Y (matrix rows) contain
@@ -305 +360 @@
         double kxx=(double)kx*dkx_;
         double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
-        double amp2 = 1./resp(kxx, kyy); 
+        double rep=resp(kxx*angscale, kyy*angscale);
+        double amp2 = (rep>1.e-19)?1./rep:1.e19; 
         hmcnt.Add(wk);
         if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
@@ -349 +405 @@
   : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
 {
+  SetFreqRange();
+  SetAngScaleConversion();
   SetPrtLevel();
 }
@@ -357 +415 @@
   tm.Nop();
   HProf hnd;
-  cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
+  cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT  
+       << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscale=" << angscale_ << endl;
   for(int igen=0; igen<NGEN; igen++) {
-    pkn3d.ComputeNoiseFourierAmp(frep);
+    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscale_);
     if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
     else pkn3d.ComputePkCumul(hnd,S2CUT);

trunk/Cosmo/RadioBeam/specpk.h

@@ -51 +51 @@
   { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
   void ComputeFourierAmp(SpectralShape& pk);
-  void ComputeNoiseFourierAmp(Four2DResponse& resp, bool crmask=false);
+// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
+// typically = ComovRadialDistance 
+  void ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale=1., bool crmask=false);
+  void ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, double angscale=1.);
+
 // Return the array size 
   inline sa_size_t NCells() { return fourAmp.Size(); }
@@ -65 +69 @@
   void  ComputePkCumul(HProf& hp, double s2cut=0.);
 
-  HProf ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt, 
+// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
+// typically = ComovRadialDistance 
+  HProf ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt, double angscale=1.,
                        double s2cut=0., int nbin=256, double kmin=0., double kmax=-1.);
 
@@ -85 +91 @@
   PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut=100., int ngen=1, const char* tit="PkNoise");
   
+  inline void SetFreqRange(double freq0=835.,double dfreq=0.5)
+  { freq0_=freq0;  dfreq_=dfreq; }
+  inline void SetAngScaleConversion(double angscale=1.)
+  { angscale_=angscale; }
   inline void SetS2Cut(double s2cut=100.)
   {  S2CUT=s2cut; }
@@ -95 +105 @@
   Four3DPk& pkn3d;
   Four2DResponse& frep;
+  double freq0_,dfreq_;
+  double angscale_; 
   double S2CUT;
   int NGEN;