source: Sophya/trunk/Cosmo/SimLSS/pkspectrum.cc@ 3196

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

#include "pexceptions.h"

#include "constcosmo.h"
#include "geneutils.h"
#include "pkspectrum.h"


///////////////////////////////////////////////////////////
//******************** InitialSpectrum ******************//
///////////////////////////////////////////////////////////

InitialSpectrum::InitialSpectrum(double n,double a)
  : n_(n), A_(a)
{
}

InitialSpectrum::InitialSpectrum(InitialSpectrum& pkinf)
  : n_(pkinf.n_), A_(pkinf.A_)
{
}

InitialSpectrum::~InitialSpectrum(void)
{
}

void InitialSpectrum::SetNorm(double a)
{
  A_ = a;
}

void InitialSpectrum::SetSlope(double n)
{
  n_ = n;
}


///////////////////////////////////////////////////////////
//****************** TransfertEisenstein ****************//
///////////////////////////////////////////////////////////

// From Eisenstein & Hu ApJ 496:605-614 1998 April 1
TransfertEisenstein::TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon)
  : Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb)
  , nobaryon_(nobaryon) , nooscenv_(0)
{
  if(nobaryon_ == false && Ob_>0.) Init_With_Baryons();
  else Init_Without_Baryon();
}

TransfertEisenstein::TransfertEisenstein(TransfertEisenstein& tf)
  : Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , tcmb_(tf.tcmb_)
  , nobaryon_(tf.nobaryon_) , nooscenv_(tf.nooscenv_)
{
  if(nobaryon_ == false && Ob_>0.) Init_With_Baryons();
  else Init_Without_Baryon();
}

void TransfertEisenstein::Init_With_Baryons(void)
{
  int lp = 1;

  nobaryon_ = false;
  O0_ = Oc_ + Ob_;
    if(lp) cout<<"Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;

  double H0 = 100. * h_, h2 = h_*h_;

  th2p7_ = tcmb_/2.7;
  double th2p7P4 = th2p7_*th2p7_*th2p7_*th2p7_;

  // Formule 2 p 606
  zeq_ = 2.50e4 * O0_ * h2 / th2p7P4;
    if(lp) cout<<"zeq = "<<zeq_<<endl;

  // Formule 3 p 607
  // (attention ici C=1 : H0 -> H0/C si on utilise la premiere formule)
  //  keq_ = sqrt(2.*O0_*H0*H0*zeq_) / SpeedOfLight_Cst;
  keq_ = 7.46e-2 * O0_ * h2 / (th2p7_*th2p7_);
    if(lp) cout<<"keq = "<<keq_<<" Mpc^-1"<<endl;

  // Formule 4 p 607
  double b1_eq4 = 0.313*pow(O0_*h2,-0.419)*(1. + 0.607*pow(O0_*h2,0.674));
  double b2_eq4 = 0.238*pow(O0_*h2,0.223);
  double zd_ = 1291. * pow(O0_*h2,0.251) / (1.+0.659* pow(O0_*h2,0.828))
             * (1. + b1_eq4*pow(Ob_*h2,b2_eq4));
    if(lp) cout<<"zd = "<<zd_<<endl;

  // Formule 5 page 607
  Req_ = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zeq_);
  Rd_  = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zd_);
    if(lp) cout<<"Req = "<<Req_<<" Rd = "<<Rd_<<endl;

  // Formule 6 p 607
  s_ = 2./(3.*keq_) * sqrt(6./Req_)
      * log( (sqrt(1.+Rd_) + sqrt(Rd_+Req_)) / (1.+sqrt(Req_)) );
    if(lp) cout<<"s = "<<s_<<" Mpc"<<endl;

  // Formule 7 page 607
  ksilk_ = 1.6*pow(Ob_*h2,0.52)*pow(O0_*h2,0.73) * (1. + pow(10.4*O0_*h2,-0.95));
    if(lp) cout<<"ksilk = "<<ksilk_<<" Mpc^-1"<<endl;

  // Formules 10 page 608
  double a1 = pow(46.9*O0_*h2,0.670) * (1. + pow(32.1*O0_*h2,-0.532));
  double a2 = pow(12.0*O0_*h2,0.424) * (1. + pow(45.0*O0_*h2,-0.582));
  alphac_ = pow(a1,-Ob_/O0_) * pow(a2,-pow(Ob_/O0_,3.));
  double b1 = 0.944 / (1. + pow(458.*O0_*h2,-0.708));
  double b2 = pow(0.395*O0_*h2,-0.0266);
  betac_ = 1 / ( 1. + b1*(pow(Oc_/O0_,b2) - 1.) );

  // Formule 23 page 610
  bnode_ = 8.41 * pow(O0_*h2,0.435);
    if(lp) cout<<"bnode = "<<bnode_<<endl;

  // Formule 14 page 608
  double y = (1.+zeq_)/(1.+zd_);
  double s1py = sqrt(1.+y);
  double Gy = y*( -6.*s1py + (2.+3.*y)*log((s1py+1.)/(s1py-1.)) );
  alphab_ = 2.07*keq_*s_*pow(1.+Rd_,-3./4.)*Gy;

  // Formule 24 page 610
  betab_ = 0.5 + Ob_/O0_
         + (3.-2.*Ob_/O0_) * sqrt(pow(17.2*O0_*h2,2.) + 1.);

  // Formule 31 page 612
  alphag_ = 1.
          - 0.328*log(431.*O0_*h2)*Ob_/O0_
          + 0.38*log(22.3*O0_*h2)*pow(Ob_/O0_,2.);

 // --- Pour la positon du premier pic acoustique
 // Formule 26 page 611
 double s_peak = 44.5*log(9.83/(O0_*h2)) / sqrt(1.+10.*pow(Ob_*h2,3./4.));  // Mpc
 // Formule 25 page 611
 kpeak_ = 5*M_PI/(2.*s_peak) * (1.+0.217*O0_*h2);  // 1/Mpc
   if(lp) cout<<"for 1er peak: s="<<s_peak<<" kpeak="<<kpeak_<<endl;

  return;
}

void TransfertEisenstein::Init_Without_Baryon(void)
{
  int lp = 1;

  nobaryon_ = true;
  O0_ = Oc_;
  if(lp) cout<<"Omatter="<<O0_<<" Ocdm="<<Oc_<<endl;
  th2p7_ = tcmb_/2.7;
}

TransfertEisenstein::~TransfertEisenstein(void)
{
}

void TransfertEisenstein::SetNoOscEnv(unsigned short nooscenv)
// To obtain the non-oscillatory part of the spectrum
// nooscenv = 0 : use the baryon oscillatory part of transfert function
// nooscenv = 1 : use approx. paragraph 3.3 p610 (middle of right column)
// nooscenv = 2 : use formulae 29+30+31 page 612
{
  if(nooscenv>2) {
    cout<<"TransfertEisenstein::SetNoOscEnv: Error bad nooscenv value: "<<nooscenv<<endl;
    throw ParmError("TransfertEisenstein::SetNoOscEnv: Error bad nooscenv value");
  }
  nooscenv_ = nooscenv;
}

double TransfertEisenstein::T0tild(double k,double alphac,double betac)
{
  // Formule 10 p 608
  //double q = k*th2p7_*th2p7_/(O0_*h_*h_);
  double q = k/(13.41*keq_);
  // Formule 20 p 610
  double C = (14.2/alphac) + 386./(1.+69.9*pow(q,1.08));
  // Formule 19 p 610
  double x = log(M_E+1.8*betac*q);
  return x / (x + C*q*q);
}

double TransfertEisenstein::operator() (double k)
{

  // --- Pour zero baryon
  //  OU Pour function lissee sans oscillation baryon
  if(nobaryon_  || nooscenv_ == 2) {
    double gamma = O0_*h_;
    // Formule 30 page 612 (pour fct lissee)
    if( nobaryon_==false && nooscenv_ == 2 )
      gamma = O0_*h_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*s_,4.)));
    // Formule 28 page 612
    double q = k / h_ * th2p7_*th2p7_/gamma;  // Mpc^-1
    // Formules 29 page 612
    double l0 = log(2.*M_E + 1.8*q);
    double c0 = 14.2 + 731./(1.+62.5*q);
    return l0 / (l0 + c0*q*q);
  }

  // --- CDM
  double f = 1. / (1. + pow(k*s_/5.4,4.));
  double Tc = f*T0tild(k,1.,betac_) + (1.-f)*T0tild(k,alphac_,betac_);

  // --- Baryons
  // Formule 22 page 610
  double stk, ksbnode = k*s_/bnode_;
  if(ksbnode<0.001) stk =s_ * ksbnode;
     else   stk = s_ / pow(1. + pow(1./ksbnode,3.), 1/.3);
  // Formule 21 page 610
  double j0kst = 0.;
  if(nooscenv_ == 1) j0kst = pow(1.+pow(k*stk,4.) , -1./4.); //lissee sans oscillation baryon
  else {
    double x = k*stk;
    if(x<0.01) j0kst = 1. - x*x/6.*(1.-x*x/20.);
      else j0kst = sin(x)/x;
    ////if(k>0.038&&k<0.07) cout<<"k="<<k<<" stk="<<stk<<" x="<<x<<endl;
  }
  double Tb = T0tild(k,1.,1.) / (1. + pow(k*s_/5.2,2.));
  Tb += alphab_/(1.+pow(betab_/(k*s_),3.)) * exp(-pow(k/ksilk_,1/4.));
  Tb *= j0kst;

  // --- Total
  double T = (Ob_/O0_)*Tb + (Oc_/O0_)*Tc;

  return T;
}

double TransfertEisenstein::KPeak(void)
// Position du premier pic acoustic
{
 if(nobaryon_) return -1.;
 return kpeak_;
}


///////////////////////////////////////////////////////////
//********************* GrowthFactor ********************//
///////////////////////////////////////////////////////////

// From Eisenstein & Hu ApJ 496:605-614 1998 April 1
// Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
  : O0_(OmegaMatter0) , Ol_(OmegaLambda0) , Ok_(1.-OmegaMatter0-OmegaLambda0)
{
  if(OmegaMatter0==0.) {
    cout<<"GrowthFactor::GrowthFactor: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
    throw ParmError("GrowthFactor::GrowthFactor: Error badOmegaMatter0  value");
  }
  norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
  norm_ = (*this)(0.);
  cout<<"GrowthFactor::GrowthFactor : norm="<<norm_<<endl;
}

GrowthFactor::GrowthFactor(GrowthFactor& d1)
  : O0_(d1.O0_) , Ol_(d1.Ol_) , Ok_(d1.Ok_) , norm_(d1.norm_)
{
}

GrowthFactor::~GrowthFactor(void)
{
}

double GrowthFactor::operator() (double z)
// see Formulae A4 + A5 + A6 page 614
{
 z += 1.;
 double z2 = z*z, z3 = z2*z;
 double den = Ol_ + Ok_*z2 + O0_*z3;
 double o0z = O0_ *z3 / den;
 double olz = Ol_ / den;

 // 4./7. = 0.571429
 double D1z = pow(o0z,0.571429) - olz + (1.+o0z/2.)*(1.+olz/70.);
 D1z = 2.5*o0z / z / D1z;

 return D1z / norm_;
}


///////////////////////////////////////////////////////////
//************** PkSpectrum0 et PkSpectrumZ *************//
///////////////////////////////////////////////////////////

PkSpectrum0::PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf)
  : pkinf_(pkinf) , tf_(tf)
{
}

PkSpectrum0::PkSpectrum0(PkSpectrum0& pk0)
  : pkinf_(pk0.pkinf_) , tf_(pk0.tf_)
{
}

PkSpectrum0::~PkSpectrum0(void)
{
}

double PkSpectrum0::operator() (double k)
{
  double tf = tf_(k);
  double pkinf = pkinf_(k);
  return pkinf *tf*tf;
}

//------------------------------------
PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
  : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(0)
  , zold_(-1.) , d1old_(1.)
{
}

PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
  : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(0)
  , zold_(pkz.zold_) , d1old_(pkz.d1old_)
{
}

PkSpectrumZ::~PkSpectrumZ(void) 
{
}

void PkSpectrumZ::SetTypSpec(unsigned short typspec)
  // typsec = 0 : compute Pk(k)
  //        = 1 : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
{
  if(typspec>1) {
    cout<<"PkSpectrumZ::SetTypSpec: Error bad typspec value: "<<typspec<<endl;
    throw ParmError("PkSpectrumZ::SetTypSpec: Error bad typspec value");
  }
  typspec_ = typspec;
}

double PkSpectrumZ::operator() (double k)
{
  return (*this)(k,zref_);
}

double PkSpectrumZ::operator() (double k,double z)
{
  double d1;
  if(z == zold_) d1 = d1old_;
    else {d1 = d1old_ = d1_(z); zold_ = z;}

  double v = pk0_(k) * d1*d1;
  if(typspec_) v *= k*k*k/(2.*M_PI*M_PI);

  return scale_ * v;
}



///////////////////////////////////////////////////////////
//******************* VarianceSpectrum ******************//
///////////////////////////////////////////////////////////

VarianceSpectrum::VarianceSpectrum(GenericFunc& pk,unsigned short typfilter=0)
  : pk_(pk) , R_(0.)
{
  SetFilter(typfilter);
}

VarianceSpectrum::VarianceSpectrum(VarianceSpectrum& vpk)
  : pk_(vpk.pk_) , R_(vpk.R_)
{
  SetFilter(vpk.typfilter_);
}

VarianceSpectrum::~VarianceSpectrum(void)
{
}

//------------------------------------
void VarianceSpectrum::SetFilter(unsigned short typfilter)
// typfilter = 0 : spherical 3D top-hat
//           = 1 : spherical 3D gaussian
//           = 2 : no filter juste integrate spectrum)
{
  if(typfilter>2) {
    cout<<"VarianceSpectrum::SetFilter: Error bad value for type of filter"<<endl;
    throw ParmError("VarianceSpectrum::SetFilter: Error bad value for type of filter");
  }
  typfilter_ = typfilter;
}

void VarianceSpectrum::SetInteg(double dperc,double dlogkinc,double dlogkmax,unsigned short glorder)
// ATTENTION: on n'integre pas f(k)*dk mais k*f(k)*d(log10(k))
// see argument details in function IntegrateFuncLog (geneutils.cc)
{
  dperc_ = dperc;  if(dperc_<=0.) dperc_ = 0.1;
  dlogkinc_ = dlogkinc;
  dlogkmax_ = dlogkmax;
  glorder_ = glorder;
}


//------------------------------------
double VarianceSpectrum::Filter2(double x)
// ATTENTION: c'est le filtre au carre qui est renvoye
{
  // Just integrate the spectrum without filtering
  if(typfilter_ == 2) return 1.;

  double x2 = x*x;
  // Filtre gaussien G(x) = exp(-x^2/2)
  //        remarque G(x)^2 = exp(-x^2)
  // on prend le DL de G(x)^2 pour x->0 a l'ordre O(x^6)
  //             DL(x) = 1-x^2*(1-x^2/2)
  //             pour x<0.01  |DL(x)-G(X)^2|<2.0e-13
  if(typfilter_ == 1)
    if(x<0.01) return 1.-x2*(1.-x2/2.); else return exp(-x2);

  // Filtre top-hat T(x) = 3*(sin(x)-x*cos(x))/x^3
  // --- Gestion de la pseudo-divergence pour x->0
  // on prend le DL de T(x)^2 pour x->0 a l'ordre O(x^7)
  //             DL(x) = 1-x^2/5*(1-3*x^2/35*(1-4*x^2/81))
  //             pour x<0.1  |DL(x)-T(X)^2|<2.5e-13
  double f2=0.;
  if(x<0.1) {
    f2 = 1.-x2/5.*(1.-3.*x2/35.*(1.-4.*x2/81.));
  } else {
    f2 = 3.*(sin(x)-x*cos(x))/(x2*x);
    f2 *= f2;
  }
  return f2;
  
}

double VarianceSpectrum::Variance(double R,double kmin,double kmax)
// Compute variance of spectrum pk_ by integration
// Input:
//     R = taille du filter top-hat ou gaussien
//     kmin,kmax = bornes en k de l'integrale pour calculer la variance
// Return:
//     valeur de la variance (sigma^2)
// Remarque:
//   la meilleure approximation du filtre top-hat (R) est un filtre gaussien avec (Rg=R/sqrt(5))
//   la variance renvoyee est la variance de la masse
{
  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
    cout<<"VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax"<<endl;
    throw ParmError("VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax");
  }

  R_ = R;
  double lkmin = log10(kmin), lkmax = log10(kmax);

  double var = IntegrateFuncLog(*this,lkmin,lkmax,dperc_,dlogkinc_,dlogkmax_,glorder_);

  return var;
}

//------------------------------------
double VarianceSpectrum::FindMaximum(double R,double kmin,double kmax,double eps)
// Retourne le maximum de la fonction a integrer
// La recherche a lieu entre [kmin,kmax] par pas logarithmiques
// Input:
//     R : taille du filter top-hat ou gaussien
//     kmin,kmax : intervalle de recherche
//     eps : precision requise sur les valeurs
// Return:
//     position (en k) du maximum
{
  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
    cout<<"VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<<endl;
    throw ParmError("VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
  }

  R_ = R;

  int n = 10; // toujours >2
  double lkmin = log10(kmin), lkmax = log10(kmax), dlk = (lkmax-lkmin)/n;

  double lkfind=lkmin, pkfind=-1.;
  while(1) {
    for(int i=0; i<=n; i++) {
      double lk = lkmin  + i*dlk;
      double v = (*this)(pow(10.,lk));
      if(v<pkfind) continue;
      pkfind = v; lkfind = lk;
    }
    //cout<<"VarianceSpectrum::FindMaximum: lkfind="<<lkfind<<" pkfind="<<pkfind
    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
    // --- Convergence si l'encadrement de "kfind" est tel que "dk/kfind<eps"
    // On a dk = 10^(lkfind+dlk) - 10^(lkfind-dlk) = kfind * (10^(dlk) - 10^(-dlk))
    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
    dlk = (lkmax-lkmin)/n;
  }

  return pow(10.,lkfind);
}

int VarianceSpectrum::FindLimits(double R,double high,double &kmin,double &kmax,double eps)
// Retourne "[kmin,kmax]" tel que la fonction a integrer soit "f(k) <= high"
//   La recherche a lieu entre [kmin,kmax] par pas logarithmiques
// Input:
//     R : taille du filter top-hat ou gaussien
//     kmin,kmax : intervalle de recherche
//     eps : precision requise sur les valeurs kmin et kmax
// Output:
//     kmin,kmax telles que "f(k) <= high"
// Return:
//     rc  = 0 si OK
//     rc |= 1 "f(kmin) >= high"   (bit0 =1)
//     rc |= 2 "f(kmax) >= high"   (bit1 =1)
//     rc |= 4 "f(k) < high pour tout k"   (bit2 =1)
{
  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax  || eps<=0.) {
    cout<<"VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<<endl;
    throw ParmError("VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
  }

  R_ = R;
  int n = 10; // toujours >2

  int rc = 0;
  double lkmin,lkmax,dlk,lkfind;

  // --- Find kmin
  lkmin=log10(kmin); lkmax=log10(kmax); dlk=(lkmax-lkmin)/n;
  while(1) {
    lkfind = lkmin;
    for(int i=0;i<=n;i++) {
      if( (*this)(pow(10,lkfind)) >= high ) break;   //cmvbug
      lkfind = lkmin + i*dlk;
    }
    //cout<<"VarianceSpectrum::FindLimits[kmin]: lkfind="<<lkfind
    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
    if(fabs(lkfind-lkmax)<dlk/2.) {rc |= 4; return rc;} // protect against f(k)<high for all k
    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
    dlk = (lkmax-lkmin)/n;
  }
  if(lkfind-lkmin<dlk/2.) rc |= 1;  // f(kmin) >= high
  else kmin = pow(10.,lkmin);
  //cout<<"rc="<<rc<<" lkmin="<<lkmin<<"  pk="<<(*this)(pow(10.,lkmin))<<endl;

  // --- Find kmax
  lkmin=log10(kmin); lkmax=log10(kmax); dlk=(lkmax-lkmin)/n;
  while(1) {
    lkfind=lkmax;
    for(int i=0;i<=n;i++) {
      if( (*this)(pow(10,lkfind)) >= high ) break;   //cmvbug
      lkfind -= dlk;
      lkfind = lkmax - i*dlk;
    }
    //cout<<"VarianceSpectrum::FindLimits[kmax]: lkfind="<<lkfind
    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
    dlk = (lkmax-lkmin)/n;
  }
  if(lkmax-lkfind<dlk/2.) rc |= 2;  // f(kmax) >= high
  else kmax = pow(10.,lkmax);
  //cout<<"rc="<<rc<<" lkmax="<<lkmax<<"  pk="<<(*this)(pow(10.,lkmax))<<endl;

  return rc;
}