Changeset 3155 in Sophya for trunk/Cosmo/SimLSS/genefluct3d.cc

Timestamp:

Jan 19, 2007, 6:18:38 PM (19 years ago)

Author:

cmv

Message:

PrtTim deplaces de genefluct3d.cc -> cmvobserv3d.cc cmv 19/01/2007

File:

• trunk/Cosmo/SimLSS/genefluct3d.cc (modified) (42 diffs)

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -8 +8 @@
 #include <unistd.h>
 
-#include "timing.h"
 #include "tarray.h"
 #include "pexceptions.h"
@@ -74 +73 @@
 void GeneFluct3D::setsize(long nx,long ny,long nz,double dx,double dy,double dz)
 {
+ if(lp_>1) cout<<"--- GeneFluct3D::setsize: N="<<nx<<","<<ny<<","<<nz
+                <<" D="<<dx<<","<<dy<<","<<dz<<endl;
  if(nx<=0 || dx<=0.) {
-   cout<<"GeneFluct3D::SetSize_Error: bad value(s)"<<endl;
-   throw ParmError("GeneFluct3D::SetSize_Error: bad value(s)");
+   cout<<"GeneFluct3D::setsize_Error: bad value(s)"<<endl;
+   throw ParmError("GeneFluct3D::setsize_Error: bad value(s)");
  }
 
@@ -112 +113 @@
 void GeneFluct3D::setalloc(void)
 {
+ if(lp_>1) cout<<"--- GeneFluct3D::setalloc ---"<<endl;
  // Dimensionnement du tableau complex<r_8>
  // ATTENTION: TArray adresse en memoire a l'envers du C
@@ -120 +122 @@
    array_allocated_ = true;
  } catch (...) {
-   cout<<"GeneFluct3D::SetSize_Error: Problem allocating T_"<<endl;
+   cout<<"GeneFluct3D::setalloc_Error: Problem allocating T_"<<endl;
  }
  T_.SetMemoryMapping(BaseArray::CMemoryMapping);
@@ -127 +129 @@
 void GeneFluct3D::setpointers(bool from_real)
 {
+ if(lp_>1) cout<<"--- GeneFluct3D::setpointers ---"<<endl;
  if(from_real) T_ = ArrCastR2C(R_);
    else        R_ = ArrCastC2R(T_);
@@ -147 +150 @@
 {
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
- if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: before fftw_plan ---");
+ if(lp_>1) cout<<"--- GeneFluct3D::init_fftw ---"<<endl;
 #ifdef FFTW_THREAD
  if(nthread_>0) {
-   cout<<"Computing with "<<nthread_<<" threads"<<endl;
+   cout<<"...Computing with "<<nthread_<<" threads"<<endl;
    fftw_init_threads();
    fftw_plan_with_nthreads(nthread_);
  }
 #endif
+ if(lp_>1) cout<<"...forward plan"<<endl;
  pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
+ if(lp_>1) cout<<"...backward plan"<<endl;
  pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
- if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: after fftw_plan ---");
 }
 
@@ -164 +168 @@
 void GeneFluct3D::WriteFits(string cfname,int bitpix)
 {
- cout<<"GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;
+ cout<<"--- GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;
  try {
    FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5);
@@ -191 +195 @@
 void GeneFluct3D::ReadFits(string cfname)
 {
- cout<<"GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;
+ cout<<"--- GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;
  try {
    FitsImg3DRead fimg(cfname.c_str(),0,5);
@@ -220 +224 @@
 // On ecrit soit le TArray<r_8> ou le TArray<complex <r_8> >
 {
- cout<<"GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;
+ cout<<"--- GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;
  try {
    R_.Info()["NX"] = (int_8)Nx_;
@@ -245 +249 @@
 void GeneFluct3D::ReadPPF(string cfname)
 {
- cout<<"GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;
+ cout<<"--- GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;
  try {
    bool from_real = true;
@@ -306 +310 @@
 
  // --- realisation d'un tableau de tirage gaussiens
- if(lp_>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
+ if(lp_>0) cout<<"--- ComputeFourier0: before gaussian filling ---"<<endl;
  // On tient compte du pb de normalisation de FFTW3
  double sntot = sqrt((double)NRtot_);
@@ -313 +317 @@
    data_[ip] = NorRand()/sntot;
  }
- if(lp_>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");
 
  // --- realisation d'un tableau de tirage gaussiens
- if(lp_>1) PrtTim("--- ComputeFourier0: before fft real ---");
+ if(lp_>0) cout<<"...before fft real ---"<<endl;
  fftw_execute(pf_);
- if(lp_>1) PrtTim("--- ComputeFourier0: after fft real ---");
 
  // --- On remplit avec une realisation
- if(lp_>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
+ if(lp_>0) cout<<"...before Fourier realization filling ---"<<endl;
  T_(0,0,0) = complex<r_8>(0.);  // on coupe le continue et on l'initialise
  long lmod = Nx_/10; if(lmod<1) lmod=1;
@@ -327 +329 @@
    long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_;  kx *= kx;
-   if(lp_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
    for(long j=0;j<Ny_;j++) {
      long jj = (j>Ny_/2) ? Ny_-j : j;
@@ -342 +344 @@
    }
  }
- if(lp_>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---");
-
+
+ if(lp_>0) cout<<"...computing power"<<endl;
  double p = compute_power_carte();
- cout<<"Puissance dans la realisation: "<<p<<endl;
- if(lp_>1) PrtTim("--- ComputeFourier0: after Computing power ---");
+ if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;
 
 }
@@ -365 +366 @@
 
  // --- On remplit avec une realisation
- if(lp_>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
+ if(lp_>0) cout<<"--- ComputeFourier ---"<<endl;
  long lmod = Nx_/10; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
    long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_;  kx *= kx;
-   if(lp_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
    for(long j=0;j<Ny_;j++) {
      long jj = (j>Ny_/2) ? Ny_-j : j;
@@ -386 +387 @@
    }
  }
- if(lp_>1) PrtTim("--- ComputeFourier: after Fourier realization filling ---");
 
  manage_coefficients();   // gros effet pour les spectres que l'on utilise !
- if(lp_>1) PrtTim("--- ComputeFourier: after managing coefficients ---");
-
+
+ if(lp_>0) cout<<"...computing power"<<endl;
  double p = compute_power_carte();
- cout<<"Puissance dans la realisation: "<<p<<endl;
- if(lp_>1) PrtTim("--- ComputeFourier: after before Computing power ---");
+ if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;
 
 }
@@ -401 +400 @@
 // because we want a realization of a real data in real space
 {
+ if(lp_>1) cout<<"...managing coefficients"<<endl;
  check_array_alloc();
 
@@ -421 +421 @@
    }
  }
- cout<<"Number of forced real number ="<<nreal<<endl;
+ if(lp_>1) cout<<"Number of forced real number ="<<nreal<<endl;
 
  // 2./ Les elements complexe conjugues (tous dans le plan k=0,Nyquist)
@@ -451 +451 @@
    }
  }
- cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;
+ if(lp_>1) cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;
 
  // b./ les lignes et colonnes hors continu et de nyquist
@@ -469 +469 @@
    }
  }
- cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
-
- cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
+ if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
+
+ if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
 
  return nreal+nconj1+nconj2;
@@ -507 +507 @@
 //                          avec y = k_x*dx/2
 {
- check_array_alloc();
-
- if(lp_>1) PrtTim("--- FilterByPixel: before ---");
+ if(lp_>0) cout<<"--- FilterByPixel ---"<<endl;
+ check_array_alloc();
 
  for(long i=0;i<Nx_;i++) {
@@ -527 +526 @@
  }
 
- if(lp_>1) PrtTim("--- FilterByPixel: after ---");
 }
 
@@ -534 +532 @@
 // Calcule une realisation dans l'espace reel
 {
+ if(lp_>0) cout<<"--- ComputeReal ---"<<endl;
  check_array_alloc();
 
  // On fait la FFT
- if(lp_>1) PrtTim("--- ComputeReal: before fftw backward ---");
  fftw_execute(pb_);
- if(lp_>1) PrtTim("--- ComputeReal: after fftw backward ---");
 }
 
@@ -545 +542 @@
 void GeneFluct3D::ReComputeFourier(void)
 {
+ if(lp_>0) cout<<"--- ReComputeFourier ---"<<endl;
  check_array_alloc();
 
  // On fait la FFT
- if(lp_>1) PrtTim("--- ComputeReal: before fftw forward ---");
  fftw_execute(pf_);
  // On corrige du pb de la normalisation de FFTW3
@@ -556 +553 @@
      for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
 
- if(lp_>1) PrtTim("--- ComputeReal: after fftw forward ---");
 }
 
@@ -565 +561 @@
 // cad T(kz,ky,kx) avec  0<kz<kz_nyq  -ky_nyq<ky<ky_nyq  -kx_nyq<kx<kx_nyq
 {
- check_array_alloc();
- if(lp_>1) PrtTim("--- ComputeSpectrum: before ---");
+ if(lp_>0) cout<<"--- ComputeSpectrum ---"<<endl;
+ check_array_alloc();
 
  if(herr.NBins()<0) return -1;
@@ -592 +588 @@
  herr *= norm;
 
- if(lp_>1) PrtTim("--- ComputeSpectrum: after ---");
-
  return 0;
 }
@@ -599 +593 @@
 int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr)
 {
- check_array_alloc();
- if(lp_>1) PrtTim("--- ComputeSpectrum2D: before ---");
+ if(lp_>0) cout<<"--- ComputeSpectrum2D ---"<<endl;
+ check_array_alloc();
 
  if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
@@ -626 +620 @@
  herr *= norm;
 
- if(lp_>1) PrtTim("--- ComputeSpectrum2D: after ---");
-
  return 0;
 }
@@ -635 +627 @@
 // Recompute MASS variance in spherical top-hat (rayon=R)
 {
- check_array_alloc();
-
- if(lp_>1) PrtTim("--- VarianceFrReal: before computation ---");
+ if(lp_>0) cout<<"--- VarianceFrReal ---"<<endl;
+ check_array_alloc();
 
  long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1;
  long dny = long(R/Dy_+0.5); if(dny<=0) dny = 1;
  long dnz = long(R/Dz_+0.5); if(dnz<=0) dnz = 1;
- cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
+ if(lp_>0) cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
 
  double sum=0., sum2=0., r2 = R*R; int_8 nsum=0;
@@ -673 +664 @@
  sum /= nsum;
  sum2 = sum2/nsum - sum*sum;
- if(lp_) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
+ if(lp_>0) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
 
  var = sum2/(sum*sum);  // <dM>^2/<M>^2
- if(lp_) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
-
- if(lp_>1) PrtTim("--- VarianceFrReal: after computation ---");
+ if(lp_>0) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
+
  return nsum;
 }
@@ -745 +735 @@
 // d_rho/rho -> Mass  (add one!)
 {
- check_array_alloc();
-
- if(lp_>1) PrtTim("--- TurnFluct2Mass: before computation ---");
+ if(lp_>0) cout<<"--- TurnFluct2Mass ---"<<endl;
+ check_array_alloc();
+
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
@@ -753 +743 @@
    data_[ip] += 1.;
  }
-
- if(lp_>1) PrtTim("--- TurnFluct2Mass: after computation ---");
 }
 
@@ -760 +748 @@
 // do NOT treate negative or nul values
 {
- if(lp_>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
+ if(lp_>0) cout<<"--- TurnMass2MeanNumber ---"<<endl;
 
  double m,s2;
  int_8 ngood = MeanSigma2(m,s2,0.,1e+200);
- if(lp_) cout<<"TurnMass2MeanNumber: ngood="<<ngood
-            <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;
+ if(lp_>0) cout<<"...ngood="<<ngood
+               <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;
 
  // On doit mettre n*Vol galaxies dans notre survey
@@ -775 +763 @@
  //   (rappel sur tout les pixels <1+d_rho/rho>=1)
  double dn = n_by_mpc3*Vol_/m /(double)ngood;  // nb de gal a mettre ds 1 px
- if(lp_) cout<<"galaxy density move from "
-            <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
+ if(lp_>0) cout<<"...galaxy density move from "
+               <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
  double sum = 0.;
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
@@ -783 +771 @@
    if(data_[ip]>0.) sum += data_[ip];  // mais on ne les compte pas
  }
- if(lp_) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
-
- if(lp_>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
+ if(lp_>0) cout<<sum<<"...galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
+
  return sum;
 }
@@ -792 +779 @@
 // do NOT treate negative or nul mass  -> let it as it is
 {
- check_array_alloc();
-
- if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
+ if(lp_>0) cout<<"--- ApplyPoisson ---"<<endl;
+ check_array_alloc();
 
  double sum = 0.;
@@ -806 +792 @@
    }
  }
- if(lp_) cout<<sum<<" galaxies put into survey"<<endl;
-
- if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
+ if(lp_>0) cout<<sum<<" galaxies put into survey"<<endl;
+
  return sum;
 }
@@ -820 +805 @@
 // RETURN la masse totale
 {
- check_array_alloc();
-
- if(lp_>1) PrtTim("--- TurnNGal2Mass: before computation ---");
+ if(lp_>0) cout<<"--- TurnNGal2Mass ---"<<endl;
+ check_array_alloc();
 
  double sum = 0.;
@@ -839 +823 @@
    }
  }
- if(lp_) cout<<sum<<" MSol HI mass put into survey"<<endl;
-
- if(lp_>1) PrtTim("--- TurnNGal2Mass: after computation ---");
+ if(lp_>0) cout<<sum<<" MSol HI mass put into survey"<<endl;
+
  return sum;
 }
@@ -848 +831 @@
 // add noise to every pixels (meme les <=0 !)
 {
- check_array_alloc();
-
- if(lp_>1) PrtTim("--- AddNoise2Real: before computation ---");
+ if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<endl;
+ check_array_alloc();
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
@@ -856 +838 @@
    data_[ip] += snoise*NorRand();
  }
-
- if(lp_>1) PrtTim("--- AddNoise2Real: after computation ---");
-}
+}