Changeset 3314 in Sophya for trunk/Cosmo

Timestamp:

Aug 24, 2007, 5:00:00 PM (18 years ago)

Author:

cmv

Message:

correction bug (petit) dans fct transfert cmv 24/08/2007

Location:

trunk/Cosmo/SimLSS

Files:

• pkspectrum.cc (modified) (10 diffs)
• pkspectrum.h (modified) (2 diffs)

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -48 +48 @@
 ///////////////////////////////////////////////////////////
 
-// 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();
+// From Eisenstein & Hu ApJ 496:605-614 1998 April 1 (ou astro-ph/9709112)
+TransfertEisenstein::TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon,int lp)
+  : lp_(lp)
+  , Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb)
+  , nobaryon_(nobaryon) , nooscenv_(0), retpart_(0)
+{
+ zero_();
+ Init_();
 }
 
 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;
+  : lp_(tf.lp_)
+  ,Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , tcmb_(tf.tcmb_)
+  , nobaryon_(tf.nobaryon_) , nooscenv_(tf.nooscenv_), retpart_(tf.retpart_)
+{
+ zero_();
+ Init_();
+}
+
+void TransfertEisenstein::Init_(void)
+{
+
   O0_ = Oc_ + Ob_;
-  if(lp) cout<<"h100="<<h_<<" Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
+  if(nobaryon_) {O0_ = Oc_; Ob_ = 0.;}
+  if(lp_) cout<<"h100="<<h_<<" Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
 
   double H0 = 100. * h_, h2 = h_*h_;
 
+  if(tcmb_<0.) tcmb_ = T_CMB_Par;
   th2p7_ = tcmb_/2.7;
   double th2p7P4 = th2p7_*th2p7_*th2p7_*th2p7_;
+  if(lp_) cout<<"tcmb = "<<tcmb_<<" K = "<<th2p7_<<" *2.7K "<<endl;
 
   // Formule 2 p 606
   zeq_ = 2.50e4 * O0_ * h2 / th2p7P4;
-    if(lp) cout<<"zeq = "<<zeq_<<endl;
+  if(lp_) cout<<"zeq = "<<zeq_<<" (redshift of matter-radiation equality)"<<endl;
 
   // Formule 3 p 607
@@ -86 +89 @@
   //  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;
+  if(lp_) cout<<"keq = "<<keq_<<" Mpc^-1 (scale of equality)"<<endl;
+
+  // On s'arrete ici si pas de baryons
+  if(nobaryon_) return;
 
   // 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
+  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_<<" (Redshift of drag epoch)"<<endl;
+
+  // Formule 5 page 607    (R = 3*rho_baryon/4*rho_gamma)
   Req_ = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zeq_);
+  //WARNING: W.Hu code (tf_fit.c) en des-accord avec l'article: zd -> (1+zd)
   Rd_  = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zd_);
-    if(lp) cout<<"Req = "<<Req_<<" Rd = "<<Rd_<<endl;
+  //in tf_fit.c: Rd_  = 31.5*Ob_*h2 / th2p7P4 * (1.e3/(1.+zd_));
+  if(lp_) {
+    cout<<"Req = "<<Req_<<" Rd = "<<Rd_
+        <<" (Photon-baryon ratio at equality/drag epoch)"<<endl;
+    cout<<"Sound speed at equality "<<1./sqrt(3.*(1.+Req_))
+        <<", at drag "<<1./sqrt(3.*(1.+Rd_))<<" in unit of C"<<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;
+  if(lp_) cout<<"s = "<<s_<<" Mpc (sound horizon at drag epoch)"<<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;
+  if(lp_) cout<<"ksilk = "<<ksilk_<<" Mpc^-1 (silk damping scale)"<<endl;
 
   // Formules 10 page 608
@@ -116 +129 @@
   double b2 = pow(0.395*O0_*h2,-0.0266);
   betac_ = 1 / ( 1. + b1*(pow(Oc_/O0_,b2) - 1.) );
+  if(lp_) cout<<"alphac = "<<alphac_<<" betac = "<<betac_
+               <<" (CDM suppression/log shift)"<<endl;
 
   // Formule 23 page 610
   bnode_ = 8.41 * pow(O0_*h2,0.435);
-    if(lp) cout<<"bnode = "<<bnode_<<endl;
+  if(lp_) cout<<"bnode = "<<bnode_<<" (sound horizon shift)"<<endl;
 
   // Formule 14 page 608
+  //WARNING: W.Hu code (tf_fit.c) en des-accord avec l'article: (1+zeq) -> zeq
   double y = (1.+zeq_)/(1.+zd_);
+  //in tf_fit.c: double y = zeq_/(1.+zd_);
   double s1py = sqrt(1.+y);
   double Gy = y*( -6.*s1py + (2.+3.*y)*log((s1py+1.)/(s1py-1.)) );
@@ -130 +147 @@
   betab_ = 0.5 + Ob_/O0_
          + (3.-2.*Ob_/O0_) * sqrt(pow(17.2*O0_*h2,2.) + 1.);
+  if(lp_) cout<<"alphab = "<<alphab_<<" betab = "<<betab_
+               <<" (Baryon suppression/envelope shift)"<<endl;
 
   // Formule 31 page 612
@@ -135 +154 @@
           - 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;
+  if(lp_) cout<<"alphag = "<<alphag_<<" (gamma suppression in approximate TF)"<<endl;
+
+ // The approximate value of the sound horizon, formule 26 page 611
+ sfit_ = 44.5*log(9.83/(O0_*h2)) / sqrt(1.+10.*pow(Ob_*h2,3./4.));  // Mpc
+ if(lp_) cout<<"sfit="<<sfit_<<" Mpc (fit to sound horizon)"<<endl;
+
+ // La positoin du premier pic acoustique, formule 25 page 611
+ kpeak_ = 5*M_PI/(2.*sfit_) * (1.+0.217*O0_*h2);  // 1/Mpc
+ if(lp_) cout<<"kpeak="<<kpeak_<<" Mpc^-1 (fit to wavenumber of first peak)"<<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)
 {
@@ -161 +172 @@
 
 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
+// To obtain an approximate form of the non-oscillatory part of the transfert function
+// nooscenv = 0 : use the baryon oscillatory part of transfert function (full tf)
 // nooscenv = 1 : use approx. paragraph 3.3 p610 (middle of right column)
+//                Replace  j0(k*stilde)  ->  [1+(k*stilde)^4]^(-1/4)
 // 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;
+//                The value of an approximate transfer function that captures
+//                the non-oscillatory part of a partial baryon transfer function.
+//                In other words, the baryon oscillations are left out,
+//                but the suppression of power below the sound horizon is included.
+{
+ if(nooscenv!=1 && nooscenv!=2) nooscenv = 0;
+ nooscenv_ = nooscenv;
+}
+
+void TransfertEisenstein::SetReturnPart(unsigned short retpart)
+// To return only baryon or CDM part part of transfert function
+// retpart = 1 : return only CDM part of transfert function
+// retpart = 2 : return only Baryon part of transfert function
+// retpart = anything else: return only full transfert function
+// WARNING: only relevant for nobaryon_=false AND nooscenv!=2
+{
+ if(retpart!=1 && retpart!=2) retpart = 0;
+ retpart_ = retpart;
 }
 
@@ -192 +216 @@
   if(nobaryon_  || nooscenv_ == 2) {
     double gamma = O0_*h_;
-    // Formule 30 page 612 (pour fct lissee)
+    // Calcul de Gamma_eff, 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
+      gamma = O0_*h_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*sfit_,4.))); // Gamma_eff
+    // Formule 28 page 612 : qui est est equivalent a:
+    //    q = k / h_ * th2p7_*th2p7_ / gamma;
+    // qui est est equivalent a:
+    //    q = k / (13.41 * keq)                   pour Ob=0
+    //    q = k / (13.41 * keq) * (O0*h/Gamma)    pour le spectre lisse
+    // Les resultats sont legerement differents a cause des valeurs approx.
+    // des constantes numeriques: on prend comme W.Hu (tf_fit.c)
+    //double q = k / h_ * th2p7_*th2p7_ / gamma;  // Mpc^-1
+    double q = k/(13.41*keq_) * (O0_*h_/gamma);  // Mpc^-1
     // Formules 29 page 612
     double l0 = log(2.*M_E + 1.8*q);
@@ -203 +234 @@
   }
 
+  // --- Pour CDM + Baryons
   // --- CDM
   double f = 1. / (1. + pow(k*s_/5.4,4.));
   double Tc = f*T0tild(k,1.,betac_) + (1.-f)*T0tild(k,alphac_,betac_);
+  if(retpart_ == 1) return Tc; 
 
   // --- 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);
+  double stilde, ksbnode = k*s_/bnode_;
+  if(ksbnode<0.001) stilde =s_ * ksbnode;
+     else   stilde = 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(nooscenv_ == 1) {
+    j0kst = pow(1.+pow(k*stilde,4.) , -1./4.); //lissee sans oscillation baryon
+  } else {
+    double x = k*stilde;
     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;
+    //cout<<"DEBUG: k="<<k<<" stilde="<<stilde<<" x="<<x<<" j0kst="<<j0kst<<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 += alphab_/(1.+pow(betab_/(k*s_),3.)) * exp(-pow(k/ksilk_,1.4));
   Tb *= j0kst;
+  if(retpart_ == 2) return Tb;
 
   // --- Total
@@ -528 +563 @@
     lkfind = lkmin;
     for(int i=0;i<=n;i++) {
-      if( (*this)(pow(10,lkfind)) >= high ) break;   //cmvbug
+      if( (*this)(pow(10,lkfind)) >= high ) break;
       lkfind = lkmin + i*dlk;
     }
@@ -548 +583 @@
     lkfind=lkmax;
     for(int i=0;i<=n;i++) {
-      if( (*this)(pow(10,lkfind)) >= high ) break;   //cmvbug
+      if( (*this)(pow(10,lkfind)) >= high ) break;
       lkfind -= dlk;
       lkfind = lkmax - i*dlk;

trunk/Cosmo/SimLSS/pkspectrum.h

@@ -23 +23 @@
 class TransfertEisenstein : public GenericFunc {
 public:
-  TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon=false);
+  TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon=false,int lp=0);
   TransfertEisenstein(TransfertEisenstein& tf);
   virtual ~TransfertEisenstein(void);
@@ -29 +29 @@
   double KPeak(void);
   void SetNoOscEnv(unsigned short nooscenv=0);
+  void SetReturnPart(unsigned short retpart=0);
 protected:
-  double O0_,Oc_,Ob_,h_,tcmb_,th2p7_;
+  int lp_;
+  double O0_,Oc_,Ob_,h_,tcmb_;
+  double th2p7_;
   double zeq_,keq_,zd_,Req_,Rd_,s_,ksilk_,alphac_,betac_,bnode_,alphab_,betab_;
   double alphag_;
-  double kpeak_;
+  double sfit_,kpeak_;
 
   bool nobaryon_;
-  unsigned short nooscenv_;
+  unsigned short nooscenv_, retpart_;
+
   double T0tild(double k,double alphac,double betac);
-  void Init_With_Baryons(void);
-  void Init_Without_Baryon(void);
+  void Init_(void);
+  inline void zero_(void)
+    {th2p7_=zeq_=keq_=zd_=Req_=Rd_=s_=ksilk_=alphac_=betac_=bnode_=
+     alphab_=betab_=alphag_=sfit_=kpeak_=0.;}
 };