Index: trunk/Cosmo/RadioBeam/specpk.cc
===================================================================
--- trunk/Cosmo/RadioBeam/specpk.cc	(revision 4027)
+++ trunk/Cosmo/RadioBeam/specpk.cc	(revision 4028)
@@ -213,4 +213,5 @@
 // typically = ComovRadialDistance 
 {
+  uint_8 nmodeok=0;
   if (angscales.Size() != fourAmp.SizeZ()) 
     throw SzMismatchError("ComputeNoiseFourierAmp(): angscales.Size()!=fourAmp.SizeZ()");
@@ -231,15 +232,17 @@
 	kxx=(double)kx*dkx_;
 	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.);
-	  fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
-	}
-      }
-    }
-  }
-
-  if (prtlev_>1)
-    cout << " Four3DPk::ComputeNoiseFourierAmp(...) Computing FFT along frequency ..." << endl;
+	if (rep<1.e-19)  rep=1.e-19;
+	else  nmodeok++;
+	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 ... \n" 
+	 << " NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size() 
+	 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
+  }
   TVector< complex<TF> >  veczin(fourAmp.SizeZ()), veczout(fourAmp.SizeZ());
   FFTWServer ffts(true);                     
@@ -257,5 +260,5 @@
   }
 
-  if (prtlev_>2)  fourAmp.Show();
+  //  if (prtlev_>2)  fourAmp.Show();
   if (prtlev_>0)
     cout << " Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df) done ..." << endl;
@@ -315,10 +318,10 @@
   // 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++) {
+  for(sa_size_t kz=0; 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++) {
+    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=1; kx<fourAmp.SizeX(); kx++) {  // ignore the 0th coefficient (constant term)
-	double kxx=(double)kx*dkx_;
+      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {  // ignore the 0th coefficient (constant term)
+	kxx=(double)kx*dkx_;
 	complex<TF> za = fourAmp(kx, ky, kz);
 	//	if (za.real()>8.e19) continue;
@@ -354,4 +357,5 @@
   // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
   //  if (nbin < 1) nbin = nbh/2;
+  double kmax2d=0.;
   if ((kmax<0.)||(kmax<kmin)) {
     kmin=0.;
@@ -360,4 +364,5 @@
     double maxz=fourAmp.SizeZ()*dkz_/2;
     kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
+    kmax2d=sqrt(maxx*maxx+maxy*maxy);
   }
   if (nbin<2) nbin=256; 
@@ -374,10 +379,11 @@
   // 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++) {
+  for(sa_size_t kz=0; 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++) {
+    uint_8 nmodeok_kz=0;
+    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=1; kx<fourAmp.SizeX(); kx++) {  // ignore the 0th coefficient (constant term)
-	double kxx=(double)kx*dkx_;
+      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {  // ignore the 0th coefficient (constant term)
+	kxx=(double)kx*dkx_;
 	double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
 	double rep=resp(kxx*angscale, kyy*angscale);
@@ -387,10 +393,15 @@
 	hmcntok_p_->Add(wk);
 	hp_pk_p_->Add(wk, amp2);
-	nmodeok++;
-      }
-    }
-  }
+	nmodeok++;    	nmodeok_kz++;
+      }
+    }
+    if (prtlev_>2)
+	cout << " Four3DPk::ComputeNoisePk(kz="<<kz<<") ModeOK_Kz=" << nmodeok_kz 
+	     << " FracOK=" << 100.*(double)nmodeok_kz/(double)fourAmp.SizeX()/(double)fourAmp.SizeY() << endl;
+   }
   if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
-    cout << " Four3DPk::ComputeNoisePk/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size() 
+    cout << " Four3DPk::ComputeNoisePk/Info : angscale=" << angscale 
+	 << " kmin,kmax=" << kmin << "," << kmax << " kmax2d_ang=" << kmax2d*angscale 
+	 << " \n ... NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size() 
 	 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
   }
