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

Timestamp:

Dec 23, 2010, 12:49:22 AM (15 years ago)

Author:

ansari

Message:

ajout config double hexagon + programme calcul pknoise a partir du resultat de repicon.cc, Reza 22/12/2010

File:

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

trunk/Cosmo/RadioBeam/specpk.cc

@@ -1 @@
-//  Classes to compute 3D power spectrum 
-// R. Ansari - Nov 2008, May 2010 
+
+/*  ------------------------ Projet BAORadio -------------------- 
+    Classes to compute 3D power spectrum and noise power spectrum
+    R. Ansari - Nov 2008 ... Dec 2010 
+---------------------------------------------------------------  */
 
 #include "specpk.h"
 #include "randr48.h"      
+#include "ctimer.h"      
 
 //------------------------------------
@@ -154 +158 @@
     }
   }
-  if (prtlev_>0)
+  if (prtlev_>2)
     cout << " Four3DPk::ComputeFourierAmp() done ..." << endl;
 }
@@ -182 +186 @@
     }
   }
-  if (prtlev_>1)  fourAmp.Show();
+  if (prtlev_>2)  fourAmp.Show();
   if (crmask) {
     POutPersist po("mask.ppf");
@@ -233 +237 @@
 void Four3DPk::ComputePkCumul(HProf& hp, double s2cut)
 {
-
+  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
@@ -251 +255 @@
         if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
         hp.Add(wk, amp2);
+        nmodeok++;
       }
     }
   }
+  if (prtlev_>1) {
+    cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size() 
+         << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
+  }
   return;
 }
 
+//-----------------------------------------------------
+// -- MassDist2D class :  2D mass distribution 
+// --- PkNoiseCalculator : Classe de calcul du spectre de bruit PNoise(k) 
+// determine par une reponse 2D de l'instrument
+//-----------------------------------------------------
+PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen, 
+                                     const char* tit)
+  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
+{
+  SetPrtLevel();
+}
+
+HProf PkNoiseCalculator::Compute()
+{
+  Timer tm(title.c_str());
+  tm.Nop();
+  HProf hnd;
+  cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
+  for(int igen=0; igen<NGEN; igen++) {
+    pkn3d.ComputeNoiseFourierAmp(frep);
+    if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
+    else pkn3d.ComputePkCumul(hnd,S2CUT);
+    if ((prtlev_>0)&&(((igen+1)%prtmodulo_)==0)) 
+      cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
+  }
+  return hnd;
+}