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

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

#include "fabtcolread.h"
#include "fabtwriter.h"
#include "fioarr.h"
#include "ntuple.h"

#include "arrctcast.h"

#include "constcosmo.h"
#include "geneutils.h"
#include "schechter.h"

#include "genefluct3d.h"

#if defined(GEN3D_FLOAT)
#define GEN3D_FFTW_INIT_THREADS       fftwf_init_threads
#define GEN3D_FFTW_CLEANUP_THREADS    fftwf_cleanup_threads
#define GEN3D_FFTW_PLAN_WITH_NTHREADS fftwf_plan_with_nthreads
#define GEN3D_FFTW_PLAN_DFT_R2C_3D    fftwf_plan_dft_r2c_3d
#define GEN3D_FFTW_PLAN_DFT_C2R_3D    fftwf_plan_dft_c2r_3d
#define GEN3D_FFTW_DESTROY_PLAN       fftwf_destroy_plan
#define GEN3D_FFTW_EXECUTE            fftwf_execute
#else
#define GEN3D_FFTW_INIT_THREADS       fftw_init_threads
#define GEN3D_FFTW_CLEANUP_THREADS    fftw_cleanup_threads
#define GEN3D_FFTW_PLAN_WITH_NTHREADS fftw_plan_with_nthreads
#define GEN3D_FFTW_PLAN_DFT_R2C_3D    fftw_plan_dft_r2c_3d
#define GEN3D_FFTW_PLAN_DFT_C2R_3D    fftw_plan_dft_c2r_3d
#define GEN3D_FFTW_DESTROY_PLAN       fftw_destroy_plan
#define GEN3D_FFTW_EXECUTE            fftw_execute
#endif

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

namespace SOPHYA {

//-------------------------------------------------------
GeneFluct3D::GeneFluct3D(long nx,long ny,long nz,double dx,double dy,double dz
                        ,unsigned short nthread,int lp)
{
 init_default();

 lp_ = lp;
 nthread_ = nthread;

 setsize(nx,ny,nz,dx,dy,dz);
 setalloc();
 setpointers(false);
 init_fftw();
}

GeneFluct3D::GeneFluct3D(unsigned short nthread)
{
 init_default();
 setsize(2,2,2,1.,1.,1.);
 nthread_ = nthread;
 setalloc();
 setpointers(false);
 init_fftw();
}

GeneFluct3D::~GeneFluct3D(void)
{
 delete_fftw();
}

//-------------------------------------------------------
void GeneFluct3D::init_default(void)
{
 Nx_ = Ny_ = Nz_ = 0;
 is_set_fft_plan = false;
 nthread_ = 0;
 lp_ = 0;
 array_allocated_ = false; array_type = 0;
 cosmo_ = NULL;
 growth_ = NULL;
 good_dzinc_ = 0.01;
 compute_pk_redsh_ref_ = -999.;
 redsh_ref_ = -999.;
 kredsh_ref_ = 0.;
 dred_ref_ = -999.;
 h_ref_ = -999.;
 loscom_ref_ = -999.;
 dtrc_ref_ = dlum_ref_ = dang_ref_ = -999.;
 nu_ref_ = dnu_ref_ = -999.;
 loscom_min_ = loscom_max_ = -999.;
 loscom2zred_min_ = loscom2zred_max_ = 0.;
 xobs_[0] = xobs_[1] = xobs_[2] = 0.;
 zred_.resize(0);
 loscom_.resize(0);
 loscom2zred_.resize(0);
}

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.) {
   const char *bla = "GeneFluct3D::setsize_Error: bad value(s) for nn/dx";
   cout<<bla<<endl; throw ParmError(bla);
 }

 // 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_ = (int_8)Nx_*(int_8)Ny_*(int_8)Nz_; // nombre de pixels dans le survey
 NCz_ =  Nz_/2 +1;
 NFcoef_ = (int_8)Nx_*(int_8)Ny_*(int_8)NCz_; // nombre de coeff de Fourier dans le survey
 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)
{
#if defined(GEN3D_FLOAT)
 if(lp_>1) cout<<"--- GeneFluct3D::setalloc FLOAT ---"<<endl;
#else
 if(lp_>1) cout<<"--- GeneFluct3D::setalloc DOUBLE ---"<<endl;
#endif
 // 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; array_type=0;
   if(lp_>1) cout<<"  allocating: "<<T_.Size()*sizeof(complex<GEN3D_TYPE>)/1.e6<<" Mo"<<endl;
 } catch (...) {
   cout<<"GeneFluct3D::setalloc_Error: Problem allocating T_"<<endl;
 }
}

void GeneFluct3D::setpointers(bool from_real)
{
 if(lp_>1) cout<<"--- GeneFluct3D::setpointers ---"<<endl;
 if(from_real) T_ = ArrCastR2C(R_);
   else        R_ = ArrCastC2R(T_);
 // On remplit les pointeurs
 fdata_ = (GEN3D_FFTW_COMPLEX *) (&T_(0,0,0));
 data_ = (GEN3D_TYPE *) (&R_(0,0,0));
}

void GeneFluct3D::init_fftw(void)
{
 if( is_set_fft_plan ) delete_fftw();

 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
 if(lp_>1) cout<<"--- GeneFluct3D::init_fftw ---"<<endl;
#ifdef WITH_FFTW_THREAD
 if(nthread_>0) {
   cout<<"...Computing with "<<nthread_<<" threads"<<endl;
   GEN3D_FFTW_INIT_THREADS();
   GEN3D_FFTW_PLAN_WITH_NTHREADS(nthread_);
 }
#endif
 if(lp_>1) cout<<"...forward plan"<<endl;
 pf_ = GEN3D_FFTW_PLAN_DFT_R2C_3D(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
 if(lp_>1) cout<<"...backward plan"<<endl;
 pb_ = GEN3D_FFTW_PLAN_DFT_C2R_3D(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
 is_set_fft_plan = true;
}

void GeneFluct3D::delete_fftw(void)
{
 if( !is_set_fft_plan ) return;
 GEN3D_FFTW_DESTROY_PLAN(pf_);
 GEN3D_FFTW_DESTROY_PLAN(pb_);
#ifdef WITH_FFTW_THREAD
 if(nthread_>0) GEN3D_FFTW_CLEANUP_THREADS();
#endif
 is_set_fft_plan = false;
}

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::SetObservator(double redshref,double kredshref)
// L'observateur est au redshift z=0
//               est situe sur la "perpendiculaire" a la face x,y
//                         issue du centre de cette face
// Il faut positionner le cube sur l'axe des z cad des redshifts:
//     redshref  = redshift de reference
//                 Si redshref<0 alors redshref=0
//     kredshref = indice (en double) correspondant a ce redshift
//                 Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
// Exemple: redshref=1.5 kredshref=250.75
//    -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
{
 if(redshref<0.) redshref = 0.;
 if(kredshref<0.) {
   if(Nz_<=0) {
     const char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 define cube geometry first";
     cout<<bla<<endl; throw ParmError(bla);
   }
   kredshref = Nz_/2.;
 }
 redsh_ref_  = redshref;
 kredsh_ref_ = kredshref;
 if(lp_>0)
   cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
}

void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
{
 cosmo_ = &cosmo;
 if(lp_>1) cosmo_->Print();
}

void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
{
 growth_ = &growth;
}

long GeneFluct3D::LosComRedshift(double zinc,long npoints)
// Given a position of the cube relative to the observer
// and a cosmology
// (SetObservator() and SetCosmology() should have been called !)
// This routine filled:
//   the vector "zred_" of scanned redshift (by zinc increments)
//   the vector "loscom_" of corresponding los comoving distance
// -- Input:
// zinc : redshift increment for computation
// npoints : number of points required for inverting loscom -> zred
// 
{
 if(zinc<=0.) zinc = good_dzinc_;
 if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;

 if(cosmo_ == NULL || growth_==NULL || redsh_ref_<0.) {
   const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator, Cosmology and Growth first";
   cout<<bla<<endl; throw ParmError(bla);
 }

 // La distance angulaire/luminosite/Dnu au pixel de reference
 dred_ref_ = Dz_/(cosmo_->Dhubble()/cosmo_->E(redsh_ref_));
 loscom_ref_ = cosmo_->Dloscom(redsh_ref_);
 dtrc_ref_ = cosmo_->Dtrcom(redsh_ref_);
 dlum_ref_ = cosmo_->Dlum(redsh_ref_);
 dang_ref_ = cosmo_->Dang(redsh_ref_);
 h_ref_ = cosmo_->H(redsh_ref_);
 growth_ref_ = (*growth_)(redsh_ref_);
 dsdz_growth_ref_ = growth_->DsDz(redsh_ref_,zinc);
 nu_ref_   = Fr_HyperFin_Par/(1.+redsh_ref_); // GHz
 dnu_ref_  = Fr_HyperFin_Par *dred_ref_/pow(1.+redsh_ref_,2.); // GHz
 if(lp_>0) {
   cout<<"...reference pixel redshref="<<redsh_ref_
       <<", dredref="<<dred_ref_
       <<", nuref="<<nu_ref_ <<" GHz"
       <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
       <<"   H="<<h_ref_<<" km/s/Mpc"
       <<", D="<<growth_ref_
       <<", dD/dz="<<dsdz_growth_ref_<<endl
       <<"   dlosc="<<loscom_ref_<<" Mpc com"
       <<", dtrc="<<dtrc_ref_<<" Mpc com"
       <<", dlum="<<dlum_ref_<<" Mpc"
       <<", dang="<<dang_ref_<<" Mpc"<<endl;
 }

 // On calcule les coordonnees de l'observateur dans le repere du cube
 // cad dans le repere ou l'origine est au centre du pixel i=j=l=0.
 // L'observateur est sur un axe centre sur le milieu de la face Oxy
 xobs_[0] = Nx_/2.*Dx_;
 xobs_[1] = Ny_/2.*Dy_;
 xobs_[2] = kredsh_ref_*Dz_ - loscom_ref_;

 // L'observateur est-il dans le cube?
 bool obs_in_cube = false;
 if(xobs_[2]>=0. && xobs_[2]<=Nz_*Dz_) obs_in_cube = true;

 // Find MINIMUM los com distance to the observer:
 // c'est le centre de la face a k=0
 // (ou zero si l'observateur est dans le cube)
 loscom_min_ = 0.;
 if(!obs_in_cube) loscom_min_ = -xobs_[2];

 // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED
 if(loscom_min_<=1.e-50)
   for(int i=0;i<10;i++)
     cout<<"ATTENTION TOUTES LES PARTIES DU CODE NE MARCHENT PAS POUR UN OBSERVATEUR DANS LE CUBE"<<endl;
 // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED


 // Find MAXIMUM los com distance to the observer:
 // ou que soit positionne l'observateur, la distance
 // maximal est sur un des coins du cube
 loscom_max_ = 0.;
 for(long i=0;i<=1;i++) {
   double dx2 = DXcom(i*(Nx_-1)); dx2 *= dx2;
   for(long j=0;j<=1;j++) {
     double dy2 = DYcom(j*(Ny_-1)); dy2 *= dy2;
     for(long k=0;k<=1;k++) {
       double dz2 = DZcom(k*(Nz_-1)); dz2 *= dz2;
       dz2 = sqrt(dx2+dy2+dz2);
       if(dz2>loscom_max_) loscom_max_ = dz2;
     }
   }
 }
 if(lp_>0) {
   cout<<"...zref="<<redsh_ref_<<" kzref="<<kredsh_ref_<<" losref="<<loscom_ref_<<" Mpc\n"
       <<"   xobs="<<xobs_[0]<<" , "<<xobs_[1]<<" , "<<xobs_[2]<<" Mpc "
       <<" in_cube="<<obs_in_cube
       <<" loscom_min="<<loscom_min_<<" loscom_max="<<loscom_max_<<" Mpc (com)"<<endl;
 }

 // Fill the corresponding vectors for loscom and zred
 // Be shure to have one dlc<loscom_min and one dlc>loscom_max
 double zmin = 0., dlcmin=0.;
 while(1) {
   if(lp_>0)
     cout<<"...Filling zred starting at zmin="<<zmin<<" with zinc="<<zinc
         <<", loscom_min-max=["<<loscom_min_<<","<<loscom_max_<<"]"<<endl;
   zred_.resize(0); loscom_.resize(0);
   for(double z=zmin; ; z+=zinc) {
     double dlc = cosmo_->Dloscom(z);
     if(dlc<loscom_min_) {
       dlcmin = dlc;
       zmin = z;
       zred_.resize(0); loscom_.resize(0);
     }
     zred_.push_back(z);
     loscom_.push_back(dlc);
     z += zinc;
     if(dlc>loscom_max_) {
       if(lp_>0)
         cout<<"  Min: z="<<zmin<<" dlc="<<dlcmin<<", Max: z="<<z<<" dlc="<<dlc<<endl;
       break; // on sort apres avoir stoque un dlc>dlcmax
     }
   }
   if(zred_.size()>=10) break;
   zinc /= 10.;
   cout<<"  not enough points ("<<zred_.size()<<") for zref, retry with zinc="<<zinc<<endl;
 }

 if(lp_>0) {
   long n = zred_.size();
   cout<<"...zred/loscom tables[zinc="<<zinc<<"]: n="<<n;
   if(n>0) cout<<" z="<<zred_[0]<<" -> d="<<loscom_[0];
   if(n>1) cout<<" , z="<<zred_[n-1]<<" -> d="<<loscom_[n-1];
   cout<<endl;
 }

 // Compute the parameters and tables needed for inversion loscom->zred
 if(npoints<3) npoints = zred_.size();
 InverseFunc invfun(zred_,loscom_);
 invfun.ComputeParab(npoints,loscom2zred_);
 loscom2zred_min_ = invfun.YMin();
 loscom2zred_max_ = invfun.YMax();

 if(lp_>0) {
   long n = loscom2zred_.size();
   cout<<"...loscom -> zred[npoints="<<npoints<<"]: n="<<n
       <<" los_min="<<loscom2zred_min_
       <<" los_max="<<loscom2zred_max_
       <<" -> zred=[";
   if(n>0) cout<<loscom2zred_[0];
   cout<<",";
   if(n>1) cout<<loscom2zred_[n-1];
   cout<<"]"<<endl;
   if(lp_>1 && n>0)
     for(int i=0;i<n;i++)
       if(i<2 || abs(i-n/2)<2 || i>=n-2)
         cout<<"    i="<<i
             <<"  d="<<loscom2zred_min_+i*(loscom2zred_max_-loscom2zred_min_)/(n-1.)
             <<" Mpc   z="<<loscom2zred_[i]<<endl;
 }


 // Compute the table for D'(z)/D(z) where D'(z)=dD/dEta (Eta is conformal time)
 zredmin_dpsd_ = zred_[0];
 zredmax_dpsd_ = zred_[zred_.size()-1];
 long nptd = long(sqrt(Nx_*Nx_ + Ny_*Ny_ +Nz_*Nz_));
 if(nptd<10) nptd = 10;
 good_dzinc_ = (zredmax_dpsd_ - zredmin_dpsd_)/nptd;
 if(lp_>0) cout<<"...good_dzinc changed to "<<good_dzinc_<<endl;
 dpsdfrzred_.resize(0);
 double dz = (zredmax_dpsd_ - zredmin_dpsd_)/(nptd-1);
 int nmod = nptd/5; if(nmod==0) nmod = 1;
 if(lp_>0) cout<<"...Compute table D'/D on "<<nptd<<" pts for z["
               <<zredmin_dpsd_<<","<<zredmax_dpsd_<<"]"<<endl;
 for(int i=0;i<nptd;i++) {
   double z = zredmin_dpsd_ + i*dz;
   double v = -cosmo_->H(z) * growth_->DsDz(z,good_dzinc_) / (*growth_)(z);
   dpsdfrzred_.push_back(v);
   if(lp_ && (i%nmod==0 || i==nptd-1)) cout<<"    z="<<z<<"  D'/D="<<v<<" km/s/Mpc"<<endl;
 }

 return zred_.size();
}

//-------------------------------------------------------
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");  // NAXIS3
   fwrt.WriteKey("NY",Ny_," axe transverse 2");  // NAXIS2
   fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)");  // NAXIS1
   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.WriteKey("ZREFPK",compute_pk_redsh_ref_," Pk computed redshift");
   fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
   fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
   fwrt.WriteKey("HREF",h_ref_," ref H(z)");
   fwrt.WriteKey("Growth",Dref()," growth at reference redshift");
   fwrt.WriteKey("GrowthPk",DrefPk()," growth at Pk computed redshift");
   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");  // NAXIS3
   long ny = fimg.ReadKeyL("NY");  // NAXIS2
   long nz = fimg.ReadKeyL("NZ");  // NAXIS1
   double dx = fimg.ReadKey("DX");
   double dy = fimg.ReadKey("DY");
   double dz = fimg.ReadKey("DZ");
   double pkzref = fimg.ReadKey("ZREFPK");
   double zref = fimg.ReadKey("ZREF");
   double kzref = fimg.ReadKey("KZREF");
   setsize(nx,ny,nz,dx,dy,dz);
   setpointers(true);
   init_fftw();
   SetObservator(zref,kzref);
   array_allocated_ = true;
   compute_pk_redsh_ref_ = pkzref;
 } 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_;
   R_.Info()["ZREFPK"] = (r_8)compute_pk_redsh_ref_;
   R_.Info()["ZREF"] = (r_8)redsh_ref_;
   R_.Info()["KZREF"] = (r_8)kredsh_ref_;
   R_.Info()["HREF"] = (r_8)h_ref_;
   R_.Info()["Growth"] = (r_8)Dref();
   R_.Info()["GrowthPk"] = (r_8)DrefPk();
   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 spectrum 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_;
   }
   setpointers(from_real);  // a mettre ici pour relire les DVInfo
   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"];
   r_8 pkzref = R_.Info()["ZREFPK"];
   r_8 zref = R_.Info()["ZREF"];
   r_8 kzref = R_.Info()["KZREF"];
   setsize(nx,ny,nz,dx,dy,dz);
   init_fftw();
   SetObservator(zref,kzref);
   array_allocated_ = true;
   compute_pk_redsh_ref_ = pkzref;
 } 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::WriteSlicePPF(string cfname)
// On ecrit 3 tranches du cube selon chaque axe
{
 cout<<"--- GeneFluct3D::WriteSlicePPF: Writing Cube Slices "<<cfname<<endl;
 try {

   POutPersist pos(cfname.c_str());
   TMatrix<r_4> S;
   char str[16];
   long i,j,l;

   // Tranches en Z
   for(int s=0;s<3;s++) {
     S.ReSize(Nx_,Ny_);
     if(s==0) l=0; else if(s==1) l=(Nz_+1)/2; else  l=Nz_-1;
     sprintf(str,"z%ld",l);
     for(i=0;i<Nx_;i++) for(j=0;j<Ny_;j++) S(i,j)=data_[IndexR(i,j,l)];
     pos<<PPFNameTag(str)<<S; S.RenewObjId();
   }

   // Tranches en Y
   for(int s=0;s<3;s++) {
     S.ReSize(Nz_,Nx_);
     if(s==0) j=0; else if(s==1) j=(Ny_+1)/2; else  j=Ny_-1;
     sprintf(str,"y%ld",j);
     for(i=0;i<Nx_;i++) for(l=0;l<Nz_;l++) S(l,i)=data_[IndexR(i,j,l)];
     pos<<PPFNameTag(str)<<S; S.RenewObjId();
   }

   // Tranches en X
   for(int s=0;s<3;s++) {
     S.ReSize(Nz_,Ny_);
     if(s==0) i=0; else if(s==1) i=(Nx_+1)/2; else  i=Nx_-1;
     sprintf(str,"x%ld",i);
     for(j=0;j<Ny_;j++) for(l=0;l<Nz_;l++) S(l,j)=data_[IndexR(i,j,l)];
     pos<<PPFNameTag(str)<<S; S.RenewObjId();
   }

 } 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::NTupleCheck(POutPersist &pos,string ntname,unsigned long nent)
// Remplit le NTuple "ntname" avec "nent" valeurs du cube (reel ou complex) et l'ecrit dans "pos"
{
  if(ntname.size()<=0 || nent==0) return;
  int nvar = 0;
  if(array_type==1) nvar = 3;
  else if(array_type==2) nvar = 4;
  else return;
  const char *vname[4] = {"t","z","re","im"};
  float xnt[4];
  NTuple nt(nvar,vname);

  if(array_type==1) {
    unsigned long nmod = Nx_*Ny_*Nz_/nent; if(nmod==0) nmod=1;
    unsigned long n=0;
    for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
      if(n==nmod) {
        int_8 ip = IndexR(i,j,l);
        xnt[0]=sqrt(i*i+j*j); xnt[1]=l; xnt[2]=data_[ip];
        nt.Fill(xnt);
        n=0;
      }
      n++;
    }
  } else {
    unsigned long nmod = Nx_*Ny_*NCz_/nent; if(nmod==0) nmod=1;
    unsigned long n=0;
    for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<NCz_;l++) {
      if(n==nmod) {
        xnt[0]=sqrt(i*i+j*j); xnt[1]=l; xnt[2]=T_(l,j,i).real(); xnt[3]=T_(l,j,i).imag();
        nt.Fill(xnt);
        n=0;
      }
      n++;
    }
  }

  pos.PutObject(nt,ntname);
}

//-------------------------------------------------------
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_<<" ncoeff="<<NFcoef_<<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;
 cout<<"Redshift "<<redsh_ref_<<" for z axe at k="<<kredsh_ref_<<endl;
}

//-------------------------------------------------------
void GeneFluct3D::ComputeFourier0(PkSpectrumZ& pk_at_z)
// cf ComputeFourier() mais avec autre methode de realisation du spectre
//    (attention on fait une fft pour realiser le spectre)
{
  compute_pk_redsh_ref_ = pk_at_z.GetZ();
 // --- realisation d'un tableau de tirage gaussiens
  if(lp_>0) cout<<"--- ComputeFourier0 at z="<<compute_pk_redsh_ref_
                <<": before gaussian filling ---"<<endl;
 // 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;
 }

 // --- realisation d'un tableau de tirage gaussiens
 if(lp_>0) cout<<"...before fft real ---"<<endl;
 GEN3D_FFTW_EXECUTE(pf_);

 // --- On remplit avec une realisation
 if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
 T_(0,0,0) = complex<GEN3D_TYPE>(0.);  // on coupe le continue et on l'initialise
 long lmod = Nx_/20; if(lmod<1) lmod=1;
 for(long i=0;i<Nx_;i++) {
   double kx = AbsKx(i);  kx *= kx;
   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
   for(long j=0;j<Ny_;j++) {
     double ky = AbsKy(j);  ky *= ky;
     for(long l=0;l<NCz_;l++) {
       double kz = AbsKz(l);  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);
     }
   }
 }

 array_type = 2;

 if(lp_>0) cout<<"...computing power"<<endl;
 double p = compute_power_carte();
 if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;

}

//-------------------------------------------------------
void GeneFluct3D::ComputeFourier(PkSpectrumZ& 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
{
 // --- RaZ du tableau
 T_ = complex<GEN3D_TYPE>(0.);
 compute_pk_redsh_ref_ = pk_at_z.GetZ();

 // --- On remplit avec une realisation
 if(lp_>0) cout<<"--- ComputeFourier at z="<<compute_pk_redsh_ref_<<" ---"<<endl;
 long lmod = Nx_/20; if(lmod<1) lmod=1;
 for(long i=0;i<Nx_;i++) {
   double kx = AbsKx(i);  kx *= kx;
   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
   for(long j=0;j<Ny_;j++) {
     double ky = AbsKy(j);  ky *= ky;
     for(long l=0;l<NCz_;l++) {
       double kz = AbsKz(l);  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) = ComplexGaussianRand(sqrt(pk/2.));
     }
   }
 }

 array_type = 2;
 manage_coefficients();   // gros effet pour les spectres que l'on utilise !

 if(lp_>0) cout<<"...computing power"<<endl;
 double p = compute_power_carte();
 if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;

}

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
// On ecrit que: conj(P(k_x,k_y,k_z)) = P(-k_x,-k_y,-k_z)
//               avec k_x = i, -k_x = N_x - i  etc...
{
 if(lp_>1) cout<<"...managing coefficients"<<endl;
 check_array_alloc();

 // 1./ Le Continu et Nyquist sont reels
 int_8 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++;
     }
   }
 }
 if(lp_>1) 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_8  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++;
     }
   }
 }
 if(lp_>1) cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;

 // b./ les lignes et colonnes hors continu et de nyquist
 int_8  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++;
     }
   }
 }
 if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;

 if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(NFcoef_-nconj1-nconj2)-nreal<<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
{
 if(lp_>0) cout<<"--- FilterByPixel ---"<<endl;
 check_array_alloc();

 for(long i=0;i<Nx_;i++) {
   double kx = AbsKx(i) *Dx_/2;
   double pk_x = pixelfilter(kx);
   for(long j=0;j<Ny_;j++) {
     double ky = AbsKy(j) *Dy_/2;
     double pk_y = pixelfilter(ky);
     for(long l=0;l<NCz_;l++) {
       double kz = AbsKz(l) *Dz_/2;
       double pk_z =  pixelfilter(kz);
       T_(l,j,i) *= pk_x*pk_y*pk_z;
     }
   }
 }

}

//-------------------------------------------------------------------
void GeneFluct3D::ToVelLoS(void)
// Le spectre Pk doit etre (dRho/Rho)(k)
{
 double zpk = compute_pk_redsh_ref_;
 double dpsd = -cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
 if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", D'/D="<<dpsd
               <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
 check_array_alloc();

 for(long i=0;i<Nx_;i++) {
   double kx = Kx(i);
   for(long j=0;j<Ny_;j++) {
     double ky = Ky(j);
     double kt2 = kx*kx + ky*ky;
     for(long l=0;l<NCz_;l++) {
       double kz = Kz(l);
       double k2 = kt2 + kz*kz;
       if(l==0) k2 = 0.; else k2 = dpsd*kz/k2;
       T_(l,j,i) *= complex<double>(0.,k2);
     }
   }
 }

}

//-------------------------------------------------------------------
void GeneFluct3D::ApplyGrowthFactor(int type_evol)
// Apply Growth to real space
// Using the correspondance between redshift and los comoving distance
// describe in vector "zred_" "loscom_"
// type_evol = 1 : evolution avec la distance a l'observateur
//             2 : evolution avec la distance du plan Z
//             (tous les pixels d'un plan Z sont mis au meme redshift z que celui du milieu)
{
 if(lp_>0) cout<<"--- ApplyGrowthFactor: evol="<<type_evol<<endl;
 check_array_alloc();

 if(growth_ == NULL) {
   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
   cout<<bla<<endl; throw ParmError(bla);
 }
 if(type_evol<1 || type_evol>2) {
   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: bad type_evol value";
   cout<<bla<<endl; throw ParmError(bla);
 }

 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);

 for(long i=0;i<Nx_;i++) {
   double dx2 = DXcom(i); dx2 *= dx2;
   for(long j=0;j<Ny_;j++) {
     double dy2 = DYcom(j); dy2 *= dy2;
     for(long l=0;l<Nz_;l++) {
       double dz = DZcom(l);
       if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
         else dz = fabs(dz); // tous les plans Z au meme redshift
       double z = interpinv(dz);   // interpolation par morceau
       double dzgr = (*growth_)(z);
       int_8 ip = IndexR(i,j,l);
       data_[ip] *= dzgr;
     }
   }
 }

}

//-------------------------------------------------------------------
void GeneFluct3D::ApplyDerGrowthFactor(int type_evol)
// Apply Conformal derivative of Growth to real space for transverse velocity cube
{
 if(lp_>0) cout<<"--- ApplyDerGrowthFactor: evol="<<type_evol<<endl;
 check_array_alloc();

 if(growth_ == NULL) {
   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
   cout<<bla<<endl; throw ParmError(bla);
 }
 if(type_evol<1 || type_evol>2) {
   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: bad type_evol value";
   cout<<bla<<endl; throw ParmError(bla);
 }

 double zpk = compute_pk_redsh_ref_;
 double dpsd_orig = - cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);

 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
 InterpFunc interdpd(zredmin_dpsd_,zredmax_dpsd_,dpsdfrzred_);

 for(long i=0;i<Nx_;i++) {
   double dx2 = DXcom(i); dx2 *= dx2;
   for(long j=0;j<Ny_;j++) {
     double dy2 = DYcom(j); dy2 *= dy2;
     for(long l=0;l<Nz_;l++) {
       double dz = DZcom(l);
       if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
         else dz = fabs(dz); // tous les plans Z au meme redshift
       double z = interpinv(dz);   // interpolation par morceau
       double dpsd = interdpd(z);
       int_8 ip = IndexR(i,j,l);
       data_[ip] *= dpsd / dpsd_orig;
     }
   }
 }

}

//-------------------------------------------------------------------
void GeneFluct3D::ComputeReal(void)
// Calcule une realisation dans l'espace reel
{
 if(lp_>0) cout<<"--- ComputeReal ---"<<endl;
 check_array_alloc();

 // On fait la FFT
 GEN3D_FFTW_EXECUTE(pb_);
 array_type = 1;
}

//-------------------------------------------------------------------
void GeneFluct3D::ReComputeFourier(void)
{
 if(lp_>0) cout<<"--- ReComputeFourier ---"<<endl;
 check_array_alloc();

 // On fait la FFT
 GEN3D_FFTW_EXECUTE(pf_);
 array_type = 2;

 // On corrige du pb de la normalisation de FFTW3
 complex<r_8> v((r_8)NRtot_,0.);
 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) /= v;
}

//-------------------------------------------------------------------
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
{
 if(lp_>0) cout<<"--- ComputeSpectrum ---"<<endl;
 check_array_alloc();

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

 // Attention a l'ordre
 for(long i=0;i<Nx_;i++) {
   double kx = AbsKx(i);  kx *= kx;
   for(long j=0;j<Ny_;j++) {
     double ky = AbsKy(j);  ky *= ky;
     for(long l=0;l<NCz_;l++) {
       double kz = AbsKz(l);
       double k = sqrt(kx+ky+kz*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;

 return 0;
}

int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr)
{
 if(lp_>0) cout<<"--- ComputeSpectrum2D ---"<<endl;
 check_array_alloc();

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

 // Attention a l'ordre
 for(long i=0;i<Nx_;i++) {
   double kx = AbsKx(i);  kx *= kx;
   for(long j=0;j<Ny_;j++) {
     double ky = AbsKy(j);  ky *= ky;
     double kt = sqrt(kx+ky);
     for(long l=0;l<NCz_;l++) {
       double kz = AbsKz(l);
       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;

 return 0;
}

//-------------------------------------------------------------------
int GeneFluct3D::ComputeSpectrum(HistoErr& herr,double sigma,bool pixcor)
// Compute spectrum from "T" and fill HistoErr "herr"
// AVEC la soustraction du niveau de bruit et la correction par filterpixel()
// Si on ne fait pas ca, alors on obtient un spectre non-isotrope!
//
// 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(lp_>0) cout<<"--- ComputeSpectrum: sigma="<<sigma<<endl;
 check_array_alloc();

 if(sigma<=0.) sigma = 0.;
 double sigma2 = sigma*sigma / (double)NRtot_;

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

 TVector<r_8> vfz(NCz_);
 if(pixcor)  // kz = l*Dkz_
   for(long l=0;l<NCz_;l++) {vfz(l)=pixelfilter(l*Dkz_ *Dz_/2); vfz(l)*=vfz(l);}

 // Attention a l'ordre
 for(long i=0;i<Nx_;i++) {
   double kx = AbsKx(i);
   double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
   kx *= kx; fx *= fx;
   for(long j=0;j<Ny_;j++) {
     double ky = AbsKy(j);
     double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
     ky *= ky; fy *= fy;
     for(long l=0;l<NCz_;l++) {
       double kz = AbsKz(l);
       double k = sqrt(kx+ky+kz*kz);
       double pk = MODULE2(T_(l,j,i)) - sigma2;
       double fz = (pixcor) ? vfz(l): 1.;
       double f = fx*fy*fz;
       if(f>0.) herr.Add(k,pk/f);
     }
   }
 }
 herr.ToVariance();
 for(int i=0;i<herr.NBins();i++) herr(i) += sigma2;

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

 return 0;
}

int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr,double sigma,bool pixcor)
// AVEC la soustraction du niveau de bruit et la correction par filterpixel()
{
 if(lp_>0) cout<<"--- ComputeSpectrum2D: sigma="<<sigma<<endl;
 check_array_alloc();

 if(sigma<=0.) sigma = 0.;
 double sigma2 = sigma*sigma / (double)NRtot_;

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

 TVector<r_8> vfz(NCz_);
 if(pixcor)  // kz = l*Dkz_
   for(long l=0;l<NCz_;l++) {vfz(l)=pixelfilter(l*Dkz_ *Dz_/2); vfz(l)*=vfz(l);}

 // Attention a l'ordre
 for(long i=0;i<Nx_;i++) {
   double kx = AbsKx(i);
   double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
   kx *= kx; fx *= fx;
   for(long j=0;j<Ny_;j++) {
     double ky = AbsKy(j);
     double fy = (pixcor) ? pixelfilter(ky*Dy_/2) : 1.;
     ky *= ky; fy *= fy;
     double kt = sqrt(kx+ky);
     for(long l=0;l<NCz_;l++) {
       double kz = AbsKz(l);
       double pk = MODULE2(T_(l,j,i)) - sigma2;
       double fz = (pixcor) ? vfz(l): 1.;
       double f = fx*fy*fz;
       if(f>0.) herr.Add(kt,kz,pk/f);
     }
   }
 }
 herr.ToVariance();
 for(int i=0;i<herr.NBinX();i++)
    for(int j=0;j<herr.NBinY();j++) herr(i,j) += sigma2;

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

 return 0;
}

//-------------------------------------------------------
int_8 GeneFluct3D::VarianceFrReal(double R,double& var)
// Recompute MASS variance in spherical top-hat (rayon=R)
// Par definition: SigmaR^2 = <(M-<M>)^2>/<M>^2
//                 ou M = masse dans sphere de rayon R
// --- ATTENTION: la variance calculee a une tres grande dispersion
//  (surtout si le volume du cube est petit). Pour verifier
//  que le sigmaR calcule par cette methode est en accord avec
//  le sigmaR en input, il faut faire plusieurs simulations (~100)
//  et regarder la moyenne des sigmaR reconstruits
{
  if(lp_>0) cout<<"--- VarianceFrReal R="<<R<<endl;
 check_array_alloc();

 long dnx = long(R/Dx_)+1; if(dnx<=0) dnx = 1;
 long dny = long(R/Dy_)+1; if(dny<=0) dny = 1;
 long dnz = long(R/Dz_)+1; if(dnz<=0) dnz = 1;
 if(lp_>0) cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;

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

 for(long i=dnx;i<Nx_-dnx;i+=2*dnx) {
   for(long j=dny;j<Ny_-dny;j+=2*dny) {
     for(long l=dnz;l<Nz_-dnz;l+=2*dnz) {
       double m=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);
             m += 1.+data_[ip];  // 1+drho/rho
             n++;
	   }
	 }
       }
       if(n>0) {sum += m; sum2 += m*m; nsum++; sn += n;}
       //cout<<i<<","<<j<<","<<l<<" n="<<n<<" m="<<m<<" sum="<<sum<<" sum2="<<sum2<<endl;
     }
   }
 }

 if(nsum<=1) {var=0.; return nsum;}
 sum /= nsum;
 sum2 = sum2/nsum - sum*sum;
 sn /= nsum;
 if(lp_>0) cout<<"...<n>="<<sn<<", nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
 var = sum2/(sum*sum);  // <dM^2>/<M>^2
 if(lp_>0) cout<<"...sigmaR^2 = <(M-<M>)^2>/<M>^2 = "<<var
               <<" -> sigmaR = "<<sqrt(var)<<endl;

 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++;
 }

 if(lp_>0) cout<<"--- NumberOfBad "<<nbad<<" px out of ]"<<vmin<<","<<vmax
               <<"[ i.e. frac="<<nbad/(double)NRtot_<<endl;
 return nbad;
}

int_8 GeneFluct3D::MinMax(double& xmin,double& xmax,double vmin,double vmax)
// Calcul des valeurs xmin et xmax dans le cube reel avec valeurs ]vmin,vmax[ extremites exclues
{
 bool tstval = (vmax>vmin)? true: false;
 if(lp_>0) {
   cout<<"--- MinMax";
   if(tstval) cout<<"  range=]"<<vmin<<","<<vmax<<"[";
   cout<<endl;
 }
 check_array_alloc();

 int_8 n = 0;
 xmin = xmax = data_[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 x = data_[ip];
   if(tstval && (x<=vmin || x>=vmax)) continue;
   if(x<xmin) xmin = x;
   if(x>xmax) xmax = x;
   n++;
 }

 if(lp_>0) cout<<"  n="<<n<<" min="<<xmin<<" max="<<xmax<<endl;

 return n;
}

int_8 GeneFluct3D::MeanSigma2(double& rm,double& rs2,double vmin,double vmax
                             ,bool useout,double vout)
// Calcul de mean,sigma2 dans le cube reel avec valeurs ]vmin,vmax[ extremites exclues
// useout = false: ne pas utiliser les pixels hors limites pour calculer mean,sigma2
//          true : utiliser les pixels hors limites pour calculer mean,sigma2
//                 en remplacant leurs valeurs par "vout"
{
 bool tstval = (vmax>vmin)? true: false;
 if(lp_>0) {
   cout<<"--- MeanSigma2";
   if(tstval) cout<<"  range=]"<<vmin<<","<<vmax<<"[";
   if(useout) cout<<", useout="<<useout<<" vout="<<vout;
   cout<<endl;
 }
 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(tstval) {
     if(v<=vmin || v>=vmax) {if(useout) v=vout; else continue;}
   }
   rm += v;
   rs2 += v*v;
   n++;
 }

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

 if(lp_>0) cout<<"  n="<<n<<" m="<<rm<<" s2="<<rs2<<" s="<<sqrt(fabs(rs2))<<endl;

 return n;
}

int_8 GeneFluct3D::SetToVal(double vmin, double vmax,double val0)
// set to "val0" if out of range ]vmin,vmax[ extremites exclues
// cad set to "val0" if in [vmin,vmax] -> 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++;}
 }

 if(lp_>0) cout<<"--- SetToVal "<<nbad<<" px set to="<<val0
               <<" because out of range=]"<<vmin<<","<<vmax<<"["<<endl;
 return nbad;
}

void GeneFluct3D::ScaleOffset(double scalecube,double offsetcube)
// Replace  "V"  by  "scalecube * ( V + offsetcube )"
{
 if(lp_>0) cout<<"--- ScaleCube scale="<<scalecube<<" offset="<<offsetcube<<endl;

 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] = scalecube * ( data_[ip] + offsetcube );
 }

 return;
}

//-------------------------------------------------------
void GeneFluct3D::TurnFluct2Mass(void)
// d_rho/rho -> Mass  (add one!)
{
 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++) {
   int_8 ip = IndexR(i,j,l);
   data_[ip] += 1.;
 }
}

double GeneFluct3D::TurnFluct2MeanNumber(double val_by_mpc3)
// ATTENTION: la gestion des pixels<0 proposee ici induit une perte de variance
// sur la carte, le spectre Pk reconstruit sera plus faible!
// L'effet sera d'autant plus grand que le nombre de pixels<0 sera grand.
{
 if(lp_>0) cout<<"--- TurnFluct2MeanNumber : "<<val_by_mpc3<<" quantity (gal or mass)/Mpc^3"<<endl;

 // First convert dRho/Rho into 1+dRho/Rho
 int_8  nball = 0; double sumall = 0., sumall2 = 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] += 1.;
   nball++; sumall += data_[ip]; sumall2 += data_[ip]*data_[ip];
 }
 if(nball>2) {
   sumall /= (double)nball;
   sumall2 = sumall2/(double)nball - sumall*sumall;
   if(lp_>0) cout<<"1+dRho/Rho: mean="<<sumall<<" variance="<<sumall2
                 <<" -> "<<sqrt(fabs(sumall2))<<endl;
 }

 // Find contribution for positive pixels
 int_8  nbpos = 0; double sumpos = 0. , sumpos2 = 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(data_[ip]>0.) {nbpos++; sumpos += v; sumpos2 += v*v;}
 }
 if(nbpos<1) {
   cout<<"TurnFluct2MeanNumber_Error: nbpos<1"<<endl;
   throw RangeCheckError("TurnFluct2MeanNumber_Error: nbpos<1");
 }
 sumpos2 = sumpos2/nball - sumpos*sumpos/(nball*nball);
 if(lp_>0)
   cout<<"1+dRho/Rho with v<0 set to zero: mean="<<sumpos/nball
       <<" variance="<<sumpos2<<" -> "<<sqrt(fabs(sumpos2))<<endl;
   cout<<"Sum of positive values: sumpos="<<sumpos
       <<" (n(v>0) = "<<nbpos<<" frac(v>0)="<<nbpos/(double)NRtot_<<")"<<endl;

 // - Mettre exactement val_by_mpc3*Vol galaxies (ou Msol) dans notre survey
 // - Uniquement dans les pixels de masse >0.
 // - Mise a zero des pixels <0
 double dn = val_by_mpc3 * Vol_ / sumpos;
 if(lp_>0) cout<<"...density move from "
               <<val_by_mpc3*dVol_<<" 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);
   if(data_[ip]<=0.) data_[ip] = 0.;
   else {
     data_[ip] *= dn;
     sum += data_[ip];
   }
 }

 if(lp_>0) cout<<"...quantity put into survey "<<sum<<" / "<<val_by_mpc3*Vol_<<endl;

 return sum;
}

double GeneFluct3D::ApplyPoisson(void)
// do NOT treate negative or nul mass  -> let it as it is
{
 if(lp_>0) cout<<"--- ApplyPoisson ---"<<endl;
 check_array_alloc();

 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.) {
     uint_8 dn = PoissonRand(v,10.);
     data_[ip] = (double)dn;
     sum += (double)dn;
   }
 }
 if(lp_>0) cout<<sum<<" galaxies put into survey"<<endl;

 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
{
 if(lp_>0) cout<<"--- TurnNGal2Mass ---"<<endl;
 check_array_alloc();

 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_>0) cout<<sum<<" MSol HI mass put into survey"<<endl;

 return sum;
}

double GeneFluct3D::TurnNGal2MassQuick(SchechterMassDist& schmdist)
// idem TurnNGal2Mass mais beaucoup plus rapide
{
 if(lp_>0) cout<<"--- TurnNGal2MassQuick ---"<<endl;
 check_array_alloc();

 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] = schmdist.TirMass(ngal);
     sum += data_[ip];
   }
 }
 if(lp_>0) cout<<sum<<" MSol HI mass put into survey"<<endl;

 return sum;
}

void GeneFluct3D::AddNoise2Real(double snoise,int type_evol)
// add noise to every pixels (meme les <=0 !)
// type_evol = 0 : pas d'evolution de la puissance du bruit
//             1 : evolution de la puissance du bruit avec la distance a l'observateur
//             2 : evolution de la puissance du bruit avec la distance du plan Z
//                 (tous les plans Z sont mis au meme redshift z de leur milieu)
{
 if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<type_evol<<endl;
 check_array_alloc();

 if(type_evol<0) type_evol = 0;
 if(type_evol>2) {
   const char *bla = "GeneFluct3D::AddNoise2Real_Error: bad type_evol value";
   cout<<bla<<endl; throw ParmError(bla);
 }

 vector<double> correction;
 InterpFunc *intercor = NULL;

 if(type_evol>0) {
   // Sigma_Noise(en mass) :
   //      Slim ~ 1/sqrt(DNu) * sqrt(nlobe)   en W/m^2Hz
   //      Flim ~ sqrt(DNu) * sqrt(nlobe)   en W/m^2
   //      Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe)  en Msol
   //                                    nlobe ~ 1/Dtrcom^2
   //      Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
   if(cosmo_ == NULL || redsh_ref_<0. || loscom2zred_.size()<1) {
     const char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
     cout<<bla<<endl; throw ParmError(bla);
   }
   InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
   long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
   for(long i=0;i<nsz;i++) {
     double d = interpinv.X(i);
     double zred = interpinv(d);
     double dtrc = cosmo_->Dtrcom(zred);  // pour variation angle solide
     double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
     double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
     double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
     double corr = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
     if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
       cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
           <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
     correction.push_back(corr);
   }
   intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
 } 

 double corrlim[2] = {1.,1.};
 for(long i=0;i<Nx_;i++) {
   double dx2 = DXcom(i); dx2 *= dx2;
   for(long j=0;j<Ny_;j++) {
     double dy2 = DYcom(j); dy2 *= dy2;
     for(long l=0;l<Nz_;l++) {
       double corr = 1.;
       if(type_evol>0) {
         double dz = DZcom(l);
         if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
           else dz = fabs(dz); // tous les plans Z au meme redshift
         corr = (*intercor)(dz);
         if(corr<corrlim[0]) corrlim[0]=corr; else if(corr>corrlim[1]) corrlim[1]=corr;
       }
       int_8 ip = IndexR(i,j,l);
       data_[ip] += snoise*corr*NorRand();
     }
   }
 }
 if(type_evol>0)
   cout<<"correction factor range: ["<<corrlim[0]<<","<<corrlim[1]<<"]"<<endl;

 if(intercor!=NULL) delete intercor;
}

}  // Fin namespace SOPHYA




/*********************************************************************
void GeneFluct3D::AddAGN(double lfjy,double lsigma,double powlaw)
// Add AGN flux into simulation:
// --- Procedure:
// 1. lancer "cmvdefsurv" avec les parametres du survey
//        (au redshift de reference du survey)
//    et recuperer l'angle solide "angsol sr" du pixel elementaire
//    au centre du cube.
// 2. lancer "cmvtstagn" pour cet angle solide -> cmvtstagn.ppf
// 3. regarder l'histo "hlfang" et en deduire un equivalent gaussienne
//    cad une moyenne <log10(S)> et un sigma "sig"
//    Attention: la distribution n'est pas gaussienne les "mean,sigma"
//               de l'histo ne sont pas vraiment ce que l'on veut
// --- Limitations actuelle du code:
// . les AGN sont supposes evoluer avec la meme loi de puissance pour tout theta,phi
// . le flux des AGN est mis dans une colonne Oz (indice k) et pas sur la ligne de visee
// . la distribution est approximee a une gaussienne
// ... C'est une approximation pour un observateur loin du centre du cube
//     et pour un cube peu epais / distance observateur 
// --- Parametres de la routine:
// llfy : c'est le <log10(S)> du flux depose par les AGN
//        dans l'angle solide du pixel elementaire de reference du cube
// lsigma : c'est le sigma de la distribution des log10(S)
// powlaw : c'est la pente de la distribution cad que le flux "lmsol"
//          et considere comme le flux a 1.4GHz et qu'on suppose une loi
//             F(nu) = (1.4GHz/nu)^powlaw * F(1.4GHz)
// - Comme on est en echelle log10():
//   on tire log10(Msol) + X
//   ou X est une realisation sur une gaussienne de variance "sig^2"
//   La masse realisee est donc: Msol*10^X
// - Pas de probleme de pixel negatif car on a une multiplication!
{
  if(lp_>0) cout<<"--- AddAGN: <log10(S Jy)> = "<<lfjy<<" , sigma = "<<lsigma<<endl;
  check_array_alloc();

 if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
   char *bla = "GeneFluct3D::AddAGN_Error: set Observator and Cosmology first";
   cout<<bla<<endl; throw ParmError(bla);
 }

 // Le flux des AGN en Jy et en mass solaire
 double fagnref = pow(10.,lfjy)*(dnu_ref_*1.e9); // Jy.Hz = W/m^2
 double magnref = FluxHI2Msol(fagnref*Jansky2Watt_cst,dlum_ref_); // Msol
 if(lp_>0)
   cout<<"Au pixel de ref: fagnref="<<fagnref
       <<" Jy.Hz (a 1.4GHz), magnref="<<magnref<<" Msol"<<endl;

 if(powlaw!=0.) {
   // F(nu) = F(1.4GHz)*(nu GHz/1.4 Ghz)^p = F(1.4GHz)*(1/(1+z))^p , car nu = 1.4 GHz/(1+z)
   magnref *= pow(1/(1.+redsh_ref_),powlaw);
   if(lp_>0) cout<<" powlaw="<<powlaw<<"  -> change magnref to "<<magnref<<" Msol"<<endl;
 }

 // Les infos en fonction de l'indice "l" selon Oz
 vector<double> correction;
 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
 long nzmod = ((Nz_>10)?Nz_/10:1);
 for(long l=0;l<Nz_;l++) {
   double z = fabs(DZcom(l));
   double zred = interpinv(z);
   double dtrc = cosmo_->Dtrcom(zred);  // pour variation angle solide
   double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
   double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
   double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
   // on a: Mass ~ DNu * Dlum^2 / Dtrcom^2
   double corr = dnu/dnu_ref_*pow(dtrc_ref_/dtrc*dlum/dlum_ref_,2.);
   // F(nu) = F(1.4GHz)*(nu GHz/1.4 Ghz)^p = F(1.4GHz)*(1/(1+z))^p , car nu = 1.4 GHz/(1+z)
   if(powlaw!=0.) corr *= pow((1.+redsh_ref_)/(1.+zred),powlaw);
   correction.push_back(corr);
   if(lp_>0 && (l==0 || l==Nz_-1 || l%nzmod==0)) {
     cout<<"l="<<l<<" z="<<z<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
         <<" dtrc="<<dtrc<<" dlum="<<dlum
         <<" -> cor="<<corr<<endl;
   }
 }

 double sum=0., sum2=0., nsum=0.;
 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
   double a = lsigma*NorRand();
   a = magnref*pow(10.,a);
   // On met le meme tirage le long de Oz (indice k)
   for(long l=0;l<Nz_;l++) {
     int_8 ip = IndexR(i,j,l);
     data_[ip] += a*correction[l];
   }
   sum += a; sum2 += a*a; nsum += 1.;
 }

 if(lp_>0 && nsum>1.) {
   sum /= nsum;
   sum2 = sum2/nsum - sum*sum;
   cout<<"...Mean mass="<<sum<<" Msol , s^2="<<sum2<<" s="<<sqrt(fabs(sum2))<<endl;
 }
 
}
*********************************************************************/