Changeset 4026 in Sophya for trunk/Cosmo/RadioBeam/specpk.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/specpk.cc (modified) (8 diffs)

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