Changeset 3952 in Sophya for trunk/Cosmo/SimLSS/cosmocalc.cc

Timestamp:

Feb 22, 2011, 4:32:25 PM (15 years ago)

Author:

cmv

Message:

ajout evolution de w, cmv 22/02/2011

File:

• trunk/Cosmo/SimLSS/cosmocalc.cc (modified) (9 diffs)

trunk/Cosmo/SimLSS/cosmocalc.cc

@@ -67 +67 @@
  _Olambda0 = univ._Olambda0;
  _W0 = univ._W0;
+ _Wa = univ._Wa;
  _Omatter0 = univ._Omatter0;
  _Obaryon0 = univ._Obaryon0;
@@ -149 +150 @@
 }
 
-void CosmoCalc::SetDynParam(double h100,double om0,double or0,double ol0,double w0)
+void CosmoCalc::SetDynParam(double h100,double om0,double or0,double ol0,double w0,double wa)
 {
   _computespl=true;
@@ -161 +162 @@
   _Olambda0 = ol0;
   _W0 = w0;
+  _Wa = wa;
   _Otot0 = _Olambda0 + _Omatter0 + _Orelat0;
   _Ocurv0 = 1. - _Otot0;
@@ -209 +211 @@
 
   double h100 = 0.71;
-  double ol0 = 0.73, w0 = -1.;
+  double ol0 = 0.73, w0 = -1., wa = 0.;
   double om0 = 0.135/(h100*h100);
   // Relat = photons (2.725 K) + neutrinos (1.9 K)
   double or0 = 2.47e-5 * (1. + 21./8.*pow(T_NU_Par/T_CMB_Par,4.)) / (h100*h100);
-  SetDynParam(h100,om0,or0,ol0,w0);
+  SetDynParam(h100,om0,or0,ol0,w0,wa);
   _Obaryon0 = 0.0224/(h100*h100);
 
@@ -224 +226 @@
   if(z<0.) z = 0.;
 
-  printf("CosmoCalc::Print(spl=%d,zmax=%g,flat=%u)  for z=%g\n",int(_usespline),ZMax(),Flat(),z);
+  printf("CosmoCalc::Print(spl=%d,zmax=%g,flat=%u)  for z=%g (a=%g)\n"
+         ,int(_usespline),ZMax(),Flat(),z,1./(1.+z));
   printf("h100=%g H0=%g Dhub=%g  H(z)=%g  Rhoc=%g g/cm^3 =%g Msol/Mpc^3\n"
         ,h100(),H0(),Dhubble(),H(z),Rhoc(z),Rhoc(z)*GCm3toMsolMpc3_Cst);
-  printf("Olambda=%g  W0=%g\n",Olambda(z),_W0);
+  printf("Olambda=%g  W=%g   (%g+(%g)*(1-a))\n",Olambda(z),W(z),_W0,_Wa);
   printf("Omatter=%g  Obaryon=%g\n",Omatter(z),Obaryon(z));
   printf("Orelat=%g\n",Orelat(z));
@@ -240 +243 @@
 //----------------------------------------------------------
 double CosmoCalc::Olambda(double z)
-{
-  double v = _Olambda0;
-  if(_W0 != -1.) v *= pow((1.+z),3.*(1+_W0));
+// Equation d'etat de l'energie noire: P=W(z)*Rho*C^2
+//   avec W(z) = W0 + Wa*(1-a) = W0 + Wa * z/(1+z) et a=1/(1+z)
+{
+  double zp1 = 1. + z;
+  double v = _Olambda0 * pow(zp1,3.*(1.+_W0+_Wa));
+  if(_Wa!=0.) v *= exp(-3.*_Wa*z/zp1);
   return v / E2(z);
 }
@@ -248 +254 @@
 double CosmoCalc::Omatter(double z)
 {
-  return _Omatter0 * (1.+z)*(1.+z)*(1.+z) /E2(z);
+  double zp1 = 1. + z;
+  return _Omatter0 * zp1*zp1*zp1 /E2(z);
 }
 
 double CosmoCalc::Obaryon(double z)
 {
-  return _Obaryon0 * (1.+z)*(1.+z)*(1.+z) /E2(z);
+  double zp1 = 1. + z;
+  return _Obaryon0 * zp1*zp1*zp1 /E2(z);
 }
 
@@ -269 +277 @@
 double CosmoCalc::Otot(double z)
 {
-  return Olambda(z)+Omatter(z)+Orelat(z);
+  return Olambda(z) + Omatter(z) + Orelat(z);
 }
 
@@ -369 +377 @@
 double CosmoCalc::E2(double z) const
 {
-  z += 1.;
+  double zp1 = 1. + z;
+
   double oldum = _Olambda0;
-
-  if(oldum>0. && _W0!=-1.) oldum *= pow(z,3.*(1.+_W0));
-
-  return oldum + z*z*(_Ocurv0 + z*(_Omatter0+z*_Orelat0));
+  if(oldum>0.) {
+    oldum *= pow(zp1,3.*(1.+_W0+_Wa));
+    if(_Wa!=0.) oldum *= exp(-3.*_Wa*z/zp1);
+  }
+
+  return oldum + zp1*zp1*(_Ocurv0 + zp1*(_Omatter0+zp1*_Orelat0));
 }