source: Sophya/trunk/Cosmo/SimLSS/pkspectrum.h@ 3805

#ifndef PKSPECTRUM_SEEN
#define PKSPECTRUM_SEEN

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

namespace SOPHYA {

//-----------------------------------------------------------------------------------
class InitialSpectrum : public GenericFunc {
public:
  InitialSpectrum(void) {};
  InitialSpectrum(InitialSpectrum& pkinf) {};
  virtual ~InitialSpectrum(void) {};
};

//-----------------------------------------------------------------------------------
class InitialPowerLaw : public InitialSpectrum {
public:
  InitialPowerLaw(double n,double a=1.);
  InitialPowerLaw(InitialPowerLaw& pkinf);
  virtual ~InitialPowerLaw(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 TransfertFunction : public GenericFunc {
public:
  TransfertFunction(void) {};
  virtual ~TransfertFunction(void) {};
};

//-----------------------------------------------------------------------------------
class TransfertEisenstein : public TransfertFunction {
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 TransfertFunction {
public:
  TransfertTabulate(void);
  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,double h100,double OmegaCDM0,double OmegaBaryon0);
  int ReadCAMB(string filename, double h100=0.71);
protected:
  double kmin_,kmax_;
  int interptyp_;
  vector<double> k_, tf_;
};


//-----------------------------------------------------------------------------------
class GrowthFactor : public GenericFunc {
public:
  GrowthFactor(void) {};
  virtual ~GrowthFactor(void) {};
  virtual double DsDz(double z, double);
};

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


//-----------------------------------------------------------------------------------
class PkSpectrum : public GenericFunc {
public:
  // typsec = PK : compute Pk(k)
  //        = DELTA : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
  typedef enum {PK=0, DELTA=1} ReturnSpectrum;

  PkSpectrum(void);
  PkSpectrum(PkSpectrum& pk);
  virtual ~PkSpectrum(void) {};
  virtual void   SetZ(double z) {zref_ = z;}
  virtual double GetZ(void) {return zref_;}
  virtual void SetScale(double scale=1.) {scale_ = scale;}
  virtual double GetScale(void) {return scale_;}
  virtual void SetTypSpec(ReturnSpectrum typspec=PK) {typspec_ = typspec;}
  virtual ReturnSpectrum GetTypSpec(void) {return typspec_ ;}
protected:
  double zref_, scale_;
  ReturnSpectrum typspec_;
};

//-----------------------------------------------------------------------------------
class PkSpecCalc : public PkSpectrum {
public:
  PkSpecCalc(InitialSpectrum& pkinf,TransfertFunction& tf,GrowthFactor& d1,double zref=0.);
  PkSpecCalc(PkSpecCalc& pkz);
  virtual ~PkSpecCalc(void);
  virtual double operator() (double k);
  virtual double operator() (double k,double z);
  InitialSpectrum& GetPkIni(void) {return pkinf_;}
  TransfertFunction& GetTransfert(void) {return tf_;}
  GrowthFactor& GetGrowthFactor(void) {return d1_;}
protected:
  InitialSpectrum& pkinf_;
  TransfertFunction& tf_;
  GrowthFactor& d1_;
};

//-----------------------------------------------------------------------------------
class PkTabulate : public PkSpectrum {
public:
  PkTabulate(void);
  PkTabulate(PkTabulate& pkz);
  virtual ~PkTabulate(void);
  virtual double operator() (double k);
  virtual double operator() (double k,double z);
  int NPoints(void) {return k_.size();}
  void SetInterpTyp(int typ=0);
  int ReadCAMB(string filename, double h100=0.71, double zreftab = 0.);
protected:
  double kmin_,kmax_;
  int interptyp_;
  vector<double> k_, pk_;
};

//-----------------------------------------------------------------------------------
class PkEisenstein : public PkSpectrum {
public:
  PkEisenstein(InitialPowerLaw& pkinf,TransfertEisenstein& tf,GrowthEisenstein& d1,double zref=0.);
  PkEisenstein(PkEisenstein& pkz);
  virtual ~PkEisenstein(void);
  virtual double operator() (double k);
  virtual double operator() (double k,double z);
  InitialPowerLaw& GetPkIni(void) {return pkinf_;}
  TransfertEisenstein& GetTransfert(void) {return tf_;}
  GrowthEisenstein& GetGrowthFactor(void) {return d1_;}
protected:
  InitialPowerLaw& pkinf_;
  TransfertEisenstein& tf_;
  GrowthEisenstein& d1_;
};

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

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

  VarianceSpectrum(GenericFunc& pk,double R,TypeFilter typfilter);
  VarianceSpectrum(VarianceSpectrum& vpk);
  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