Changeset 3947 in Sophya for trunk/Cosmo/RadioBeam/specpk.cc

Timestamp:

Feb 14, 2011, 12:58:29 AM (15 years ago)

Author:

ansari

Message:

Amelioration et modifs diverses lors de la redaction du papier, Reza 13/02/2011

File:

• trunk/Cosmo/RadioBeam/specpk.cc (modified) (1 diff)

trunk/Cosmo/RadioBeam/specpk.cc

@@ -266 +266 @@
 }
 
+// Compute noise power spectrum as a function of wave number k 
+// No random generation, computes profile of noise sigma^2
+// cells with amp^2=re^2+im^2>s2cut are ignored
+// Output : noise power spectrum (profile histogram)
+HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt, 
+                               double s2cut, int nbin, double kmin, double kmax)
+{
+  // The second half of the array along Y (matrix rows) contain
+  // negative frequencies
+  //  int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
+  // The profile histogram will contain the mean value of noise sigma^2
+  // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
+  //  if (nbin < 1) nbin = nbh/2;
+  if ((kmax<0.)||(kmax<kmin)) {
+    kmin=0.;
+    double maxx=fourAmp.SizeX()*dkx_;
+    double maxy=fourAmp.SizeY()*dky_/2;
+    double maxz=fourAmp.SizeZ()*dkz_/2;
+    kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
+  }
+  if (nbin<2) nbin=128; 
+  HProf hp(kmin, kmax, nbin);
+  hp.SetErrOpt(false);
+  Histo hmcnt(kmin, kmax, nbin);
+  Histo hmcntok(kmin, kmax, nbin);
+
+  uint_8 nmodeok=0;
+  // 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;
+  // sa_size_t is large integer type  
+  // We ignore 0th term in all frequency directions ...
+  for(sa_size_t kz=1; kz<fourAmp.SizeZ(); kz++) {
+    kzz =  (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_; 
+    for(sa_size_t ky=1; ky<fourAmp.SizeY(); ky++) {
+      kyy =  (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
+      for(sa_size_t kx=1; kx<fourAmp.SizeX(); kx++) {  // ignore the 0th coefficient (constant term)
+        double kxx=(double)kx*dkx_;
+        double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
+        double amp2 = 1./resp(kxx, kyy); 
+        hmcnt.Add(wk);
+        if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
+        hmcntok.Add(wk);
+        hp.Add(wk, amp2);
+        nmodeok++;
+      }
+    }
+  }
+  if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
+    cout << " Four3DPk::ComputeNoisePk/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size() 
+         << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
+  }
+
+  fracmodok=hmcntok/hmcnt;
+
+  char* nomcol[5] = {"k","pnoise","nmode","nmodok","fracmodok"};
+  dt.Clear();
+  dt.AddDoubleColumn(nomcol[0]);    
+  dt.AddDoubleColumn(nomcol[1]);    
+  dt.AddIntegerColumn(nomcol[2]);    
+  dt.AddIntegerColumn(nomcol[3]);    
+  dt.AddFloatColumn(nomcol[4]);    
+  DataTableRow  dtr = dt.EmptyRow();
+  for(int_4 ib=0; ib<hp.NBins(); ib++) {
+    dtr[0]=hp.BinCenter(ib);
+    dtr[1]=hp(ib);
+    dtr[2]=hmcnt(ib);
+    dtr[3]=hmcntok(ib);
+    dtr[4]=fracmodok(ib);
+    dt.AddRow(dtr);
+  }
+  return hp;
+}
+
 //-----------------------------------------------------
 // -- MassDist2D class :  2D mass distribution