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

#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 "fabtcolread.h"
#include "fabtwriter.h"
#include "fioarr.h"

#include "arrctcast.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)
  : T_(T) , array_allocated_(false) , lp_(0)
{
 SetNThread();
}

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

//-------------------------------------------------------
void GeneFluct3D::SetSize(long nx,long ny,long nz,double dx,double dy,double dz)
{
  setsize(nx,ny,nz,dx,dy,dz);
  setalloc();
  setpointers(false);
}

void GeneFluct3D::setsize(long nx,long ny,long 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)");
 }

 // Les tailles des tableaux
 Nx_ = nx;
 Ny_ = ny;  if(Ny_ <= 0) Ny_ = Nx_;
 Nz_ = nz;  if(Nz_ <= 0) Nz_ = Nx_;
 N_.resize(0); N_.push_back(Nx_); N_.push_back(Ny_); N_.push_back(Nz_);
 NRtot_ = Nx_*Ny_*Nz_; // nombre de pixels dans le survey
 NCz_ =  Nz_/2 +1;
 NTz_ = 2*NCz_;

 // le pas dans l'espace (Mpc)
 Dx_ = dx;
 Dy_ = dy; if(Dy_ <= 0.) Dy_ = Dx_;
 Dz_ = dz; if(Dz_ <= 0.) Dz_ = Dx_;
 D_.resize(0); D_.push_back(Dx_); D_.push_back(Dy_); D_.push_back(Dz_);
 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_);
 Dk_.resize(0); Dk_.push_back(Dkx_); Dk_.push_back(Dky_); Dk_.push_back(Dkz_);
 Dk3_ = Dkx_*Dky_*Dkz_;
 
 // La frequence de Nyquist en k (Mpc^-1)
 Knyqx_ = M_PI/Dx_;
 Knyqy_ = M_PI/Dy_;
 Knyqz_ = M_PI/Dz_;
 Knyq_.resize(0); Knyq_.push_back(Knyqx_); Knyq_.push_back(Knyqy_); Knyq_.push_back(Knyqz_);
}

void GeneFluct3D::setalloc(void)
{
 // Dimensionnement du tableau complex<r_8>
 // ATTENTION: TArray adresse en memoire a l'envers du C
 //            Tarray(n1,n2,n3) == Carray[n3][n2][n1]
 sa_size_t SzK_[3] = {NCz_,Ny_,Nx_}; // a l'envers
 try {
   T_.ReSize(3,SzK_);
   array_allocated_ = true;
 } catch (...) {
   cout<<"GeneFluct3D::SetSize_Error: Problem allocating T_"<<endl;
 }
 T_.SetMemoryMapping(BaseArray::CMemoryMapping);
}

void GeneFluct3D::setpointers(bool from_real)
{
 if(from_real) T_ = ArrCastR2C(R_);
   else        R_ = ArrCastC2R(T_);
 // On remplit les pointeurs
 fdata_ = (fftw_complex *) (&T_(0,0,0));
 data_ = (double *) (&R_(0,0,0));
}

void GeneFluct3D::check_array_alloc(void)
// Pour tester si le tableau T_ est alloue
{
 if(array_allocated_) return;
 char bla[90];
 sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
 cout<<bla<<endl;
 throw ParmError(bla);
}


//-------------------------------------------------------
void GeneFluct3D::WriteFits(string cfname,int bitpix)
{
 cout<<"GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;
 try {
   FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5);
   fwrt.WriteKey("NX",Nx_," axe transverse 1");
   fwrt.WriteKey("NY",Ny_," axe transverse 2");
   fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)");
   fwrt.WriteKey("DX",Dx_," Mpc");
   fwrt.WriteKey("DY",Dy_," Mpc");
   fwrt.WriteKey("DZ",Dz_," Mpc");
   fwrt.WriteKey("DKX",Dkx_," Mpc^-1");
   fwrt.WriteKey("DKY",Dky_," Mpc^-1");
   fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
   fwrt.Write(R_);
 } catch (PThrowable & exc) {
   cout<<"Exception : "<<(string)typeid(exc).name()
       <<" - Msg= "<<exc.Msg()<<endl;
   return;
 } catch (...) {
   cout<<" some other exception was caught !"<<endl;
   return;
 }
}

void GeneFluct3D::ReadFits(string cfname)
{
 cout<<"GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;
 try {
   FitsImg3DRead fimg(cfname.c_str(),0,5);
   fimg.Read(R_);
   long nx = fimg.ReadKeyL("NX");
   long ny = fimg.ReadKeyL("NY");
   long nz = fimg.ReadKeyL("NZ");
   double dx = fimg.ReadKey("DX");
   double dy = fimg.ReadKey("DY");
   double dz = fimg.ReadKey("DZ");
   setsize(nx,ny,nz,dx,dy,dz);
   setpointers(true);
 } catch (PThrowable & exc) {
   cout<<"Exception : "<<(string)typeid(exc).name()
       <<" - Msg= "<<exc.Msg()<<endl;
   return;
 } catch (...) {
   cout<<" some other exception was caught !"<<endl;
   return;
 }
}

void GeneFluct3D::WritePPF(string cfname,bool write_real)
// On ecrit soit le TArray<r_8> ou le TArray<complex <r_8> >
{
 cout<<"GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;
 try {
   R_.Info()["NX"] = (int_8)Nx_;
   R_.Info()["NY"] = (int_8)Ny_;
   R_.Info()["NZ"] = (int_8)Nz_;
   R_.Info()["DX"] = (r_8)Dx_;
   R_.Info()["DY"] = (r_8)Dy_;
   R_.Info()["DZ"] = (r_8)Dz_;
   POutPersist pos(cfname.c_str());
   if(write_real) pos << PPFNameTag("rgen")  << R_;
     else         pos << PPFNameTag("pkgen") << T_;
 } catch (PThrowable & exc) {
   cout<<"Exception : "<<(string)typeid(exc).name()
       <<" - Msg= "<<exc.Msg()<<endl;
   return;
 } catch (...) {
   cout<<" some other exception was caught !"<<endl;
   return;
 }
}

void GeneFluct3D::ReadPPF(string cfname)
{
 cout<<"GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;
 try {
   bool from_real = true;
   PInPersist pis(cfname.c_str());
   string name_tag_k = "pkgen";
   bool found_tag_k = pis.GotoNameTag("pkgen");
   if(found_tag_k) {
     cout<<"           ...reading spectrun into TArray<complex <r_8> >"<<endl;
     pis >> PPFNameTag("pkgen")  >> T_;
     from_real = false;
   } else {
     cout<<"           ...reading space into TArray<r_8>"<<endl;
     pis >> PPFNameTag("rgen")  >> R_;
   }
   int_8 nx = R_.Info()["NX"];
   int_8 ny = R_.Info()["NY"];
   int_8 nz = R_.Info()["NZ"];
   r_8 dx = R_.Info()["DX"];
   r_8 dy = R_.Info()["DY"];
   r_8 dz = R_.Info()["DZ"];
   setsize(nx,ny,nz,dx,dy,dz);
   setpointers(from_real);
 } catch (PThrowable & exc) {
   cout<<"Exception : "<<(string)typeid(exc).name()
       <<" - Msg= "<<exc.Msg()<<endl;
   return;
 } catch (...) {
   cout<<" some other exception was caught !"<<endl;
   return;
 }
}

//-------------------------------------------------------
void GeneFluct3D::Print(void)
{
 cout<<"GeneFluct3D(T_alloc="<<array_allocated_<<"):"<<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(GenericFunc& pk_at_z)
// cf ComputeFourier() mais avec autre methode de realisation du spectre
//    (attention on fait une fft pour realiser le spectre)
{
 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
 if(lp_>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
#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(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
   int_8 ip = IndexR(i,j,l);
   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
 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;
   for(long j=0;j<Ny_;j++) {
     long jj = (j>Ny_/2) ? Ny_-j : j;
     double ky = jj*Dky_;  ky *= ky;
     for(long 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 = pk_at_z(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 ---");

}

//-------------------------------------------------------
void GeneFluct3D::ComputeFourier(GenericFunc& pk_at_z)
// Calcule une realisation du spectre "pk_at_z"
// 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
{
 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
 if(lp_>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
#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 ---");
 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;
   for(long j=0;j<Ny_;j++) {
     long jj = (j>Ny_/2) ? Ny_-j : j;
     double ky = jj*Dky_;  ky *= ky;
     for(long 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 = pk_at_z(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 ---");

}

long 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
{
 check_array_alloc();

 // 1./ Le Continu et Nyquist sont reels
 long nreal = 0;
 for(long kk=0;kk<2;kk++) {
   long k=0;  // continu
   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
   for(long jj=0;jj<2;jj++) {
     long j=0;
     if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
     for(long ii=0;ii<2;ii++) {
       long i=0;
       if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
       int_8 ip = IndexC(i,j,k);
       //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
 long nconj1 = 0;
 for(long kk=0;kk<2;kk++) {
   long k=0;  // continu
   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
   for(long jj=0;jj<2;jj++) { // selon j
     long j=0;
     if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
     for(long i=1;i<(Nx_+1)/2;i++) {
       int_8 ip = IndexC(i,j,k);
       int_8 ip1 = IndexC(Nx_-i,j,k);
       fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
       nconj1++;
     }
   }
   for(long ii=0;ii<2;ii++) {
     long i=0;
     if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
     for(long j=1;j<(Ny_+1)/2;j++) {
       int_8 ip = IndexC(i,j,k);
       int_8 ip1 = IndexC(i,Ny_-j,k);
       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
 long nconj2 = 0;
 for(long kk=0;kk<2;kk++) {
   long k=0;  // continu
   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
   for(long j=1;j<(Ny_+1)/2;j++) {
     if(Ny_%2==0 && j==Ny_/2) continue; // on ne retraite pas nyquist en j
     for(long i=1;i<Nx_;i++) {
       if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
       int_8 ip = IndexC(i,j,k);
       int_8 ip1 = IndexC(Nx_-i,Ny_-j,k);
       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
{
 check_array_alloc();

 double s2 = 0.;
 for(long l=0;l<NCz_;l++)
   for(long j=0;j<Ny_;j++)
     for(long i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));

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

 double s2n = 0.;
 if(Nz_%2==0)
   for(long j=0;j<Ny_;j++)
     for(long 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[{-dx/2,dx/2},{-dy/2,dy/2},{-dz/2,dz/2}]
//                   e^(ik_x*x) e^(ik_y*y) e^(ik_z*z) dxdydz
//    = 2/(k_x*dx) * sin(k_x*dx/2)  * (idem y) * (idem z)
// Gestion divergence en 0: sin(y)/y = 1 - y^2/6*(1-y^2/20)
//                          avec y = k_x*dx/2
{
 check_array_alloc();

 if(lp_>1) PrtTim("--- FilterByPixel: before ---");

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

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

//-------------------------------------------------------------------
void GeneFluct3D::ComputeReal(void)
// Calcule une realisation dans l'espace reel
{
 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 ---");
}

//-------------------------------------------------------------------
void GeneFluct3D::ReComputeFourier(void)
{
 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
 double v = (double)NRtot_;
 for(long i=0;i<Nx_;i++)
   for(long j=0;j<Ny_;j++)
     for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);

 if(lp_>1) PrtTim("--- ComputeReal: after fftw forward ---");
}

//-------------------------------------------------------------------
int GeneFluct3D::ComputeSpectrum(HistoErr& herr)
// Compute spectrum from "T" and fill HistoErr "herr"
// 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
{
 check_array_alloc();
 if(lp_>1) PrtTim("--- ComputeSpectrum: before ---");

 if(herr.NBins()<0) return -1;
 herr.Zero();

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

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

 if(lp_>1) PrtTim("--- ComputeSpectrum: after ---");

 return 0;
}

int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr)
{
 check_array_alloc();
 if(lp_>1) PrtTim("--- ComputeSpectrum2D: before ---");

 if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
 herr.Zero();

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

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

 if(lp_>1) PrtTim("--- ComputeSpectrum2D: after ---");

 return 0;
}

//-------------------------------------------------------
int_8 GeneFluct3D::VarianceFrReal(double R,double& var)
// Recompute MASS variance in spherical top-hat (rayon=R)
{
 check_array_alloc();

 if(lp_>1) PrtTim("--- VarianceFrReal: before computation ---");

 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;

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

 for(long i=dnx;i<Nx_-dnx;i+=dnx) {
   for(long j=dny;j<Ny_-dny;j+=dny) {
     for(long l=dnz;l<Nz_-dnz;l+=dnz) {
       double s=0.; int_8 n=0;
       for(long ii=i-dnx;ii<=i+dnx;ii++) {
         double x = (ii-i)*Dx_; x *= x;
         for(long jj=j-dny;jj<=j+dny;jj++) {
           double y = (jj-j)*Dy_; y *= y;
           for(long ll=l-dnz;ll<=l+dnz;ll++) {
             double z = (ll-l)*Dz_; z *= z;
             if(x+y+z>r2) continue;
             int_8 ip = IndexR(ii,jj,ll);
             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;
}

//-------------------------------------------------------
int_8 GeneFluct3D::NumberOfBad(double vmin,double vmax)
// number of pixels outside of ]vmin,vmax[ extremites exclues
//     ->  vmin and vmax are considered as bad
{
 check_array_alloc();

 int_8 nbad = 0;
 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
   int_8 ip = IndexR(i,j,l);
   double v = data_[ip];
   if(v<=vmin || v>=vmax) nbad++;
 }

 return nbad;
}

int_8 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
{
 check_array_alloc();

 int_8 n = 0;
 rm = rs2 = 0.;

 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
   int_8 ip = IndexR(i,j,l);
   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;
}

int_8 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
{
 check_array_alloc();

 int_8 nbad = 0;
 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
   int_8 ip = IndexR(i,j,l);
   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!)
{
 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++) {
   int_8 ip = IndexR(i,j,l);
   data_[ip] += 1.;
 }

 if(lp_>1) PrtTim("--- TurnFluct2Mass: after computation ---");
}

double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
// do NOT treate negative or nul values
{
 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
            <<" 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(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
   int_8 ip = IndexR(i,j,l);
   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;

 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
{
 check_array_alloc();

 if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");

 double sum = 0.;
 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
   int_8 ip = IndexR(i,j,l);
   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;

 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
{
 check_array_alloc();

 if(lp_>1) PrtTim("--- TurnNGal2Mass: before computation ---");

 double sum = 0.;
 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
   int_8 ip = IndexR(i,j,l);
   double v = data_[ip];
   if(v>0.) {
     long ngal = long(v+0.1);
     data_[ip] = 0.;
     for(long 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;

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

void GeneFluct3D::AddNoise2Real(double snoise)
// add noise to every pixels (meme les <=0 !)
{
 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++) {
   int_8 ip = IndexR(i,j,l);
   data_[ip] += snoise*NorRand();
 }

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