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

Timestamp:

Oct 26, 2011, 7:42:21 PM (14 years ago)

Author:

ansari

Message:

Corrections papiers avec les commentaires du referee (1) - Reza 26/10/2011

File:

• trunk/Cosmo/RadioBeam/specpk.cc (modified) (5 diffs)

trunk/Cosmo/RadioBeam/specpk.cc

@@ -209 +209 @@
 
 // Generate mass field Fourier Coefficient
-void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales)
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales, Vector& noisevec)
 // angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
 // typically = ComovRadialDistance 
 {
   uint_8 nmodeok=0;
-  if (angscales.Size() != fourAmp.SizeZ()) 
-    throw SzMismatchError("ComputeNoiseFourierAmp(): angscales.Size()!=fourAmp.SizeZ()");
+  if ((angscales.Size() != fourAmp.SizeZ())||(noisevec.Size() != fourAmp.SizeZ()))
+    throw SzMismatchError("ComputeNoiseFourierAmp(): angscales/noisevec.Size()!=fourAmp.SizeZ()");
   H21Conversions conv;
   // fourAmp represent 3-D fourier transform of a real input array. 
@@ -225 +225 @@
     resp.setLambda(conv.getLambda());
     double angsc=angscales(kz);
+    double noisepow=noisevec(kz);
     if (prtlev_>2)  
-      cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc << endl;
+      cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc 
+           << " NoisePow=" << noisepow << endl;
     for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
       kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
@@ -234 +236 @@
         if (rep<1.e-19)  rep=1.e-19;
         else  nmodeok++;
-        amp = 1./sqrt(rep)/sqrt(2.);
+        amp = sqrt(noisepow/rep/2.);
         fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
       }
@@ -445 +447 @@
 PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen, 
                                      const char* tit)
-  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ())
+  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ()), pnoisefac_(pk3.SizeZ())
 {
   SetFreqRange();
   SetAngScaleConversion();
+  SetPNoiseFactor();
   SetPrtLevel();
 }
@@ -461 +464 @@
        << " ... " << angscales_(angscales_.Size()-1) << endl;
   for(int igen=0; igen<NGEN; igen++) {
-    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_);
+    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_, pnoisefac_);
     if (igen==0) hnd = pkn3d.ComputePk(S2CUT,nbin,kmin,kmax,true);
     else pkn3d.ComputePkCumul();