Merge branch 'hi_class' of https://github.com/emiliobellini/hi_class_…

…teaser into emiliobellini-hi_class

* 'hi_class' of https://github.com/emiliobellini/hi_class_teaser:
  Corrected bug in the tensor equation
  Modified a bit stability test for D. Moved stability tests before the isnan test. Modified hi_class.ini

Conflicts:
	hi_class.ini

hi_class.ini

@@ -6,6 +6,8 @@
 > with detailed comments. The idea here is that the parameters used also in CLASS
 > are usually not reported, i.e. they assume their standard value. For a complete
 > list of these parameters look at the "explanatory.ini" file.
+> All the quantities related to the scalar field introduced in hi_class have been
+> marked with the "_smg" label (scalar modified gravity).
 > Only lines containing an equal sign not preceded by a sharp sign "#" are
 > considered by the code. Hence, do not write an equal sign within a comment,
 > the whole line would be interpreted as relevant input. Input files must have
@@ -16,13 +18,16 @@
 ----> background parameters:
 ----------------------------
 
-1a) Omega_smg is the fractional density of the scalar field today. It can be specified or
-    found automatically if Omega_Lambda and Omega_fld are given and Omega_smg set to -1.
-    This is to avoid someone modifying gravity by mistake. The code will then use Omega_smg
-    to satisfy the closure equation (sum_i Omega_i) equals (1 + Omega_k). If Omega_smg has
-    a positive (and smaller than 1) value and either Omega_Lambda or Omega_fld are left
-    unspecified, the code will use the closure equation to determine the unspecified
-    fractional density.
+1a) Omega_smg is the fractional density of the scalar field today. There are three possibilities
+    i)    Omega_smg unspecified or equivalently set to 0. In this case the code will
+          ignore the scalar field equations and will use the standard CLASS equations
+    ii)   Omega_smg has a value larger than 0 but smaller than 1. In this case you should
+          leave either Omega_Lambda or Omega_fld unspecified. Then, hi_class will run
+          with the scalar field equations, and Omega_Lambda or Omega_fld will be inferred
+          using the closure equation (sum_i Omega_i) equals (1 + Omega_k)
+    iii)  Omega_smg has a negqative value. In this case the equations for the scalar field
+          will be used, you have to specify both Omega_Lambda and Omega_fld, and Omega_smg will
+          be inferred by the code using the closure equation
 
 Omega_Lambda = 0 # fractional density of the cosmological constant today
 Omega_fld = 0 # fractional density of a perfect fluid today, see the "explanatory.ini" file for details on this fluid
@@ -45,8 +50,10 @@ Omega_smg = -1 # fractional density of the DE field today. If it is specified to
 
 expansion_model = lcdm
 expansion_smg = 0.5 #this value will be overwritten using the closure equation. See below
+# expansion_model = wowa
+# expansion_smg = 0.5, -1., 0.
 
-1d) Given a vector containing the values of the parameters needed to fully fix the background
+1c) Given a vector containing the values of the parameters needed to fully fix the background
     evolution you chose, you should specify which parameter you want to vary to satisfy the
     closure equation (tuning_index_smg, default: 0) and  your initial guess on how much you want 
     to vary it (tuning_dxdy_guess_smg, default: 1)
@@ -59,6 +66,7 @@ expansion_smg = 0.5 #this value will be overwritten using the closure equation.
 ----> gravity parameters:
 ------------------------------
 
+
 2a) The theory of gravity is described by the alpha functions at the linear level.
     There are four possible parametrizations of the alpha functions currently implemented
     i)    "propto_omega": all the alphas are proportional to the fractional density of the Dark
@@ -83,13 +91,13 @@ expansion_smg = 0.5 #this value will be overwritten using the closure equation.
 gravity_model = propto_omega
 parameters_smg = 1., 0., 0., 0., 1.
 
-
 2b) There are two possible initial conditions for the perturbations of the scalar field
     i)   "single_clock": single_clock IC given with respect to photons (default)
     ii)  "zero": set the initial value of the scalar field and its first derivative to 0
 
 pert_initial_conditions_smg = single_clock
 
+
 ---------------------------
 ----> stability parameters:
 ---------------------------
@@ -113,6 +121,7 @@ D_safe_smg = 0.
 ct2_safe_smg = 0.
 M2_safe_smg = 0. 
 
+
 ---------------------------
 ----> precision parameters:
 ---------------------------

source/background.c

@@ -1964,7 +1964,32 @@ int background_solve(
     }//end of has_smg   
 
   }
-  /* Yet another (third!) loop to make sure the background table makes sense
+
+  /* Horndeski stability tests
+  * only if not overriden
+  * and model is tuned!
+  * TODO: commented out till properly checked!
+  */
+  if ((pba->has_smg == _TRUE_) && 
+      (pba->parameters_tuned_smg == _TRUE_) &&
+      (pba->skip_stability_tests_smg == _FALSE_)){
+
+    class_test(pba->min_D_smg <= -fabs(pba->D_safe_smg),
+        pba->error_message,
+        "Ghost instability for scalar field perturbations with minimum D=%g \n",pba->min_D_smg);
+    class_test(pba->min_cs2_smg < -fabs(pba->cs2_safe_smg),
+        pba->error_message,
+        "Gradient instability for scalar field perturbations with minimum c_s^2=%g \n",pba->min_cs2_smg);
+    class_test(pba->min_M2_smg < -fabs(pba->M2_safe_smg),
+        pba->error_message,
+        "Ghost instability for metric tensor perturbations with minimum M*^2=%g \n",pba->min_M2_smg);
+    class_test(pba->min_ct2_smg < -fabs(pba->ct2_safe_smg),
+        pba->error_message,
+        "Gradient instability for metric tensor perturbations with minimum c_t^2=%g \n",pba->min_ct2_smg);
+
+   } 
+
+   /* Yet another (third!) loop to make sure the background table makes sense
    */
   for (i=0; i < pba->bt_size; i++) {
 
@@ -2034,30 +2059,6 @@ int background_solve(
   free(pvecback);
   free(pvecback_integration);
 
-  /* Horndeski stability tests
-   * only if not overriden
-   * and model is tuned!
-   * TODO: commented out till properly checked!
-   */
-   if ((pba->has_smg == _TRUE_) && 
-       (pba->parameters_tuned_smg == _TRUE_) &&
-       (pba->skip_stability_tests_smg == _FALSE_)){
-
-     class_test(pba->min_D_smg < -fabs(pba->D_safe_smg),
- 	       pba->error_message,
- 	       "Ghost instability for scalar field perturbations with minimum D=%g \n",pba->min_D_smg);
-     class_test(pba->min_cs2_smg < -fabs(pba->cs2_safe_smg),
- 	       pba->error_message,
- 	       "Gradient instability for scalar field perturbations with minimum c_s^2=%g \n",pba->min_cs2_smg);
-     class_test(pba->min_M2_smg < -fabs(pba->M2_safe_smg),
- 	       pba->error_message,
- 	       "Ghost instability for metric tensor perturbations with minimum M*^2=%g \n",pba->min_M2_smg);
-     class_test(pba->min_ct2_smg < -fabs(pba->ct2_safe_smg),
- 	       pba->error_message,
- 	       "Gradient instability for metric tensor perturbations with minimum c_t^2=%g \n",pba->min_ct2_smg);
-
-   } 
-
   return _SUCCESS_;
 
 }

source/perturbations.c

@@ -5347,10 +5347,11 @@ int perturb_einstein(
     /* modified version if gravity is non-standard. Note that no curvature is allowed in this case */
     else{
 
+      double M2 = ppw->pvecback[pba->index_bg_M2_smg];
       double run = ppw->pvecback[pba->index_bg_mpl_running_smg];
       double c_t2 = (1. + ppw->pvecback[pba->index_bg_tensor_excess_smg]);
 
-      ppw->pvecmetric[ppw->index_mt_gw_prime_prime] = -(2. + run)*a_prime_over_a*y[ppw->pv->index_pt_gwdot]-k2*c_t2*y[ppw->pv->index_pt_gw]+ppw->gw_source;
+      ppw->pvecmetric[ppw->index_mt_gw_prime_prime] = -(2. + run)*a_prime_over_a*y[ppw->pv->index_pt_gwdot]-k2*c_t2*y[ppw->pv->index_pt_gw]+ppw->gw_source/M2;
     }
 
   }