Skip to content

Commit

Permalink
Merge pull request miguelzuma#3 from emiliobellini/new_parametrizations
Browse files Browse the repository at this point in the history 
New parametrizations added
  • Loading branch information
@miguelzuma
miguelzuma authored May 19, 2017
2 parents d4fa0e7 + 9cc77fa commit 7286064
Show file tree
Hide file tree
Showing 4 changed files with 115 additions and 44 deletions.
27 changes: 19 additions & 8 deletions hi_class.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,18 +79,29 @@ expansion_smg = 0.5 #this value will be overwritten using the closure equation.
ii) "propto_scale": all the alphas are proportional to the scale factor. Then, you have
to provide a vector containg the factor of proportionality of each alpha and the
initial value of the Planck Mass
iii) "planck_linear": this parametrization depends on one single parameter. It was used
by the Planck collaboration in 1502.01590 (this reference provides all the details
about this parametrization)
iv) "planck_exponential": this parametrization depends on two parameters. It was used
by the Planck collaboration in 1502.01590 (this reference provides all the details
about this parametrization)
iii) "eft_alphas_power_law": delta_M, alpha_{K, B, T} are proportional to the scale factor to
some power. Then, you have to provide a vector containg eight parameters, the four
proportionality constants and the four exponents of the scale factor. This generalize
the "propto_scale" parametrization with the remarkable exception that in that case
alpha_M was proportional to the scale factor, while here it is the Planck mass that
scales with the scale factor
iv) "eft_gammas_power_law": this parametrization uses the basis of the eft framework
(Omega, gamma_i) and assumes that they are power laws of the scale factor
(g1_0*pow(a,g1_exp)). You have to provide a vector containing eight parameters,
the four proportionality constants and the four exponents of the scale factor
v) "eft_gammas_exponential": this parametrization uses the basis of the eft framework
(Omega, gamma_i) and assumes that they are exponentials in the scale factor. In a
schematic way it reads: exp(gi_0*pow(a,gi_exp))-1. You have to provide a vector
containing eight parameters, the four proportionality constants and the four
exponents of the scale factor
The parametrization you want to used is stored in the variable "gravity_model", while the
value of the parameters is stored in the vector "parameters_smg"
i) "propto_omega" -> x_k, x_b, x_m, x_t, M*^2_ini (default)
ii) "propto_scale" -> x_k, x_b, x_m, x_t, M*^2_ini
iii) "planck_linear" -> Omega_0
iv) "planck_exponential" -> alpha_M0, beta
iii) "eft_alphas_power_law" -> delta_M_0, x_k, x_b, x_t, delta_M_0_exp, x_k_exp, x_b_exp, x_t_exp
iv) "eft_gammas_power_law" -> Om_0, g_1, g_2, g_3, Om_0_exp, g_1_exp, g_2_exp, g_3_exp
v) "eft_gammas_exponential" -> Om_0, g_1, g_2, g_3, Om_0_exp, g_1_exp, g_2_exp, g_3_exp


gravity_model = propto_omega
parameters_smg = 1., 0., 0., 0., 1.
Expand Down
2 changes: 1 addition & 1 deletion include/background.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
#include "parser.h"

enum spatial_curvature {flat,open,closed};
enum gravity_model {propto_omega, propto_scale, planck_linear, planck_exponential}; //write here the different models
enum gravity_model {propto_omega, propto_scale, eft_alphas_power_law, eft_gammas_power_law, eft_gammas_exponential}; //write here the different models

// initial conditions for the perturbations
enum pert_initial_conditions {single_clock, zero};
Expand Down
105 changes: 78 additions & 27 deletions source/background.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2195,16 +2195,17 @@ int background_initial_conditions(
pvecback_integration[pba->index_bi_M_pl_smg] = pba->parameters_2_smg[4];
break;

case planck_linear:
//NOTE: The Planck collaboration decided to consider models with M_pl^2=1+\Omega. Here, even if we have an additional integration parameter (i.e. M_pl_ini), we decided to fix it to M_pl^2=1+a*Omega_0 in order to be coherent with the choice of the Planck collaboration.
pvecback_integration[pba->index_bi_M_pl_smg] = 1+a*pba->parameters_2_smg[0];
case eft_alphas_power_law:
pvecback_integration[pba->index_bi_M_pl_smg] = 1. + pba->parameters_2_smg[0]*pow(a, pba->parameters_2_smg[4]);
break;

case planck_exponential:
//NOTE: The Planck collaboration decided to consider models with M_pl^2=1+\Omega. Here, even if we have an additional integration parameter (i.e. M_pl_ini), we decided to fix M_pl^2=exp(alpha_M0*pow(a, beta)/beta) in order to be coherent with the choice of the Planck collaboration.
pvecback_integration[pba->index_bi_M_pl_smg] = exp(pba->parameters_2_smg[0]*pow(a, pba->parameters_2_smg[1])/pba->parameters_2_smg[1]);
case eft_gammas_power_law:
pvecback_integration[pba->index_bi_M_pl_smg] = 1. + pba->parameters_2_smg[0]*pow(a,pba->parameters_2_smg[4]) + pba->parameters_2_smg[3]*pow(a,pba->parameters_2_smg[7]);
break;

case eft_gammas_exponential:
pvecback_integration[pba->index_bi_M_pl_smg] = exp(pba->parameters_2_smg[0]*pow(a,pba->parameters_2_smg[4])) + exp(pba->parameters_2_smg[3]*pow(a,pba->parameters_2_smg[7])) -1.;
break;

}

if (pba->M_pl_evolution_smg == _TRUE_)
Expand Down Expand Up @@ -2558,29 +2559,66 @@ int background_gravity_functions(
pvecback[pba->index_bg_mpl_running_smg] = c_m*a;
pvecback[pba->index_bg_M2_smg] = M_pl;
}
else if (pba->gravity_model_smg == planck_linear) {
//NOTE: With this parametrization every function it is expressed analytically. Then, it is possible to choose both to take the derivative of M_pl to obtain alpha_M or to integrate alpha_M to obtain M_pl. Even if the two results are undistinguishable, we choose the latter option, since in Class integrals are more stable numerically.
else if (pba->gravity_model_smg == eft_alphas_power_law) {

double Omega = a*pba->parameters_2_smg[0];
double M_0 = pba->parameters_2_smg[0];
double c_k = pba->parameters_2_smg[1];
double c_b = pba->parameters_2_smg[2];
double c_t = pba->parameters_2_smg[3];
double M_0_exp = pba->parameters_2_smg[4];
double c_k_exp = pba->parameters_2_smg[5];
double c_b_exp = pba->parameters_2_smg[6];
double c_t_exp = pba->parameters_2_smg[7];

pvecback[pba->index_bg_tensor_excess_smg] = 0.;
pvecback[pba->index_bg_mpl_running_smg] = Omega/(1.+Omega);
pvecback[pba->index_bg_braiding_smg] = -pvecback[pba->index_bg_mpl_running_smg];
pvecback[pba->index_bg_kineticity_smg] = 3.*((2.+3.*Omega)*(pvecback[pba->index_bg_rho_smg]+pvecback[pba->index_bg_p_smg])+Omega*(pvecback[pba->index_bg_rho_tot_wo_smg]+pvecback[pba->index_bg_p_tot_wo_smg]))/2./(1.+Omega)/rho_tot;
pvecback[pba->index_bg_kineticity_smg] = c_k*pow(a, c_k_exp);
pvecback[pba->index_bg_braiding_smg] = c_b*pow(a, c_b_exp);
pvecback[pba->index_bg_tensor_excess_smg] = c_t*pow(a, c_t_exp);
pvecback[pba->index_bg_mpl_running_smg] = M_0*M_0_exp*pow(a, M_0_exp)/(1. + M_0*pow(a, M_0_exp));
pvecback[pba->index_bg_M2_smg] = M_pl;
}
else if (pba->gravity_model_smg == planck_exponential) {
//NOTE: With this parametrization every function it is expressed analytically. Then, it is possible to choose both to take the derivative of M_pl to obtain alpha_M or to integrate alpha_M to obtain M_pl. Even if the two results are undistinguishable, we choose the latter option, since in Class integrals are more stable numerically.
else if ((pba->gravity_model_smg == eft_gammas_power_law) || (pba->gravity_model_smg == eft_gammas_exponential)) {

double Omega=0., g1=0., g2=0., g3=0., Omega_p=0., Omega_pp=0., g3_p=0.;
double Omega_0=0., g1_0=0., g2_0=0., g3_0=0.;
double Omega_exp=0., g1_exp=0., g2_exp=0., g3_exp=0.;

Omega_0 = pba->parameters_2_smg[0];
g1_0 = pba->parameters_2_smg[1];
g2_0 = pba->parameters_2_smg[2];
g3_0 = pba->parameters_2_smg[3];
Omega_exp = pba->parameters_2_smg[4];
g1_exp = pba->parameters_2_smg[5];
g2_exp = pba->parameters_2_smg[6];
g3_exp = pba->parameters_2_smg[7];

if (pba->gravity_model_smg == eft_gammas_power_law) {
Omega = Omega_0*pow(a,Omega_exp);
Omega_p = Omega_0*Omega_exp*pow(a,Omega_exp-1.); // Derivative w.r.t. the scale factor
Omega_pp = Omega_0*Omega_exp*(Omega_exp-1.)*pow(a,Omega_exp-2.); // Derivative w.r.t. the scale factor
g1 = g1_0*pow(a,g1_exp);
g2 = g2_0*pow(a,g2_exp);
g3 = g3_0*pow(a,g3_exp);
g3_p = g3_0*g3_exp*pow(a,g3_exp-1.); // Derivative w.r.t. the scale factor
}
else { //(pba->gravity_model_smg == eft_gammas_exponential)
Omega = exp(Omega_0*pow(a,Omega_exp))-1.;
Omega_p = Omega_0*Omega_exp*pow(a,Omega_exp-1.)*exp(Omega_0*pow(a,Omega_exp)); // Derivative w.r.t. the scale factor
Omega_pp = Omega_0*Omega_exp*pow(a,Omega_exp-2.)*exp(Omega_0*pow(a,Omega_exp))*(Omega_exp-1.+Omega_0*Omega_exp*pow(a,Omega_exp)); // Derivative w.r.t. the scale factor
g1 = exp(g1_0*pow(a,g1_exp))-1.;
g2 = exp(g2_0*pow(a,g2_exp))-1.;
g3 = exp(g3_0*pow(a,g3_exp))-1.;
g3_p = g3_0*g3_exp*pow(a,g3_exp-1.)*exp(g3_0*pow(a,g3_exp)); // Derivative w.r.t. the scale factor
}


double alpha_M0 = pba->parameters_2_smg[0];
double beta = pba->parameters_2_smg[1];
double Omega = exp(alpha_M0*pow(a, beta)/beta)-1;
double c_over_H2 = (-pow(a,2.)*Omega_pp + 3.*(Omega + a*Omega_p/2.)*(rho_tot + p_tot)/rho_tot + 3.*(pvecback[pba->index_bg_rho_smg]+pvecback[pba->index_bg_p_smg])/rho_tot)/2.;

pvecback[pba->index_bg_tensor_excess_smg] = 0;
pvecback[pba->index_bg_mpl_running_smg] = alpha_M0*pow(a, beta);
pvecback[pba->index_bg_braiding_smg] = -pvecback[pba->index_bg_mpl_running_smg];
pvecback[pba->index_bg_kineticity_smg] = (1-beta-pvecback[pba->index_bg_mpl_running_smg])*pvecback[pba->index_bg_mpl_running_smg] + 3*(2+pvecback[pba->index_bg_mpl_running_smg])*(pvecback[pba->index_bg_rho_smg]+pvecback[pba->index_bg_p_smg])/2/rho_tot + 3*(pvecback[pba->index_bg_mpl_running_smg]+2*Omega/(1+Omega))*(pvecback[pba->index_bg_rho_tot_wo_smg]+pvecback[pba->index_bg_p_tot_wo_smg])/2/rho_tot;
pvecback[pba->index_bg_kineticity_smg] = 2.*(2.*g1*pow(pba->H0,2.)/rho_tot + c_over_H2)/(1. + Omega + g3);
pvecback[pba->index_bg_braiding_smg] = -(g2*pba->H0/sqrt(rho_tot) + a*Omega_p)/(1. + Omega + g3);
pvecback[pba->index_bg_tensor_excess_smg] = -g3/(1. + Omega + g3);
pvecback[pba->index_bg_mpl_running_smg] = a*(Omega_p + g3_p)/(1. + Omega + g3);
pvecback[pba->index_bg_M2_smg] = M_pl;

}


Expand Down Expand Up @@ -2645,10 +2683,23 @@ int background_gravity_parameters(
pba->parameters_2_smg[4]);
break;

case planck_linear:
printf("Modified gravity: planck_linear with parameters: \n");
printf("-> Omega_0 = %g \n",
pba->parameters_2_smg[0]);
case eft_alphas_power_law:
printf("Modified gravity: eft_alphas_power_law with parameters: \n");
printf("-> M_*^2_0 = %g, c_K = %g, c_B = %g, c_T = %g, M_*^2_exp = %g, c_K_exp = %g, c_B_exp = %g, c_T_exp = %g\n",
pba->parameters_2_smg[0],pba->parameters_2_smg[1],pba->parameters_2_smg[2],pba->parameters_2_smg[3],
pba->parameters_2_smg[4],pba->parameters_2_smg[5],pba->parameters_2_smg[6],pba->parameters_2_smg[7]);
break;

case eft_gammas_power_law:
printf("Modified gravity: eft_gammas_power_law with parameters: \n");
printf("-> Omega_0 = %g, gamma_1 = %g, gamma_2 = %g, gamma_3 = %g, Omega_0_exp = %g, gamma_1_exp = %g, gamma_2_exp = %g, gamma_3_exp = %g \n",
pba->parameters_2_smg[0],pba->parameters_2_smg[1],pba->parameters_2_smg[2],pba->parameters_2_smg[3],pba->parameters_2_smg[4],pba->parameters_2_smg[5],pba->parameters_2_smg[6],pba->parameters_2_smg[7]);
break;

case eft_gammas_exponential:
printf("Modified gravity: eft_gammas_exponential with parameters: \n");
printf("-> Omega_0 = %g, gamma_1 = %g, gamma_2 = %g, gamma_3 = %g, Omega_0_exp = %g, gamma_1_exp = %g, gamma_2_exp = %g, gamma_3_exp = %g \n",
pba->parameters_2_smg[0],pba->parameters_2_smg[1],pba->parameters_2_smg[2],pba->parameters_2_smg[3],pba->parameters_2_smg[4],pba->parameters_2_smg[5],pba->parameters_2_smg[6],pba->parameters_2_smg[7]);
break;


Expand Down
25 changes: 17 additions & 8 deletions source/input.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1041,28 +1041,37 @@ int input_read_parameters(
class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg);
}

if (strcmp(string1,"planck_linear") == 0) {
pba->gravity_model_smg = planck_linear;
if (strcmp(string1,"eft_alphas_power_law") == 0) {
pba->gravity_model_smg = eft_alphas_power_law;
pba->field_evolution_smg = _FALSE_;
pba->M_pl_evolution_smg = _TRUE_;
flag2=_TRUE_;
pba->parameters_2_size_smg = 1;
pba->parameters_2_size_smg = 8;
class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg);
}

if (strcmp(string1,"planck_exponential") == 0) {
pba->gravity_model_smg = planck_exponential;
if (strcmp(string1,"eft_gammas_power_law") == 0) {
pba->gravity_model_smg = eft_gammas_power_law;
pba->field_evolution_smg = _FALSE_;
pba->M_pl_evolution_smg = _TRUE_;
flag2=_TRUE_;
pba->parameters_2_size_smg = 2;
pba->parameters_2_size_smg = 8;
class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg);
}


if (strcmp(string1,"eft_gammas_exponential") == 0) {
pba->gravity_model_smg = eft_gammas_exponential;
pba->field_evolution_smg = _FALSE_;
pba->M_pl_evolution_smg = _TRUE_;
flag2=_TRUE_;
pba->parameters_2_size_smg = 8;
class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg);
}


class_test(flag2==_FALSE_,
errmsg,
"could not identify gravity_theory value, check that it is one of 'propto_omega', 'propto_scale', 'planck_linear', 'planck_exponential' ...");
"could not identify gravity_theory value, check that it is one of 'propto_omega', 'propto_scale', 'eft_alphas_power_law', 'eft_gammas_power_law', 'eft_gammas_exponential' ...");

}// end of loop over models

Expand Down

0 comments on commit 7286064

Please sign in to comment.