#ifndef GENEFLUCT3D_SEEN
#define GENEFLUCT3D_SEEN

#include "machdefs.h"
#include "genericfunc.h"
#include "tarray.h"
#include "hisprof.h"
#include "perandom.h"

#include <vector>

#include "pkspectrum.h"

#include "FFTW/fftw3.h"

namespace SOPHYA {

//-----------------------------------------------------------------------------------
class GeneFluct3D {
public:
  GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz);
  virtual ~GeneFluct3D(void);

  void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
  void SetSize(long nx,long ny,long nz,double dx,double dy,double dz);  // Mpc

  inline void SetZ(double z) {pkz_.SetZ(z);}
  double GetZ(void) {return pkz_.GetZ();}
  double GetVol(void) {return Vol_;}
  double GetDVol(void) {return dVol_;}
  size_t NPix(void) {return NRtot_;}
  double GetKmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_+Knyqz_*Knyqz_);}
  vector<r_8> GetKinc(void)
    {vector<r_8> dk; dk.push_back(Dkx_); dk.push_back(Dky_); dk.push_back(Dkz_); return dk;}
  vector<r_8> GetKnyq(void)
    {vector<r_8> kn; kn.push_back(Knyqx_); kn.push_back(Knyqy_); kn.push_back(Knyqz_); return kn;}
  vector<r_8> GetDinc(void)
    {vector<r_8> d; d.push_back(Dx_); d.push_back(Dy_); d.push_back(Dz_); return d;}
  vector<long> GetNpix(void)
    {vector<long> d; d.push_back(Nx_); d.push_back(Ny_); d.push_back(Nz_); return d;}

  void ComputeFourier0(void);
  void ComputeFourier(void);
  void FilterByPixel(void);
  int  ComputeSpectrum(HProf& hp);
  void ComputeReal(void);
  void ReComputeFourier(void);

  size_t VarianceFrReal(double R,double& var);
  size_t MeanSigma2(double& rm,double& rs2,double vmin=-1.e+150,double vmax=1.e+150);
  size_t NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
  size_t SetToVal(double vmin, double vmax,double val0=0.);

  void TurnFluct2Mass(void);
  double TurnMass2MeanNumber(double n_by_mpc3);
  double ApplyPoisson(void);
  double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
  double TurnMass2Flux(void);
  void AddNoise2Real(double snoise);

  void Print(void);

protected:
  long manage_coefficients(void);
  double compute_power_carte(void);
  inline double pixelfilter(double x)
    {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}


  long Nx_,Ny_,Nz_;
  sa_size_t SzK_[3];
  long NCz_,NTz_;
  size_t NRtot_;
  double Dx_,Dy_,Dz_;
  double dVol_, Vol_;
  double Dkx_,Dky_,Dkz_;
  double Knyqx_,Knyqy_,Knyqz_;
  double Dk3_;

  fftw_plan pf_,pb_;
  unsigned short nthread_;

  //   = 0 : T_ vide
  //   = 1 : T_ contient une realisation du spectre (espace fourier)
  //   = 2 : T_ contient une realisation dans l'espace reel
  //   = 3 : T_ contient la TF d'une realisation dans l'espace reel (d_rho/rho)
  //   = 4 : T_ contient (1 + d_rho/rho)
  //   = ... etc ...
  unsigned short Tcontent_;

  TArray< complex<r_8 > >& T_;
  PkSpectrumZ& pkz_;
};

} // Fin du namespace

#endif