Merge branch 'devel'

* devel:
  updated version number to 2.3.1
  Corrected bugs related to dcdm
  Added support for omega_dcdm
  Moved dcdm parameters around inside background structure.

explanatory.ini

@@ -38,9 +38,10 @@ N_ur = 3.046
 Omega_cdm = 0.25
 #omega_cdm = 0.110616
 
-5b) density of dcdm+dr (decaying cold dark matter and its relativistic decay radiation): 'Omega_dcdmdr' (default: 'Omega_dcdmdr' set to 0)
+5b) density of dcdm+dr (decaying cold dark matter and its relativistic decay radiation): 'Omega_dcdmdr' or 'omega_dcdmdr' (default: 'Omega_dcdmdr' set to 0)
 
 Omega_dcdmdr = 0.0
+#omega_dcdmdr = 0.0
 
 5c) decay constant of dcdm in km/s/Mpc, same unit as H0 above.
 

include/background.h

@@ -47,14 +47,6 @@ struct background
 
   double Omega0_cdm; /**< \f$ \Omega_{0 cdm} \f$ : cold dark matter */
 
-  double Omega0_dcdmdr; /**< \f$ \Omega_{0 dcdm}+\Omega_{0 dr} \f$ : dcdm decaying to dr */
-
-  double Omega0_dcdm; /**< \f$ \Omega_{0 dcdm} \f$ : decaying cold dark matter */
-
-  double Gamma_dcdm; /**< \f$ \Gamma_{dcdm} \f$ : decay constant for decaying cold dark matter */
-
-  double Omega0_dr; /**< \f$ \Omega_{0 dr} \f$ : decay radiation */
-
   double Omega0_lambda; /**< \f$ \Omega_{0_\Lambda} \f$ : cosmological constant */
 
   double Omega0_fld; /**< \f$ \Omega_{0 de} \f$ : fluid with constant
@@ -69,6 +61,10 @@ struct background
 
   double Omega0_ur; /**< \f$ \Omega_{0 \nu r} \f$ : ultra-relativistic neutrinos */
 
+  double Omega0_dcdmdr; /**< \f$ \Omega_{0 dcdm}+\Omega_{0 dr} \f$ : dcdm decaying to dr */
+
+  double Gamma_dcdm; /**< \f$ \Gamma_{dcdm} \f$ : decay constant for decaying cold dark matter */
+
   double Omega0_k; /**< \f$ \Omega_{0_k} \f$ : curvature contribution */
 
   int N_ncdm;                            /**< Number of distinguishabe ncdm species */
@@ -119,6 +115,9 @@ struct background
   int sgnK; /**< K/|K|: -1, 0 or 1 */
   double * m_ncdm_in_eV; /**< list of ncdm masses in eV (infered from M_ncdm and other parameters above) */
   double Neff; /**< so-called "effective neutrino number", computed at earliest time in interpolation table */
+  double Omega0_dcdm; /**< \f$ \Omega_{0 dcdm} \f$ : decaying cold dark matter */
+  double Omega0_dr; /**< \f$ \Omega_{0 dr} \f$ : decay radiation */
+
 
   //@}
 

include/common.h

@@ -15,7 +15,7 @@
 #ifndef __COMMON__
 #define __COMMON__
 
-#define _VERSION_ "v2.3.0"
+#define _VERSION_ "v2.3.1"
 
 #define _TRUE_ 1 /**< integer associated to true statement */
 #define _FALSE_ 0 /**< integer associated to false statement */

include/input.h

@@ -134,10 +134,10 @@
  * temporary parameters for background fzero function
  */
 
-enum target_names {theta_s, Omega_dcdmdr};
+enum target_names {theta_s, Omega_dcdmdr, omega_dcdmdr};
 enum computation_stage {cs_background, cs_thermodynamics, cs_perturbations,
                         cs_primordial, cs_nonlinear, cs_transfer, cs_spectra};
-#define _NUM_TARGETS_ 2 //Keep this number as number of target_names
+#define _NUM_TARGETS_ 3 //Keep this number as number of target_names
 
 struct input_pprpba {
   struct precision * ppr;

source/input.c

@@ -223,9 +223,9 @@ int input_init(
   struct fzerofun_workspace fzw;
   /** These two arrays must contain the strings of names to be searched
       for and the coresponding new parameter */
-  char * const target_namestrings[] = {"100*theta_s","Omega_dcdmdr"};
-  char * const unknown_namestrings[] = {"h","Omega_ini_dcdm"};
-  enum computation_stage target_cs[] = {cs_thermodynamics, cs_background};
+  char * const target_namestrings[] = {"100*theta_s","Omega_dcdmdr","omega_dcdmdr"};
+  char * const unknown_namestrings[] = {"h","Omega_ini_dcdm","Omega_ini_dcdm"};
+  enum computation_stage target_cs[] = {cs_thermodynamics, cs_background, cs_background};
 
   int input_verbose = 0, int1;
 
@@ -748,18 +748,29 @@ int input_read_parameters(
   Omega_tot += pba->Omega0_cdm;
 
   /* Omega_0_dcdmdr (DCDM) */
-  class_read_double("Omega_dcdmdr",pba->Omega0_dcdmdr);
+  class_call(parser_read_double(pfc,"Omega_dcdmdr",&param1,&flag1,errmsg),
+             errmsg,
+             errmsg);
+  class_call(parser_read_double(pfc,"omega_dcdmdr",&param2,&flag2,errmsg),
+             errmsg,
+             errmsg);
+  class_test(((flag1 == _TRUE_) && (flag2 == _TRUE_)),
+             errmsg,
+             "In input file, you can only enter one of Omega_dcdmdr or omega_dcdmdr, choose one");
+  if (flag1 == _TRUE_)
+    pba->Omega0_dcdmdr = param1;
+  if (flag2 == _TRUE_)
+    pba->Omega0_dcdmdr = param2/pba->h/pba->h;
+  Omega_tot += pba->Omega0_dcdmdr;
+
   /** Omega_ini_dcdm will never be read from input, but is instead found automatically by shooting */
   class_read_double("Omega_ini_dcdm",pba->Omega_ini_dcdm);
 
-
   /* Read Gamma in same units as H0, i.e. km/(s Mpc)*/
   class_read_double("Gamma_dcdm",pba->Gamma_dcdm);
   /* Convert to Mpc */
   pba->Gamma_dcdm *= (1.e3 / _c_);
 
-  Omega_tot += pba->Omega0_dcdm;
-
   /* non-cold relics (ncdm) */
   class_read_int("N_ncdm",N_ncdm);
   if ((flag1 == _TRUE_) && (N_ncdm > 0)){
@@ -2308,6 +2319,7 @@ int input_default_params(
   pba->Omega0_ur = 3.046*7./8.*pow(4./11.,4./3.)*pba->Omega0_g;
   pba->Omega0_b = 0.02253/0.704/0.704;
   pba->Omega0_cdm = 0.1122/0.704/0.704;
+  pba->Omega0_dcdmdr = 0.0;
   pba->Omega0_dcdm = 0.0;
   pba->Gamma_dcdm = 0.0;
   pba->N_ncdm = 0;
@@ -2325,7 +2337,7 @@ int input_default_params(
   pba->Omega0_k = 0.;
   pba->K = 0.;
   pba->sgnK = 0;
-  pba->Omega0_lambda = 1.-pba->Omega0_k-pba->Omega0_g-pba->Omega0_ur-pba->Omega0_b-pba->Omega0_cdm-pba->Omega0_ncdm_tot-pba->Omega0_dcdm;
+  pba->Omega0_lambda = 1.-pba->Omega0_k-pba->Omega0_g-pba->Omega0_ur-pba->Omega0_b-pba->Omega0_cdm-pba->Omega0_ncdm_tot-pba->Omega0_dcdmdr;
   pba->Omega0_fld = 0.;
   pba->a_today = 1.;
   pba->w0_fld=-1.;
@@ -2825,37 +2837,6 @@ int get_machine_precision(double * smallest_allowed_variation) {
 
 }
 
-int input_fzerofun_for_background(double Omega_ini_dcdm,
-                                  void* container,
-                                  double *valout,
-                                  ErrorMsg error_message){
-
-  double rho_dcdm_today, rho_dr_today;
-  struct input_pprpba * pprpba;
-  struct background *pba;
-  pprpba = (struct input_pprpba *) container;
-  pba = pprpba->pba;
-
-  pba->Omega_ini_dcdm = Omega_ini_dcdm;
-  pba->keep_ncdm_stuff = _TRUE_;
-
-  class_call(background_init(pprpba->ppr,
-                             pba),
-             pba->error_message, error_message);
-  rho_dcdm_today = pba->background_table[(pba->bt_size-1)*pba->bg_size+pba->index_bg_rho_dcdm];
-  rho_dr_today = pba->background_table[(pba->bt_size-1)*pba->bg_size+pba->index_bg_rho_dr];
-  //  rho0 = pba->H0*pba->H0*pba->Omega0_ncdm[n_ncdm]; /*Remember that rho is defined such that H^2=sum(rho_i) */
-
-  *valout = (rho_dcdm_today+rho_dr_today)/(pba->H0*pba->H0)-pba->Omega0_dcdm;
-
-  if (pba->background_verbose > 3)
-    printf("rho_dcdm_today = %e, corresponding to %e\n",rho_dcdm_today,rho_dcdm_today/(pba->H0*pba->H0));
-
-  class_call(background_free(pba), pba->error_message, error_message);
-
-  return _SUCCESS_;
-}
-
 int input_fzerofun_1d(double input,
                       void* pfzw,
                       double *output,
@@ -3075,6 +3056,12 @@ int input_try_unknown_parameters(double * unknown_parameter,
       rho_dcdm_today = ba.background_table[(ba.bt_size-1)*ba.bg_size+ba.index_bg_rho_dcdm];
       rho_dr_today = ba.background_table[(ba.bt_size-1)*ba.bg_size+ba.index_bg_rho_dr];
       output[i] = (rho_dcdm_today+rho_dr_today)/(ba.H0*ba.H0)-pfzw->target_value[i];
+      break;
+    case omega_dcdmdr:
+      rho_dcdm_today = ba.background_table[(ba.bt_size-1)*ba.bg_size+ba.index_bg_rho_dcdm];
+      rho_dr_today = ba.background_table[(ba.bt_size-1)*ba.bg_size+ba.index_bg_rho_dr];
+      output[i] = (rho_dcdm_today+rho_dr_today)/(ba.H0*ba.H0)-pfzw->target_value[i]/ba.h/ba.h;
+      break;
     }
   }
 
@@ -3157,9 +3144,12 @@ int input_get_guess(double *xguess,
     case theta_s:
       xguess[index_guess] = 3.54*pow(pfzw->target_value[index_guess],2)-5.455*pfzw->target_value[index_guess]+2.548;
       dxdy[index_guess] = (7.08*pfzw->target_value[index_guess]-5.455);
+      /** Update pb to reflect guess */
+      ba.h = xguess[index_guess];
+      ba.H0 = ba.h *  1.e5 / _c_;
       break;
     case Omega_dcdmdr:
-      Omega_M = ba.Omega0_cdm+ba.Omega0_dcdm+ba.Omega0_b;
+      Omega_M = ba.Omega0_cdm+ba.Omega0_dcdmdr+ba.Omega0_b;
       /* This formula is exact in a Matter + Lambda Universe, but only
          for Omega_dcdm, not the combined.
          sqrt_one_minus_M = sqrt(1.0 - Omega_M);
@@ -3177,7 +3167,27 @@ int input_get_guess(double *xguess,
       dxdy[index_guess] = 1./a_decay;
       //printf("x = Omega_ini_guess = %g, dxdy = %g\n",*xguess,*dxdy);
       break;
+    case omega_dcdmdr:
+      Omega_M = ba.Omega0_cdm+ba.Omega0_dcdmdr+ba.Omega0_b;
+      /* This formula is exact in a Matter + Lambda Universe, but only
+         for Omega_dcdm, not the combined.
+         sqrt_one_minus_M = sqrt(1.0 - Omega_M);
+         xguess[index_guess] = pfzw->target_value[index_guess]*
+         exp(2./3.*ba.Gamma_dcdm/ba.H0*
+         atanh(sqrt_one_minus_M)/sqrt_one_minus_M);
+         dxdy[index_guess] = 1.0;//exp(2./3.*ba.Gamma_dcdm/ba.H0*atanh(sqrt_one_minus_M)/sqrt_one_minus_M);
+      */
+      gamma = ba.Gamma_dcdm/ba.H0;
+      if (gamma < 1)
+        a_decay = 1.0;
+      else
+        a_decay = pow(1+(gamma*gamma-1.)/Omega_M,-1./3.);
+      xguess[index_guess] = pfzw->target_value[index_guess]/ba.h/ba.h/a_decay;
+      dxdy[index_guess] = 1./a_decay/ba.h/ba.h;
+        //printf("x = Omega_ini_guess = %g, dxdy = %g\n",*xguess,*dxdy);
+      break;
     }
+    //printf("xguess = %g\n",xguess[index_guess]);
   }
 
   for (i=0; i<pfzw->fc.size; i++) {
@@ -3203,6 +3213,11 @@ int input_auxillary_target_conditions(enum target_names target_name, double targ
     if (target_value == 0.)
       return _FALSE_;
   }
+  if (target_name == omega_dcdmdr){
+    /* Check that Omega_dcdmdr is nonzero: */
+    if (target_value == 0.)
+      return _FALSE_;
+  }
   /* Default is no additional checks */
   return _TRUE_;
 }