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

Timestamp:

Jan 18, 2007, 7:23:56 PM (19 years ago)

Author:

cmv

Message:

suite modifs structure code + bug HistoErr/2D dans ComputeSpectrum cmv 18/01/2007

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -32 +32 @@
 //-------------------------------------------------------
 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
-  : T_(T) , array_allocated_(false)
+  : T_(T) , array_allocated_(false) , lp_(0)
 {
  SetNThread();
@@ -257 +257 @@
 //    (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 +268 @@
  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 +278 @@
    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 +292 @@
    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 +307 @@
    }
  }
- 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 +326 @@
 //                                          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 +338 @@
  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 +344 @@
 
  // --- 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 +365 @@
    }
  }
- 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 +486 @@
 //                          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 +506 @@
  }
 
- if(lp>1) PrtTim("--- FilterByPixel: after ---");
+ if(lp_>1) PrtTim("--- FilterByPixel: after ---");
 }
 
@@ -518 +513 @@
 // 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 +524 @@
 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 +535 @@
      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 +545 @@
 {
  check_array_alloc();
+ if(lp_>1) PrtTim("--- ComputeSpectrum: before ---");
 
  if(herr.NBins()<0) return -1;
@@ -571 +565 @@
    }
  }
- herr.ToCorrel();
+ herr.ToVariance();
 
  // renormalize to directly compare to original spectrum
@@ -577 +571 @@
  herr *= norm;
 
+ if(lp_>1) PrtTim("--- ComputeSpectrum: after ---");
+
  return 0;
 }
@@ -583 +579 @@
 {
  check_array_alloc();
+ if(lp_>1) PrtTim("--- ComputeSpectrum2D: before ---");
 
  if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
@@ -602 +599 @@
    }
  }
- herr.ToCorrel();
+ herr.ToVariance();
 
  // renormalize to directly compare to original spectrum
@@ -608 +605 @@
  herr *= norm;
 
+ if(lp_>1) PrtTim("--- ComputeSpectrum2D: after ---");
+
  return 0;
 }
@@ -615 +614 @@
 // 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 +652 @@
  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 +724 @@
 // 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 +733 @@
  }
 
- if(lp>1) PrtTim("--- TurnFluct2Mass: after computation ---");
+ if(lp_>1) PrtTim("--- TurnFluct2Mass: after computation ---");
 }
 
@@ -742 +739 @@
 // 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 +754 @@
  //   (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 +762 @@
    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 +771 @@
 // 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 +785 @@
    }
  }
- 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 +799 @@
 // 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 +818 @@
    }
  }
- 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 +827 @@
 // 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 +836 @@
  }
 
- if(lp>1) PrtTim("--- AddNoise2Real: after computation ---");
-}
+ if(lp_>1) PrtTim("--- AddNoise2Real: after computation ---");
+}