Changeset 3952 in Sophya for trunk/Cosmo

Timestamp:

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

Author:

cmv

Message:

ajout evolution de w, cmv 22/02/2011

Location:

trunk/Cosmo/SimLSS

Files:

• cmvtuniv.cc (modified) (7 diffs)
• cosmocalc.cc (modified) (9 diffs)
• cosmocalc.h (modified) (4 diffs)

trunk/Cosmo/SimLSS/cmvtuniv.cc

@@ -27 +27 @@
  <<"cmvtuniv [options] z1,z2"<<endl
  <<"  -D dz"<<endl
- <<"  -U h100,om0,or0,ol0,w0"<<endl
+ <<"  -U h100,om0,or0,ol0,w0,wa"<<endl
  <<"  -F flat(0,1,2)"<<endl
  <<"  -I perc,dzinc,dzmax,glorder"<<endl
@@ -44 +44 @@
  unsigned short flat = 0;
  // -- WMAP
- double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.,wa=0.;
  // -- ouvert matter only
- //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
+ //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1,wa=0..;
  // -- plat matter only
- //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.; flat = 1;
+ //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.,wa=0.; flat = 1;
  // -- plat lambda only
- //double h100=0.71, om0=0., or0=0., ol0=1.,w0=-1.; flat = 2;
+ //double h100=0.71, om0=0., or0=0., ol0=1.,w0=-1.,wa=0.; flat = 2;
 
  // -- parametre d'integration
@@ -70 +70 @@
     break;
   case 'U' :
-    sscanf(optarg,"%lf,%lf,%lf,%lf,%lf",&h100,&om0,&or0,&ol0,&w0);
+    sscanf(optarg,"%lf,%lf,%lf,%lf,%lf,%lf",&h100,&om0,&or0,&ol0,&w0,&wa);
     break;
   case 'F' :
@@ -94 +94 @@
  cout<<"z1="<<z1<<"  z2="<<z2<<"  dz="<<dz<<endl;
  cout<<"perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;
- cout<<"h100="<<h100<<" om0="<<om0<<" or0="<<or0<<" ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
+ cout<<"h100="<<h100<<" om0="<<om0<<" or0="<<or0<<" ol0="<<ol0<<" w0="<<w0<<" wa="<<wa<<" flat="<<flat<<endl;
  cout<<"todo ="<<todo<<endl;
 
  CosmoCalc univ1(flat,true,z2);
  univ1.SetInteg(perc,dzinc,dzmax,glorder);
- univ1.SetDynParam(h100,om0,or0,ol0,w0);
+ univ1.SetDynParam(h100,om0,or0,ol0,w0,wa);
  univ1.PrtInteg();
  univ1.Print();
@@ -108 +108 @@
    CosmoCalc univ2(flat,false,z2);
    univ2.SetInteg(perc,dzinc,dzmax,glorder);
-   univ2.SetDynParam(h100,om0,or0,ol0,w0);
+   univ2.SetDynParam(h100,om0,or0,ol0,w0,wa);
    univ2.PrtInteg();
    univ2.Print();
@@ -115 +115 @@
  else if(todo==3) tstspeed(univ1,z1,z2,dz);
  else if(todo==4) tstntuple(univ1,z1,z2,dz);
- else if(todo==5)tstinterp(univ1,z1,z2,dz);
+ else if(todo==5) tstinterp(univ1,z1,z2,dz);
  else cout<<"Unknown job to be done !!!!!  todo="<<todo<<endl;
 
@@ -215 +215 @@
    nt.Fill(xnt);
  }
+
  cout<<">>>> Ecriture"<<endl;
+ nt.Info()["flat"] = univ.Flat();
+ nt.Info()["h100"] = univ.h100();
+ nt.Info()["Olambda0"] = univ.Olambda0();
+ nt.Info()["W0"] = univ.W0();
+ nt.Info()["Wa"] = univ.Wa();
+ nt.Info()["Omatter0"] = univ.Omatter0();
+ nt.Info()["Obaryon0"] = univ.Obaryon0();
+ nt.Info()["Orelat0"] = univ.Orelat0();
+ nt.Info()["Ocurv0"] = univ.Ocurv0();
+ nt.Info()["Otot0"] = univ.Otot0();
+ nt.Info()["ZMax"] = univ.ZMax();
+
  string tag = "cmvtuniv_nt.ppf";
  POutPersist pos(tag);

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));
 }
 

trunk/Cosmo/SimLSS/cosmocalc.h

@@ -25 +25 @@
   void PrtInteg(void);
 
-  void SetDynParam(double h100,double om0,double or0,double ol0,double w0=-1.);
+  void SetDynParam(double h100,double om0,double or0,double ol0,double w0=-1.,double wa=0.);
   void SetObaryon0(double v);
   void DefaultParam(void);
@@ -34 +34 @@
   inline double Olambda0(void) {return _Olambda0;}
   inline double W0(void) {return _W0;}
+  inline double Wa(void) {return _Wa;}
   inline double Omatter0(void) {return _Omatter0;}
   inline double Obaryon0(void) {return _Obaryon0;}
@@ -43 +44 @@
   inline double H(double z) {return _H0 * E(z);}
   double Olambda(double z);
+  inline double W(double z) {return _W0 + _Wa * z/(1.+z);} // W(z) = W0 + Wa*(1-a)
   double Omatter(double z);
   double Obaryon(double z);
@@ -76 +78 @@
   unsigned short _flat;
   double _h100, _H0;  /* km/s/Mpc */
-  double _Olambda0,_W0;
+  double _Olambda0,_W0,_Wa;
   double _Omatter0,_Obaryon0;
   double _Orelat0;