#ifndef GENEFLUCT3D_SEEN
#define GENEFLUCT3D_SEEN

#include "machdefs.h"
#include "genericfunc.h"
#include "tarray.h"
#include "histerr.h"
#include "hist2err.h"
#include "perandom.h"

#include "FFTW/fftw3.h"
#include "FitsIO/fitsio.h"

#include <vector>

#include "pkspectrum.h"


namespace SOPHYA {

//-----------------------------------------------------------------------------------
class GeneFluct3D {
public:
  GeneFluct3D(TArray< complex<r_8 > >& T);
  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
  void SetObservator(double redshref=0.,double kredshref=-1.);

  TArray< complex<r_8> >& GetComplexArray(void) {return T_;}
  fftw_complex* GetComplexPointer(void) {return fdata_;}
  TArray<r_8>& GetRealArray(void) {return R_;}
  r_8* GetRealPointer(void) {return data_;}

  // Pour adressage data_[ip]
  inline int_8 IndexR(long i,long j,long k) {return (int_8)(k+NTz_*(j+Ny_*i));}
  // Pour adressage fdata_[ip][0-1]
  inline int_8 IndexC(long i,long j,long k) {return (int_8)(k+NCz_*(j+Ny_*i));}
  // - On peut aussi adresser:
  // TArray< complex<r_8> >& pk = gf3d.GetComplexArray(); pk(i,j,k) = ...;
  // TArray<r_8>& rgen = gf3d.GetRealArray(); rgen(i,j,k) = ...;

  vector<long> GetNpix(void) {return N_;}
  int_8 NPix(void) {return NRtot_;}

  vector<r_8> GetDinc(void) {return D_;}
  double GetDVol(void) {return dVol_;}
  double GetVol(void) {return Vol_;}

  vector<r_8> GetKinc(void) {return Dk_;}
  vector<r_8> GetKnyq(void) {return Knyq_;}
  double GetKmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_+Knyqz_*Knyqz_);}
  double GetKTmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_);}
  double GetKincMin(void)
    {vector<r_8>::const_iterator it = min_element(Dk_.begin(), Dk_.end()); return *it;}
  double GetKTincMin(void) {return min(Dk_[0],Dk_[1]);}

  void ComputeFourier0(GenericFunc& pk_at_z);
  void ComputeFourier(GenericFunc& pk_at_z);
  void FilterByPixel(void);

  void ComputeReal(void);

  void ReComputeFourier(void);

  int  ComputeSpectrum(HistoErr& herr);
  int  ComputeSpectrum2D(Histo2DErr& herr);

  int_8 VarianceFrReal(double R,double& var);
  int_8 MeanSigma2(double& rm,double& rs2,double vmin=-1.e+150,double vmax=1.e+150);
  int_8 NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
  int_8 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 WriteFits(string cfname,int bitpix=FLOAT_IMG);
  void ReadFits(string cfname);

  void WritePPF(string cfname,bool write_real=true);
  void ReadPPF(string cfname);

  void SetPrtLevel(int lp=0) {lp_ = lp;}
  void Print(void);

protected:
  void setsize(long nx,long ny,long nz,double dx,double dy,double dz);
  void setalloc(void);
  void setpointers(bool from_real);
  void init_fftw(void);
  long manage_coefficients(void);
  double compute_power_carte(void);
  void check_array_alloc(void);
  inline double pixelfilter(double x)
    {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}

  // valeurs dans l'espace reel
  long Nx_,Ny_,Nz_;  vector<long> N_;
  long NCz_,NTz_;
  int_8 NRtot_;

  double Dx_,Dy_,Dz_;  vector<double> D_;

  // valeurs dans l'espace des K
  double Dkx_,Dky_,Dkz_;  vector<double> Dk_;
  double Knyqx_,Knyqy_,Knyqz_;  vector<double> Knyq_;
  double Dk3_;
  double dVol_, Vol_;

  // la gestion de la FFT
  fftw_plan pf_,pb_;
  unsigned short nthread_;
  int lp_;

  // le stockage du Cube de donnees et les pointeurs
  bool array_allocated_;  // true if array has been allocated
  TArray< complex<r_8> >& T_;
  fftw_complex *fdata_;
  TArray<r_8> R_;
  double *data_;

  // l'observateur
  double redshref_,kredshref_;
};

} // Fin du namespace

#endif