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

Timestamp:

Oct 17, 2011, 10:56:41 AM (14 years ago)

Author:

ansari

Message:

Implementatiom prise en compte dA(z) ds la calcul de bruit Pnoise(k) - Reza 17/10/2011

File:

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

trunk/Cosmo/RadioBeam/specpk.cc

@@ -127 +127 @@
   SetPrtLevel();
   SetCellSize();
+  hp_pk_p_=NULL;  hmcnt_p_=NULL;  hmcntok_p_=NULL;   s2cut_=0.;
 }
 // Constructor
@@ -134 +135 @@
   SetPrtLevel();
   SetCellSize();
-}
-
+  hp_pk_p_=NULL;  hmcnt_p_=NULL;  hmcntok_p_=NULL;   s2cut_=0.;
+}
+
+// Destructor
+Four3DPk::~Four3DPk()
+{
+  if (hp_pk_p_) delete hp_pk_p_;
+  if (hmcnt_p_) delete hmcnt_p_;
+  if (hmcntok_p_) delete hmcntok_p_;
+}
 
 // Generate mass field Fourier Coefficient
@@ -200 +209 @@
 
 // Generate mass field Fourier Coefficient
-void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, double angscale)
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales)
 // angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
 // typically = ComovRadialDistance 
 {
+  if (angscales.Size() != fourAmp.SizeZ()) 
+    throw SzMismatchError("ComputeNoiseFourierAmp(): angscales.Size()!=fourAmp.SizeZ()");
   H21Conversions conv;
   // fourAmp represent 3-D fourier transform of a real input array. 
@@ -212 +223 @@
     conv.setFrequency(f0+kz*df);
     resp.setLambda(conv.getLambda());
+    double angsc=angscales(kz);
+    if (prtlev_>2)  
+      cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc << endl;
     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.);
+        rep = resp(kxx*angsc, kyy*angsc);
+        if (rep<1.e-19)  fourAmp(kx, ky, kz) = complex<TF>(9.e19,0.);
         else {
           amp = 1./sqrt(rep)/sqrt(2.);
@@ -264 +278 @@
 // cells with amp^2=re^2+im^2>s2cut are ignored
 // Output : power spectrum (profile histogram)
-HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax)
+HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax, bool fgmodcnt)
 {
   // The second half of the array along Y (matrix rows) contain
@@ -279 +293 @@
     kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
   }
-  if (nbin<2) nbin=128; 
-  HProf hp(kmin, kmax, nbin);
-  hp.SetErrOpt(false);
-  ComputePkCumul(hp, s2cut);
-  return hp;
+  if (nbin<2) nbin=256; 
+  hp_pk_p_ = new HProf(kmin, kmax, nbin);
+  hp_pk_p_->SetErrOpt(false);
+  if (fgmodcnt) {
+    hmcnt_p_ = new Histo(kmin, kmax, nbin);
+    hmcntok_p_ = new Histo(kmin, kmax, nbin);
+  }
+  s2cut_=s2cut; 
+  ComputePkCumul();
+  return *hp_pk_p_;
 }
 
 // Compute power spectrum as a function of wave number k 
 // Cumul dans hp - cells with amp^2=re^2+im^2>s2cut are ignored
-void Four3DPk::ComputePkCumul(HProf& hp, double s2cut)
+void Four3DPk::ComputePkCumul()
 {
   uint_8 nmodeok=0;
@@ -303 +322 @@
         double kxx=(double)kx*dkx_;
         complex<TF> za = fourAmp(kx, ky, kz);
-        if (za.real()>8.e9) continue;
+        //      if (za.real()>8.e19) continue;
         double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
         double amp2 = za.real()*za.real()+za.imag()*za.imag();
-        if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
-        hp.Add(wk, amp2);
+        if (hmcnt_p_) hmcnt_p_->Add(wk);
+        if ((s2cut_>1.e-9)&&(amp2>s2cut_))  continue;
+        if (hmcntok_p_) hmcntok_p_->Add(wk);
+        hp_pk_p_->Add(wk, amp2);
         nmodeok++;
       }
     }
   }
-  if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
+  if ((prtlev_>1)||((prtlev_>0)&&(s2cut_>1.e-9))) {
     cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size() 
          << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
@@ -325 +346 @@
 // 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 angscale, double s2cut, int nbin, double kmin, double kmax)
+HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, double angscale, double s2cut, int nbin, double kmin, double kmax)
 {
   // The second half of the array along Y (matrix rows) contain
@@ -341 +361 @@
     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);
+  if (nbin<2) nbin=256; 
+  hp_pk_p_ = new HProf(kmin, kmax, nbin);
+  hp_pk_p_->SetErrOpt(false);
+  hmcnt_p_ = new Histo(kmin, kmax, nbin);
+  hmcntok_p_ = new Histo(kmin, kmax, nbin);
+  s2cut_=s2cut; 
 
   uint_8 nmodeok=0;
@@ -362 +383 @@
         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;
-        hmcntok.Add(wk);
-        hp.Add(wk, amp2);
+        hmcnt_p_->Add(wk);
+        if ((s2cut_>1.e-9)&&(amp2>s2cut_))  continue;
+        hmcntok_p_->Add(wk);
+        hp_pk_p_->Add(wk, amp2);
         nmodeok++;
       }
@@ -375 +396 @@
   }
 
-  fracmodok=hmcntok/hmcnt;
-
+  return *hp_pk_p_;
+}
+
+// Fills a data table from the computed P(k) profile histogram and mode count 
+Histo Four3DPk::FillPkDataTable(DataTable& dt)
+{
+  if (hp_pk_p_==NULL) throw ParmError("Four3DPk::FillPkDataTable P(k) NOT computed");
+  if ((hmcnt_p_==NULL)||(hmcntok_p_==NULL)) throw ParmError("Four3DPk::FillPkDataTable Mode count histos NOT filled");
+  HProf& hp=(*hp_pk_p_);
+  Histo& hmcnt=(*hmcnt_p_);
+  Histo& hmcntok=(*hmcntok_p_);
+  Histo fracmodok=hmcntok/hmcnt;
   char* nomcol[5] = {"k","pnoise","nmode","nmodok","fracmodok"};
   dt.Clear();
@@ -393 +424 @@
     dt.AddRow(dtr);
   }
-  return hp;
+  return fracmodok;
 }
 
@@ -403 +434 @@
 PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen, 
                                      const char* tit)
-  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
+  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ())
 {
   SetFreqRange();
@@ -410 +441 @@
 }
 
-HProf PkNoiseCalculator::Compute()
+HProf PkNoiseCalculator::Compute(int nbin, double kmin, double kmax)
 {
   Timer tm(title.c_str());
@@ -416 +447 @@
   HProf hnd;
   cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT  
-       << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscale=" << angscale_ << endl;
+       << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscales=" << angscales_(0) 
+       << " ... " << angscales_(angscales_.Size()-1) << endl;
   for(int igen=0; igen<NGEN; igen++) {
-    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscale_);
-    if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
-    else pkn3d.ComputePkCumul(hnd,S2CUT);
+    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_);
+    if (igen==0) hnd = pkn3d.ComputePk(S2CUT,nbin,kmin,kmax,true);
+    else pkn3d.ComputePkCumul();
     if ((prtlev_>0)&&(igen>0)&&(((igen-1)%prtmodulo_)==0)) 
       cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
   }
-  return hnd;
+  return pkn3d.GetPk() ;
 }