#ifndef PKSPECTRUM_SEEN
#define PKSPECTRUM_SEEN

#include "machdefs.h"
#include "genericfunc.h"

namespace SOPHYA {

//-----------------------------------------------------------------------------------
class InitialSpectrum : public GenericFunc {
public:
  InitialSpectrum(double n,double a=1.);
  InitialSpectrum(InitialSpectrum& pkinf);
  virtual ~InitialSpectrum(void);
  virtual double operator() (double k) {return A_ * pow(k,n_);}
  void SetNorm(double a) {A_ = a;}
  void SetSlope(double n) {n_ = n;}
protected:
  double n_, A_;
};

//-----------------------------------------------------------------------------------
class TransfertEisenstein : public GenericFunc {
public:

  typedef enum{ALL=0, CDM=1, BARYON=2} ReturnPart;

  TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon=false,int lp=0);
  TransfertEisenstein(TransfertEisenstein& tf);
  virtual ~TransfertEisenstein(void);
  bool SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0);
  virtual double operator() (double k);
  double KPeak(void);
  void SetNoOscEnv(unsigned short nooscenv=0);
  void SetReturnPart(ReturnPart retpart=ALL);
  void SetPrintLevel(int lp=0);
protected:
  int lp_;
  double O0_,Oc_,Ob_,h100_,tcmb_;
  double th2p7_;
  double zeq_,keq_,zd_,Req_,Rd_,s_,ksilk_,alphac_,betac_,bnode_,alphab_,betab_;
  double alphag_;
  double sfit_,kpeak_;

  bool nobaryon_;
  unsigned short nooscenv_;
  ReturnPart retpart_;

  double T0tild(double k,double alphac,double betac);
  void Init_(void);
  void zero_(void);
};

//-----------------------------------------------------------------------------------
class TransfertTabulate : public GenericFunc {
public:
  TransfertTabulate(double h100,double OmegaCDM0,double OmegaBaryon0);
  TransfertTabulate(TransfertTabulate& tf);
  virtual ~TransfertTabulate(void);
  virtual double operator() (double k);
  int NPoints(void) {return k_.size();}
  void SetInterpTyp(int typ=0);
  int ReadCMBFast(string filename);
protected:
  double Oc_,Ob_,h100_;
  double kmin_,kmax_;
  int interptyp_;
  vector<double> k_, tf_;
};


//-----------------------------------------------------------------------------------
class GrowthFactor : public GenericFunc {
public:
  GrowthFactor(double OmegaMatter0,double OmegaLambda0);
  GrowthFactor(GrowthFactor& d1);
  virtual ~GrowthFactor(void);
  virtual double operator() (double z);
  void SetParTo(double OmegaMatter0,double OmegaLambda0);
  bool SetParTo(double OmegaMatter0);
protected:
  double O0_,Ol_;
};


//-----------------------------------------------------------------------------------
class PkSpectrum0 : public GenericFunc {
public:
  PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf);
  PkSpectrum0(PkSpectrum0& pk0);
  virtual ~PkSpectrum0(void);
  virtual double operator() (double z);
  InitialSpectrum& GetPkIni(void) {return pkinf_;}
  TransfertEisenstein& GetTransfert(void) {return tf_;}
protected:
  InitialSpectrum& pkinf_;
  TransfertEisenstein& tf_;
};

//-----------------------------------------------------------------------------------
class PkSpectrumZ : public GenericFunc {
public:
  typedef enum {PK=0, DELTA=1} ReturnSpectrum;
  PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref=0.);
  PkSpectrumZ(PkSpectrumZ& pkz);
  virtual ~PkSpectrumZ(void);
  virtual double operator() (double k);
  virtual double operator() (double k,double z);
  void   SetZ(double z) {zref_ = z;}
  double GetZ(void) {return zref_;}
  void SetTypSpec(ReturnSpectrum typspec=PK);
  void SetScale(double scale=1.) {scale_=scale;}
  double GetScale(void) {return scale_;}
  PkSpectrum0& GetPk0(void) {return pk0_;}
  GrowthFactor& GetGrowthFactor(void) {return d1_;}
protected:
  PkSpectrum0& pk0_;
  GrowthFactor& d1_;
  double zref_, scale_;
  ReturnSpectrum typspec_;
};

//-----------------------------------------------------------------------------------
class VarianceSpectrum : public GenericFunc {
public:

  typedef enum {TOPHAT=0, GAUSSIAN=1, NOFILTER=2} TypeFilter;

  VarianceSpectrum(GenericFunc& pk,double R,TypeFilter typfilter);
  VarianceSpectrum(VarianceSpectrum& pkinf);
  virtual ~VarianceSpectrum(void);

  void SetRadius(double R);
  void SetFilter(TypeFilter typfilter=TOPHAT);
  void SetInteg(double dperc=0.1,double dlogkinc=-1.,double dlogkmax=-1.,unsigned short glorder=4);

  double Variance(double kmin,double kmax);

  // ATTENTION: La fonction a integrer est : f(k)dk = k^3*Pk(k)/(2Pi^2) *filter2(k*R) *dk/k
  virtual double operator() (double k) {return k*k*pk_(k)*Filter2(k*R_)/(2.*M_PI*M_PI);}
  double Filter2(double x);

  // Aide a l'integration
  double FindMaximum(double kmin,double kmax,double eps=1.e-3);
  int FindLimits(double high,double &kmin,double &kmax,double eps=1.e-3);

protected:

  GenericFunc& pk_;
  TypeFilter typfilter_;
  double R_;

  double dperc_,dlogkinc_,dlogkmax_;
  unsigned short glorder_;

};

} // Fin du namespace SOPHYA

#endif