Index: trunk/Cosmo/RadioBeam/specpk.cc
===================================================================
--- trunk/Cosmo/RadioBeam/specpk.cc	(revision 3947)
+++ trunk/Cosmo/RadioBeam/specpk.cc	(revision 4026)
@@ -6,4 +6,5 @@
 
 #include "specpk.h"
+#include "radutil.h"
 #include "randr48.h"      
 #include "ctimer.h"      
@@ -162,6 +163,9 @@
 }
 
+
 // Generate mass field Fourier Coefficient
-void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, bool crmask)
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale, bool crmask)
+// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
+// typically = ComovRadialDistance 
 {
   TMatrix<r_4> mask(fourAmp.SizeY(), fourAmp.SizeX());
@@ -176,5 +180,5 @@
       for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
 	kxx=(double)kx*dkx_;
-	rep = resp(kxx, kyy);
+	rep = resp(kxx*angscale, kyy*angscale);
 	if (crmask&&(kz==0))  mask(ky,kx)=((rep<1.e-8)?9.e9:(1./rep));
 	if (rep<1.e-8)  fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
@@ -193,4 +197,53 @@
   if (prtlev_>0)
     cout << " Four3DPk::ComputeNoiseFourierAmp() done ..." << endl;
+}
+
+// Generate mass field Fourier Coefficient
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, double angscale)
+// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
+// typically = ComovRadialDistance 
+{
+  H21Conversions conv;
+  // 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, rep, amp;
+  // sa_size_t is large integer type  
+  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
+    conv.setFrequency(f0+kz*df);
+    resp.setLambda(conv.getLambda());
+    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.);
+	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;
+  TVector< complex<TF> >  veczin(fourAmp.SizeZ()), veczout(fourAmp.SizeZ());
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+
+  for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+    for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+      //      veczin=fourAmp(Range(kx), Range(ky), Range::all());
+      for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++)  veczin(kz)=fourAmp(kx,ky,kz);
+      ffts.FFTBackward(veczin,veczout);
+      veczout /= (TF)sqrt((double)fourAmp.SizeZ());
+      //      fourAmp(Range(kx), Range(ky), Range::all())=veczout;
+      for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++)  fourAmp(kx,ky,kz)=veczout(kz);
+    }
+  }
+
+  if (prtlev_>2)  fourAmp.Show();
+  if (prtlev_>0)
+    cout << " Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df) done ..." << endl;
 }
 
@@ -270,6 +323,8 @@
 // cells with amp^2=re^2+im^2>s2cut are ignored
 // Output : noise power spectrum (profile histogram)
+// 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 s2cut, int nbin, double kmin, double kmax)
+			       double angscale, double s2cut, int nbin, double kmin, double kmax)
 {
   // The second half of the array along Y (matrix rows) contain
@@ -305,5 +360,6 @@
 	double kxx=(double)kx*dkx_;
 	double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
-	double amp2 = 1./resp(kxx, kyy); 
+	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;
@@ -349,4 +405,6 @@
   : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
 {
+  SetFreqRange();
+  SetAngScaleConversion();
   SetPrtLevel();
 }
@@ -357,7 +415,8 @@
   tm.Nop();
   HProf hnd;
-  cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
+  cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT  
+       << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscale=" << angscale_ << endl;
   for(int igen=0; igen<NGEN; igen++) {
-    pkn3d.ComputeNoiseFourierAmp(frep);
+    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscale_);
     if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
     else pkn3d.ComputePkCumul(hnd,S2CUT);
