Changeset 4028 in Sophya for trunk

Timestamp:

Oct 18, 2011, 10:56:28 AM (14 years ago)

Author:

ansari

Message:

Debug de calcul Pnoise3D(k) avec D_A(z), Reza 18/10/2011

Location:

trunk/Cosmo/RadioBeam

Files:

• pknoise.cc (modified) (4 diffs)
• plpkn.pic (modified) (1 diff)
• specpk.cc (modified) (8 diffs)

trunk/Cosmo/RadioBeam/pknoise.cc

@@ -90 +90 @@
   double H_z=120.5;  // Hubble_param(z) 
   int box3dsz[3]={512,512,256};
-  double cellsz[3]={3.,3.,1.5};
+  double cellsz[3]={3.,3.,3.};
   int NBINPK=256;
   int prtlev=0;
@@ -159 +159 @@
     redshift2da.ReadXYFromFile(interpfilename);
     SLinInterp1D redshift2hz;
-    interpfilename="redshift_hz.txt";
+    interpfilename="redshift_Hz.txt";
     redshift2hz.ReadXYFromFile(interpfilename);
     comov_dA_z=redshift2da(z_Redshift);
@@ -211 +211 @@
     Vector angscales(box3dsz[2]);
     angscales=angscale;
-    if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[3]=128; }
+    if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[2]=128; }
     else {
       angscale=comov_dA_z;
@@ -225 +225 @@
          << box3dsz[2] << ")" << endl;
     Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
-    cout << " pknoise[2.d]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc << endl;
+    cout << " pknoise[2.d]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc 
+         << ") cell size (Mpc) : " << cellsz[0] << "x" << cellsz[1] << "x" << cellsz[2] << endl;
     m3d.SetCellSize(dkxmpc, dkympc, dkzmpc);
     m3d.SetPrtLevel(prtlev,prtmod);

trunk/Cosmo/RadioBeam/plpkn.pic

@@ -8 +8 @@
 delobjs *
 # openppf ../../PkNoise/cmvhpkz.ppf 
-# openppf ../pkz0p25.ppf
-openppf ../pkz0p7.ppf
+# openppf pkz0p25.ppf
+openppf pkz0p7.ppf
 rename hpkz hpkz0p7
-openppf ../pkz1p0.ppf
+openppf pkz1p0.ppf
 rename hpkz hpkz1p0
 

trunk/Cosmo/RadioBeam/specpk.cc

@@ -213 +213 @@
 // typically = ComovRadialDistance 
 {
+  uint_8 nmodeok=0;
   if (angscales.Size() != fourAmp.SizeZ()) 
     throw SzMismatchError("ComputeNoiseFourierAmp(): angscales.Size()!=fourAmp.SizeZ()");
@@ -231 +232 @@
         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 +260 @@
   }
 
-  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 +318 @@
   // 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 +357 @@
   // 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 +364 @@
     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 +379 @@
   // 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 +393 @@
         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;
   }