source: Sophya/trunk/Cosmo/SimLSS/genefluct3d.cc@ 3115

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

#include "timing.h"
#include "tarray.h"
#include "pexceptions.h"
#include "perandom.h"
#include "srandgen.h"

#include "constcosmo.h"
#include "integfunc.h"
#include "geneutils.h"

#include "genefluct3d.h"

//#define FFTW_THREAD

#define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))

//-------------------------------------------------------
GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz)
  : T_(T) , pkz_(pkz)
{
 SetNThread();
 SetSize(1,-1,1,1.,-1.,1.);
}

GeneFluct3D::~GeneFluct3D(void)
{
 fftw_destroy_plan(pf_);
 fftw_destroy_plan(pb_);
#ifdef FFTW_THREAD
 if(nthread_>0) fftw_cleanup_threads();
#endif
}

//-------------------------------------------------------
void GeneFluct3D::SetSize(int nx,int ny,int nz,double dx,double dy,double dz)
{
 if(nx<=0 || dx<=0. ) {
   cout<<"GeneFluct3D::SetSize_Error: bad value(s)"<<endl;
   throw ParmError("GeneFluct3D::SetSize_Error: bad value(s)");
 }

 Tcontent_ = 0;

 // Les taille des tableaux
 Nx_ = nx;
 Ny_ = ny;  if(Ny_ <= 0) Ny_ = Nx_;
 Nz_ = nz;  if(Nz_ <= 0) Nz_ = Nx_;
 NRtot_ = Nx_*Ny_*Nz_; // nombre de pixels dans le survey
 NCz_ =  Nz_/2 +1;
 NTz_ = 2*NCz_;
 SzK_[2] = Nx_; SzK_[1] = Ny_; SzK_[0] = NCz_; // a l'envers

 // le pas dans l'espace (Mpc)
 Dx_ = dx;
 Dy_ = dy; if(Dy_ <= 0.) Dy_ = Dx_;
 Dz_ = dz; if(Dz_ <= 0.) Dz_ = Dx_;
 dVol_ = Dx_*Dy_*Dz_;
 Vol_ = (Nx_*Dx_)*(Ny_*Dy_)*(Nz_*Dz_);

 // Le pas dans l'espace de Fourier (Mpc^-1)
 Dkx_ = 2.*M_PI/(Nx_*Dx_);
 Dky_ = 2.*M_PI/(Ny_*Dy_);
 Dkz_ = 2.*M_PI/(Nz_*Dz_);
 Dk3_ = Dkx_*Dky_*Dkz_;
 
 // La frequence de Nyquist en k (Mpc^-1)
 Knyqx_ = M_PI/Dx_;
 Knyqy_ = M_PI/Dy_;
 Knyqz_ = M_PI/Dz_;

}

//-------------------------------------------------------
void GeneFluct3D::Print(void)
{
 cout<<"GeneFluct3D(T_typ="<<Tcontent_<<"): z="<<pkz_.GetZ()<<endl;
 cout<<"Space Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<Nz_<<" ("<<NTz_<<")  size="
     <<NRtot_<<endl;
 cout<<"      Resol: dx="<<Dx_<<" dy="<<Dy_<<" dz="<<Dz_<<" Mpc"
     <<", dVol="<<dVol_<<", Vol="<<Vol_<<" Mpc^3"<<endl;
 cout<<"Fourier Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<NCz_<<endl;
 cout<<"        Resol: dkx="<<Dkx_<<" dky="<<Dky_<<" dkz="<<Dkz_<<" Mpc^-1"
     <<", Dk3="<<Dk3_<<" Mpc^-3"<<endl;
 cout<<"          (2Pi/k: "<<2.*M_PI/Dkx_<<" "<<2.*M_PI/Dky_<<" "<<2.*M_PI/Dkz_<<" Mpc)"<<endl;
 cout<<"      Nyquist: kx="<<Knyqx_<<" ky="<<Knyqy_<<" kz="<<Knyqz_<<" Mpc^-1"
     <<", Kmax="<<GetKmax()<<" Mpc^-1"<<endl;
 cout<<"          (2Pi/k: "<<2.*M_PI/Knyqx_<<" "<<2.*M_PI/Knyqy_<<" "<<2.*M_PI/Knyqz_<<" Mpc)"<<endl;
}

//-------------------------------------------------------
void GeneFluct3D::ComputeFourier0(void)
// cf ComputeFourier() mais avec autre methode de realisation du spectre
//    (attention on fait une fft pour realiser le spectre)
{
 int lp=2;

 T_.ReSize(3,SzK_);
 T_.SetMemoryMapping(BaseArray::CMemoryMapping);

 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
 if(lp>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
 fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
 double *data = (double *) (&T_(0,0,0));
#ifdef FFTW_THREAD
 if(nthread_>0) {
   cout<<"Computing with "<<nthread_<<" threads"<<endl;
   fftw_init_threads();
   fftw_plan_with_nthreads(nthread_);
 }
#endif
 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 ---");

 // --- realisation d'un tableau de tirage gaussiens
 if(lp>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
 // On tient compte du pb de normalisation de FFTW3
 double sntot = sqrt((double)NRtot_);
 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   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 ---");
 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 ---");
 T_(0,0,0) = complex<r_8>(0.);  // on coupe le continue et on l'initialise
 int lmod = Nx_/10; if(lmod<1) lmod=1;
 for(int i=0;i<Nx_;i++) {
   int 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;
   for(int j=0;j<Ny_;j++) {
     int jj = (j>Ny_/2) ? Ny_-j : j;
     double ky = jj*Dky_;  ky *= ky;
     for(int l=0;l<NCz_;l++) {
       double kz = l*Dkz_;  kz *= kz;
       if(i==0 && j==0 && l==0) continue; // Suppression du continu
       double k = sqrt(kx+ky+kz);
       // cf normalisation: Peacock, Cosmology, formule 16.38 p504
       double pk = pkz_(k)/Vol_;
       // ici pas de "/2" a cause de la remarque ci-dessus
       T_(l,j,i) *= sqrt(pk);
     }
   }
 }
 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 ---");

 Tcontent_ = 1;

}

//-------------------------------------------------------
void GeneFluct3D::ComputeFourier(void)
// Calcule une realisation du spectre "pkz_"
// Attention: dans TArray le premier indice varie le + vite
// Explication normalisation: see Coles & Lucchin, Cosmology, p264-265
// FFTW3: on note N=Nx*Ny*Nz
// f  --(FFT)-->  F = TF(f)  --(FFT^-1)-->  fb = TF^-1(F) = TF^-1(TF(f))
// sum(f(x_i)^2) = S
//                sum(F(nu_i)^2) = S*N
//                                          sum(fb(x_i)^2) = S*N^2
{
 int lp=2;
 
 // --- Dimensionnement du tableau
 // ATTENTION: TArray adresse en memoire a l'envers du C
 //            Tarray(n1,n2,n3) == Carray[n3][n2][n1]
 T_.ReSize(3,SzK_);
 T_.SetMemoryMapping(BaseArray::CMemoryMapping);

 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
 if(lp>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
 fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
 double *data = (double *) (&T_(0,0,0));
#ifdef FFTW_THREAD
 if(nthread_>0) {
   cout<<"Computing with "<<nthread_<<" threads"<<endl;
   fftw_init_threads();
   fftw_plan_with_nthreads(nthread_);
 }
#endif
 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);
 //fftw_print_plan(pb_); cout<<endl;
 if(lp>1) PrtTim("--- ComputeFourier: after fftw_plan ---");

 // --- RaZ du tableau
 T_ = complex<r_8>(0.);

 // --- On remplit avec une realisation
 if(lp>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
 int lmod = Nx_/10; if(lmod<1) lmod=1;
 for(int i=0;i<Nx_;i++) {
   int 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;
   for(int j=0;j<Ny_;j++) {
     int jj = (j>Ny_/2) ? Ny_-j : j;
     double ky = jj*Dky_;  ky *= ky;
     for(int l=0;l<NCz_;l++) {
       double kz = l*Dkz_;  kz *= kz;
       if(i==0 && j==0 && l==0) continue; // Suppression du continu
       double k = sqrt(kx+ky+kz);
       // cf normalisation: Peacock, Cosmology, formule 16.38 p504
       double pk = pkz_(k)/Vol_;
       // Explication de la division par 2: voir perandom.cc
       // ou egalement Coles & Lucchin, Cosmology formula 13.7.2 p279
       T_(l,j,i) = ComplexGaussRan(sqrt(pk/2.));
     }
   }
 }
 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 ---");

 double p = compute_power_carte();
 cout<<"Puissance dans la realisation: "<<p<<endl;
 if(lp>1) PrtTim("--- ComputeFourier: after before Computing power ---");

 Tcontent_ = 1;

}

int GeneFluct3D::manage_coefficients(void)
// Take into account the real and complexe conjugate coefficients
// because we want a realization of a real data in real space
{
 fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));

 // 1./ Le Continu et Nyquist sont reels
 int nreal = 0;
 for(int kk=0;kk<2;kk++) {
   int k=0;  // continu
   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
   for(int jj=0;jj<2;jj++) {
     int j=0;
     if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
     for(int ii=0;ii<2;ii++) {
       int i=0;
       if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
       size_t ip = k+NCz_*(j+Ny_*i);
       //cout<<"i="<<i<<" j="<<j<<" k="<<k<<" = ("<<fdata[ip][0]<<","<<fdata[ip][1]<<")"<<endl;
       fdata[ip][1] = 0.; fdata[ip][0] *= M_SQRT2;
       nreal++;
     }
   }
 }
 cout<<"Number of forced real number ="<<nreal<<endl;

 // 2./ Les elements complexe conjugues (tous dans le plan k=0,Nyquist)

 // a./ les lignes et colonnes du continu et de nyquist
 int nconj1 = 0;
 for(int kk=0;kk<2;kk++) {
   int k=0;  // continu
   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
   for(int jj=0;jj<2;jj++) { // selon j
     int j=0;
     if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
     for(int i=1;i<(Nx_+1)/2;i++) {
       size_t ip = k+NCz_*(j+Ny_*i);
       size_t ip1 = k+NCz_*(j+Ny_*(Nx_-i));
       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
       nconj1++;
     }
   }
   for(int ii=0;ii<2;ii++) {
     int i=0;
     if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
     for(int j=1;j<(Ny_+1)/2;j++) {
       size_t ip = k+NCz_*(j+Ny_*i);
       size_t ip1 = k+NCz_*((Ny_-j)+Ny_*i);
       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
       nconj1++;
     }
   }
 }
 cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;

 // b./ les lignes et colonnes hors continu et de nyquist
 int nconj2 = 0;
 for(int kk=0;kk<2;kk++) {
   int k=0;  // continu
   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
   for(int j=1;j<(Ny_+1)/2;j++) {
     if(Ny_%2==0 && j==Ny_/2) continue; // on ne retraite pas nyquist en j
     for(int i=1;i<Nx_;i++) {
       if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
       size_t ip = k+NCz_*(j+Ny_*i);
       size_t ip1 = k+NCz_*((Ny_-j)+Ny_*(Nx_-i));
       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
       nconj2++;
     }
   }
 }
 cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;

 cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;

 return nreal+nconj1+nconj2;
}

double GeneFluct3D::compute_power_carte(void)
// Calcul la puissance de la realisation du spectre Pk
{
 double s2 = 0.;
 for(int l=0;l<NCz_;l++)
   for(int j=0;j<Ny_;j++)
     for(int i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));

 double s20 = 0.;
 for(int j=0;j<Ny_;j++)
   for(int i=0;i<Nx_;i++) s20 += MODULE2(T_(0,j,i));

 double s2n = 0.;
 if(Nz_%2==0)
   for(int j=0;j<Ny_;j++)
     for(int i=0;i<Nx_;i++) s2n += MODULE2(T_(NCz_-1,j,i));

 return 2.*s2 -s20 -s2n;
}

//-------------------------------------------------------------------
void GeneFluct3D::FilterByPixel(void)
// Filtrage par la fonction fenetre du pixel (parallelepipede)
// TF = 1/(dx*dy*dz)*Int[{0,dx},{0,dy},{0,dz}]
//                   e^(ik_x*x) e^(ik_y*y) e^(ik_z*z) dxdydz
// |TF|^2 = 2*(1-cos(k_x*dx)/dx^2/k_x^2 * (idem dy) * (idem dz)
// et on traite la divergence en K_x = 0:
// 2*(1-cos(k*d)/d^2/k^2 ~= 1 - (k*d)^2/12 *[1 - (k*d)^2/30]
{
 int lp=2;
 if(lp>1) PrtTim("--- FilterByPixel: before ---");

 for(int i=0;i<Nx_;i++) {
   int ii = (i>Nx_/2) ? Nx_-i : i;
   double kx = ii*Dkx_ *Dx_;
   double pkx2 = pixelfilter(kx); 
   for(int j=0;j<Ny_;j++) {
     int jj = (j>Ny_/2) ? Ny_-j : j;
     double ky = jj*Dky_ *Dy_;
     double pky2 =  pixelfilter(ky); 
     for(int l=0;l<NCz_;l++) {
       double kz = l*Dkz_ *Dz_;
       double pkz2 =  pixelfilter(kz); 
       T_(l,j,i) *= pkx2*pky2*pkz2;
     }
   }
 }

 if(lp>1) PrtTim("--- FilterByPixel: after ---");
}

//-------------------------------------------------------------------
void GeneFluct3D::ComputeReal(void)
// Calcule une realisation dans l'espace reel
{
 int lp=2;

 if( Tcontent_==0 ) {
   cout<<"GeneFluct3D::ComputeReal_Error: empty array"<<endl;
   throw ParmError("GeneFluct3D::ComputeReal_Error: empty array");
 }

 // On fait la FFT
 if(lp>1) PrtTim("--- ComputeReal: before fftw backward ---");
 fftw_execute(pb_);
 if(lp>1) PrtTim("--- ComputeReal: after fftw backward ---");

 Tcontent_ = 2;
}

//-------------------------------------------------------------------
void GeneFluct3D::ReComputeFourier(void)
{
 int lp=2;

 // On fait la FFT
 if(lp>1) PrtTim("--- ComputeReal: before fftw forward ---");
 fftw_execute(pf_);
 // On corrige du pb de la normalisation de FFTW3
 double v = (double)NRtot_;
 for(int i=0;i<Nx_;i++)
   for(int j=0;j<Ny_;j++)
     for(int l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);

 Tcontent_ = 3;
 if(lp>1) PrtTim("--- ComputeReal: after fftw forward ---");
}

//-------------------------------------------------------------------
int GeneFluct3D::ComputeSpectrum(HProf& hp)
// Compute spectrum from "T" and fill HProf "hp"
// T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
// cad T(kz,ky,kx) avec  0<kz<kz_nyq  -ky_nyq<ky<ky_nyq  -kx_nyq<kx<kx_nyq
{

 if(hp.NBins()<0) return -1;
 hp.Zero();
 hp.SetErrOpt(false);

 // Attention a l'ordre
 for(int i=0;i<Nx_;i++) {
   int ii = (i>Nx_/2) ? Nx_-i : i;
   double kx = ii*Dkx_;  kx *= kx;
   for(int j=0;j<Ny_;j++) {
     int jj = (j>Ny_/2) ? Ny_-j : j;
     double ky = jj*Dky_;  ky *= ky;
     for(int l=0;l<NCz_;l++) {
       double kz = l*Dkz_;  kz *= kz;
       double k = sqrt(kx+ky+kz);
       double pk = MODULE2(T_(l,j,i));
       hp.Add(k,pk);
     }
   }
 }
 hp.UpdateHisto();

 // renormalize to directly compare to original spectrum
 double norm = Vol_;
 hp *= norm;

 return 0;
}

//-------------------------------------------------------
size_t GeneFluct3D::VarianceFrReal(double R,double& var)
// Recompute MASS variance in spherical top-hat (rayon=R)
{
 int lp=2;
 if(lp>1) PrtTim("--- VarianceFrReal: before computation ---");

 double *data = (double *) (&T_(0,0,0));
 int dnx = int(R/Dx_+0.5); if(dnx<=0) dnx = 1;
 int dny = int(R/Dy_+0.5); if(dny<=0) dny = 1;
 int dnz = int(R/Dz_+0.5); if(dnz<=0) dnz = 1;
 cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;

 double sum=0., sum2=0., r2 = R*R; size_t nsum=0;

 for(int i=dnx;i<Nx_-dnx;i+=dnx) {
   for(int j=dny;j<Ny_-dny;j+=dny) {
     for(int l=dnz;l<Nz_-dnz;l+=dnz) {
       double s=0.; size_t n=0;
       for(int ii=i-dnx;ii<=i+dnx;ii++) {
         double x = (ii-i)*Dx_; x *= x;
         for(int jj=j-dny;jj<=j+dny;jj++) {
           double y = (jj-j)*Dy_; y *= y;
           for(int ll=l-dnz;ll<=l+dnz;ll++) {
             double z = (ll-l)*Dz_; z *= z;
             if(x+y+z>r2) continue;
             size_t ip = ll+NTz_*(jj+Ny_*ii);
             s += 1.+data[ip];
             n++;
           }
         }
       }
       if(n>0) {sum += s; sum2 += s*s; nsum++;}
       //cout<<i<<","<<j<<","<<l<<" n="<<n<<" s="<<s<<" sum="<<sum<<" sum2="<<sum2<<endl;
     }
   }
 }

 if(nsum<=1) {var=0.; return nsum;}

 sum /= nsum;
 sum2 = sum2/nsum - sum*sum;
 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 ---");
 return nsum;
}

//-------------------------------------------------------
size_t GeneFluct3D::NumberOfBad(double vmin,double vmax)
// number of pixels outside of ]vmin,vmax[ extremites exclues
//     ->  vmin and vmax are considered as bad
{
 double *data = (double *) (&T_(0,0,0));

 size_t nbad = 0;
 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   double v = data[ip];
   if(v<=vmin || v>=vmax) nbad++;
 }

 return nbad;
}

size_t GeneFluct3D::MeanSigma2(double& rm,double& rs2,double vmin,double vmax)
// mean,sigma^2 pour pixels avec valeurs ]vmin,vmax[ extremites exclues
//   -> mean and sigma^2 are NOT computed with pixels values vmin and vmax
{
 double *data = (double *) (&T_(0,0,0));

 size_t n = 0;
 rm = rs2 = 0.;

 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   double v = data[ip];
   if(v<=vmin || v>=vmax) continue;
   rm += v;
   rs2 += v*v;
   n++;
 }

 if(n>1) {
   rm /= (double)n;
   rs2 = rs2/(double)n - rm*rm;
 }

 return n;
}

size_t GeneFluct3D::SetToVal(double vmin, double vmax,double val0)
// set to "val0" if out of range ]vmin,vmax[ extremites exclues
//     ->  vmin and vmax are set to val0
{
 double *data = (double *) (&T_(0,0,0));

 size_t nbad = 0;
 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   double v = data[ip];
   if(v<=vmin || v>=vmax) {data[ip] = val0; nbad++;}
 }

 return nbad;
}

//-------------------------------------------------------
void GeneFluct3D::TurnFluct2Mass(void)
// d_rho/rho -> Mass  (add one!)
{
 int lp=2;
 if(lp>1) PrtTim("--- TurnFluct2Mass: before computation ---");
 double *data = (double *) (&T_(0,0,0));

 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   data[ip] += 1.;
 }

 Tcontent_ = 4;
 if(lp>1) PrtTim("--- TurnFluct2Mass: after computation ---");
}

double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
// do NOT treate negative or nul values
{
 int lp=2;
 if(lp>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");

 double *data = (double *) (&T_(0,0,0));

 double m,s2;
 size_t ngood = MeanSigma2(m,s2,0.,1e+200);
 if(lp) cout<<"TurnMass2MeanNumber: ngood="<<ngood
            <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;

 // On doit mettre n*Vol galaxies dans notre survey
 // On en met uniquement dans les pixels de masse >0.
 // On NE met PAS a zero les pixels <0
 // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies
 //   comme on ne prend que les pixels >0, on doit normaliser
 //   a la moyenne de <1+d_rho/rho> sur ces pixels
 //   (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;
 double sum = 0.;
 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   data[ip] *= dn; // par coherence on multiplie aussi les <=0
   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;

 Tcontent_ = 6;
 if(lp>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
 return sum;
}

double GeneFluct3D::ApplyPoisson(void)
// do NOT treate negative or nul mass  -> let it as it is
{
 int lp=2;
 if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");

 double *data = (double *) (&T_(0,0,0));

 double sum = 0.;
 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   double v = data[ip];
   if(v>0.) {
     unsigned long dn = PoissRandLimit(v,10.);
     data[ip] = (double)dn;
     sum += (double)dn;
   }
 }
 if(lp) cout<<sum<<" galaxies put into survey"<<endl;
 Tcontent_ = 8;

 if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
 return sum;
}

double GeneFluct3D::TurnNGal2Mass(FunRan& massdist,bool axeslog)
// do NOT treate negative or nul mass  -> let it as it is
// INPUT:
//   massdist : distribution de masse (m*dn/dm)
//   axeslog = false : retourne la masse
//           = true  : retourne le log10(mass)
// RETURN la masse totale
{
 int lp=2;
 if(lp>1) PrtTim("--- TurnNGal2Mass: before computation ---");

  double *data = (double *) (&T_(0,0,0));

 double sum = 0.;
 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   double v = data[ip];
   if(v>0.) {
     int ngal = int(v+0.1);
     data[ip] = 0.;
     for(int i=0;i<ngal;i++) {
       double m = massdist.RandomInterp();  // massdist.Random();
       if(axeslog) m = pow(10.,m);
       data[ip] += m;
     }
     sum += data[ip];
   }
 }
 if(lp) cout<<sum<<" MSol HI mass put into survey"<<endl;
 Tcontent_ = 10;

 if(lp>1) PrtTim("--- TurnNGal2Mass: after computation ---");
 return sum;
}

void GeneFluct3D::AddNoise2Real(double snoise)
// add noise to every pixels (meme les <=0 !)
{
 int lp=2;
 if(lp>1) PrtTim("--- AddNoise2Real: before computation ---");

 double *data = (double *) (&T_(0,0,0));

 for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
   size_t ip = l+NTz_*(j+Ny_*i);
   data[ip] += snoise*NorRand();
 }
 Tcontent_ = 12;

 if(lp>1) PrtTim("--- AddNoise2Real: after computation ---");
}