Index: trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3141)
+++ trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3150)
@@ -32,5 +32,5 @@
 //-------------------------------------------------------
 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
-  : T_(T) , array_allocated_(false)
+  : T_(T) , array_allocated_(false) , lp_(0)
 {
  SetNThread();
@@ -257,8 +257,6 @@
 //    (attention on fait une fft pour realiser le spectre)
 {
- int lp=2;
-
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
- if(lp>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
 #ifdef FFTW_THREAD
  if(nthread_>0) {
@@ -270,8 +268,8 @@
  pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
  pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
- if(lp>1) PrtTim("--- ComputeFourier0: after fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after fftw_plan ---");
 
  // --- realisation d'un tableau de tirage gaussiens
- if(lp>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
  // On tient compte du pb de normalisation de FFTW3
  double sntot = sqrt((double)NRtot_);
@@ -280,13 +278,13 @@
    data_[ip] = NorRand()/sntot;
  }
- if(lp>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");
+ 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_>1) PrtTim("--- ComputeFourier0: before fft real ---");
  fftw_execute(pf_);
- if(lp>1) PrtTim("--- ComputeFourier0: after fft real ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after fft real ---");
 
  // --- On remplit avec une realisation
- if(lp>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
  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;
@@ -294,5 +292,5 @@
    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_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
    for(long j=0;j<Ny_;j++) {
      long jj = (j>Ny_/2) ? Ny_-j : j;
@@ -309,9 +307,9 @@
    }
  }
- if(lp>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---");
 
  double p = compute_power_carte();
  cout<<"Puissance dans la realisation: "<<p<<endl;
- if(lp>1) PrtTim("--- ComputeFourier0: after Computing power ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after Computing power ---");
 
 }
@@ -328,8 +326,6 @@
 //                                          sum(fb(x_i)^2) = S*N^2
 {
- int lp=2;
-
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
- if(lp>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
 #ifdef FFTW_THREAD
  if(nthread_>0) {
@@ -342,5 +338,5 @@
  pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
  //fftw_print_plan(pb_); cout<<endl;
- if(lp>1) PrtTim("--- ComputeFourier: after fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: after fftw_plan ---");
 
  // --- RaZ du tableau
@@ -348,10 +344,10 @@
 
  // --- On remplit avec une realisation
- if(lp>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
  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_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
    for(long j=0;j<Ny_;j++) {
      long jj = (j>Ny_/2) ? Ny_-j : j;
@@ -369,12 +365,12 @@
    }
  }
- if(lp>1) PrtTim("--- ComputeFourier: after Fourier realization filling ---");
+ 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_>1) PrtTim("--- ComputeFourier: after managing coefficients ---");
 
  double p = compute_power_carte();
  cout<<"Puissance dans la realisation: "<<p<<endl;
- if(lp>1) PrtTim("--- ComputeFourier: after before Computing power ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: after before Computing power ---");
 
 }
@@ -490,8 +486,7 @@
 //                          avec y = k_x*dx/2
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- FilterByPixel: before ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- FilterByPixel: before ---");
 
  for(long i=0;i<Nx_;i++) {
@@ -511,5 +506,5 @@
  }
 
- if(lp>1) PrtTim("--- FilterByPixel: after ---");
+ if(lp_>1) PrtTim("--- FilterByPixel: after ---");
 }
 
@@ -518,11 +513,10 @@
 // Calcule une realisation dans l'espace reel
 {
- int lp=2;
  check_array_alloc();
 
  // On fait la FFT
- if(lp>1) PrtTim("--- ComputeReal: before fftw backward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: before fftw backward ---");
  fftw_execute(pb_);
- if(lp>1) PrtTim("--- ComputeReal: after fftw backward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: after fftw backward ---");
 }
 
@@ -530,9 +524,8 @@
 void GeneFluct3D::ReComputeFourier(void)
 {
- int lp=2;
  check_array_alloc();
 
  // On fait la FFT
- if(lp>1) PrtTim("--- ComputeReal: before fftw forward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: before fftw forward ---");
  fftw_execute(pf_);
  // On corrige du pb de la normalisation de FFTW3
@@ -542,5 +535,5 @@
      for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
 
- if(lp>1) PrtTim("--- ComputeReal: after fftw forward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: after fftw forward ---");
 }
 
@@ -552,4 +545,5 @@
 {
  check_array_alloc();
+ if(lp_>1) PrtTim("--- ComputeSpectrum: before ---");
 
  if(herr.NBins()<0) return -1;
@@ -571,5 +565,5 @@
    }
  }
- herr.ToCorrel();
+ herr.ToVariance();
 
  // renormalize to directly compare to original spectrum
@@ -577,4 +571,6 @@
  herr *= norm;
 
+ if(lp_>1) PrtTim("--- ComputeSpectrum: after ---");
+
  return 0;
 }
@@ -583,4 +579,5 @@
 {
  check_array_alloc();
+ if(lp_>1) PrtTim("--- ComputeSpectrum2D: before ---");
 
  if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
@@ -602,5 +599,5 @@
    }
  }
- herr.ToCorrel();
+ herr.ToVariance();
 
  // renormalize to directly compare to original spectrum
@@ -608,4 +605,6 @@
  herr *= norm;
 
+ if(lp_>1) PrtTim("--- ComputeSpectrum2D: after ---");
+
  return 0;
 }
@@ -615,8 +614,7 @@
 // Recompute MASS variance in spherical top-hat (rayon=R)
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- VarianceFrReal: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- VarianceFrReal: before computation ---");
 
  long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1;
@@ -654,10 +652,10 @@
  sum /= nsum;
  sum2 = sum2/nsum - sum*sum;
- if(lp) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
+ if(lp_) 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_) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
+
+ if(lp_>1) PrtTim("--- VarianceFrReal: after computation ---");
  return nsum;
 }
@@ -726,8 +724,7 @@
 // d_rho/rho -> Mass  (add one!)
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- TurnFluct2Mass: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- TurnFluct2Mass: before computation ---");
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
@@ -736,5 +733,5 @@
  }
 
- if(lp>1) PrtTim("--- TurnFluct2Mass: after computation ---");
+ if(lp_>1) PrtTim("--- TurnFluct2Mass: after computation ---");
 }
 
@@ -742,10 +739,9 @@
 // do NOT treate negative or nul values
 {
- int lp=2;
- if(lp>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
+ if(lp_>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
 
  double m,s2;
  int_8 ngood = MeanSigma2(m,s2,0.,1e+200);
- if(lp) cout<<"TurnMass2MeanNumber: ngood="<<ngood
+ if(lp_) cout<<"TurnMass2MeanNumber: ngood="<<ngood
             <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;
 
@@ -758,5 +754,5 @@
  //   (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 "
+ if(lp_) cout<<"galaxy density move from "
             <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
  double sum = 0.;
@@ -766,7 +762,7 @@
    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_) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
+
+ if(lp_>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
  return sum;
 }
@@ -775,8 +771,7 @@
 // do NOT treate negative or nul mass  -> let it as it is
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
 
  double sum = 0.;
@@ -790,7 +785,7 @@
    }
  }
- if(lp) cout<<sum<<" galaxies put into survey"<<endl;
-
- if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
+ if(lp_) cout<<sum<<" galaxies put into survey"<<endl;
+
+ if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
  return sum;
 }
@@ -804,8 +799,7 @@
 // RETURN la masse totale
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- TurnNGal2Mass: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- TurnNGal2Mass: before computation ---");
 
  double sum = 0.;
@@ -824,7 +818,7 @@
    }
  }
- if(lp) cout<<sum<<" MSol HI mass put into survey"<<endl;
-
- if(lp>1) PrtTim("--- TurnNGal2Mass: after computation ---");
+ if(lp_) cout<<sum<<" MSol HI mass put into survey"<<endl;
+
+ if(lp_>1) PrtTim("--- TurnNGal2Mass: after computation ---");
  return sum;
 }
@@ -833,8 +827,7 @@
 // add noise to every pixels (meme les <=0 !)
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- AddNoise2Real: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- AddNoise2Real: before computation ---");
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
@@ -843,4 +836,4 @@
  }
 
- if(lp>1) PrtTim("--- AddNoise2Real: after computation ---");
-}
+ if(lp_>1) PrtTim("--- AddNoise2Real: after computation ---");
+}
