Skip to content

Commit

Permalink
Merge remote-tracking branch 'miguelzuma/hi_class' into hi_class
Browse files Browse the repository at this point in the history
  • Loading branch information
Emilio Bellini committed Mar 3, 2020
2 parents 073b359 + 78fcc86 commit 502b835
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 21 deletions.
2 changes: 1 addition & 1 deletion include/perturbations.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -970,7 +970,7 @@ extern "C" {

int calc_extfld_ampl(int n, double kin, double bra, double dbra, double run, double ten, double DelM2,
double Omx, double wx, double l1, double l2, double l3, double l4,
double l5, double l6,double l7,double l8, double cs2num, double Dd,
double l5, double l6,double l7,double l8, double cs2num, double Dd, double ic_regulator_smg,
double * amplitude);
#ifdef __cplusplus
}
Expand Down
6 changes: 3 additions & 3 deletions source/background.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3477,7 +3477,7 @@ int background_gravity_functions(
if (pba->initial_conditions_set_smg == _FALSE_ || pba->hubble_evolution == _FALSE_){
class_test_except(E3*pow(E0,1./2.) > 1e-10 && pba->initial_conditions_set_smg == _FALSE_,
pba->error_message,
background_free(pba);free(pvecback);free(pvecback_B),
free(pvecback);free(pvecback_B),
" E3=%e is large in Friedmann constraint when setting ICs ", E3);
/* Use Newton's method */
x = sqrt(E0);
Expand Down Expand Up @@ -3508,7 +3508,7 @@ int background_gravity_functions(

class_test_except(isnan(pvecback[pba->index_bg_H]),
pba->error_message,
background_free(pba);free(pvecback);free(pvecback_B),
free(pvecback);free(pvecback_B),
" H=%e is not a number at a = %e. phi = %e, phi_prime = %e, E0=%e, E1=%e, E2=%e, E3=%e, M_*^2 = %e ",
pvecback[pba->index_bg_H],a,phi,phi_prime,E0,E1,E2,E3,pvecback[pba->index_bg_M2_smg]);

Expand Down Expand Up @@ -3542,7 +3542,7 @@ int background_gravity_functions(

class_test_except((A*M - B*F) == 0 ,
pba->error_message,
background_free(pba);free(pvecback);free(pvecback_B),
free(pvecback);free(pvecback_B),
"scalar field mixing with metric has degenerate denominator at a = %e, phi = %e, phi_prime = %e \n with A = %e, M =%e, B=%e, F=%e, \n H=%e, E0=%e, E1=%e, E2=%e, E3=%e \n M_*^2 = %e Kineticity = %e, Braiding = %e",
a,phi,phi_prime, A, M, B, F,
pvecback[pba->index_bg_H],E0,E1,E2,E3,
Expand Down
61 changes: 44 additions & 17 deletions source/perturbations.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5027,11 +5027,11 @@ int perturb_initial_conditions(struct precision * ppr,


calc_extfld_ampl(nexpo, kin, bra, dbra, run, ten, DelM2, Omx, wx,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd, ppr->pert_ic_regulator_smg,
&amplitude);

ppw->pv->y[ppw->pv->index_pt_vx_smg] = amplitude*ktau_two*tau*(ppr->curvature_ini);
ppw->pv->y[ppw->pv->index_pt_vx_prime_smg] = nexpo*a*ppw->pvecback[pba->index_bg_H]*ppw->pv->y[ppw->pv->index_pt_vx_smg];
ppw->pv->y[ppw->pv->index_pt_vx_prime_smg] = (nexpo+1)*a*ppw->pvecback[pba->index_bg_H]*ppw->pv->y[ppw->pv->index_pt_vx_smg];


if(ppt->perturbations_verbose > 5)
Expand Down Expand Up @@ -5219,7 +5219,7 @@ int perturb_initial_conditions(struct precision * ppr,
coeff_isocurv_smg = ppr->entropy_ini*fraccdm*om;

class_call(calc_extfld_ampl(nexpo, kin, bra, dbra, run, ten, DelM2, Omx, wx,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd, ppr->pert_ic_regulator_smg,
&amplitude),
ppt->error_message,ppt->error_message);
amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2
Expand Down Expand Up @@ -5274,7 +5274,7 @@ int perturb_initial_conditions(struct precision * ppr,
coeff_isocurv_smg = ppr->entropy_ini*fracb*om;

class_call(calc_extfld_ampl(nexpo, kin, bra, dbra, run, ten, DelM2, Omx, wx,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd, ppr->pert_ic_regulator_smg,
&amplitude),
ppt->error_message,ppt->error_message);
amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2.
Expand Down Expand Up @@ -5331,7 +5331,7 @@ int perturb_initial_conditions(struct precision * ppr,
coeff_isocurv_smg = ppr->entropy_ini * fracb*fracnu/fracg/10.*k*k * om/4;

class_call(calc_extfld_ampl(nexpo, kin, bra, dbra, run, ten, DelM2, Omx, wx,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd, ppr->pert_ic_regulator_smg,
&amplitude),
ppt->error_message,ppt->error_message);
amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2
Expand All @@ -5345,7 +5345,7 @@ int perturb_initial_conditions(struct precision * ppr,
(-32.*k*k/(15.+4.*fracnu)- 9.*fracb*(fracb+fracg)*om*om/fracg/fracg);

class_call(calc_extfld_ampl(nexpo, kin, bra, dbra, run, ten, DelM2, Omx, wx,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd, ppr->pert_ic_regulator_smg,
&amplitude),
ppt->error_message,ppt->error_message);
amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2
Expand Down Expand Up @@ -5402,7 +5402,7 @@ int perturb_initial_conditions(struct precision * ppr,
coeff_isocurv_smg = ppr->entropy_ini * 9./32. *k*om*fracnu*fracb/fracg;

class_call(calc_extfld_ampl(nexpo, kin, bra, dbra, run, ten, DelM2, Omx, wx,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd, ppr->pert_ic_regulator_smg,
&amplitude),
ppt->error_message,ppt->error_message);
amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2
Expand All @@ -5416,7 +5416,7 @@ int perturb_initial_conditions(struct precision * ppr,
( -3.*om*om*k/160.*fracb*(3*fracb+5*fracg)/fracg/fracg -4*k*k*k/15./(5.+4.*fracnu) );

class_call(calc_extfld_ampl(nexpo, kin, bra, dbra, run, ten, DelM2, Omx, wx,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd,
l1, l2, l3, l4, l5, l6,l7,l8, cs2num, Dd, ppr->pert_ic_regulator_smg,
&amplitude),
ppt->error_message,ppt->error_message);
amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2
Expand Down Expand Up @@ -10701,7 +10701,7 @@ int perturb_test_ini_extfld_ic_smg(struct precision * ppr,
int calc_extfld_ampl(int nexpo, double kin, double bra, double dbra, double run, double ten, double DelM2,
double Omx, double wx, double l1, double l2, double l3, double l4,
double l5, double l6,double l7,double l8, double cs2num, double Dd,
double * amplitude){
double ic_regulator_smg, double * amplitude){

/* Solutions assuming the alphas are small, i.e. Vx does not gravitate but moves
* on an attractor provided by collapsing radiation. (w!=1/3 terms included properly here!)
Expand Down Expand Up @@ -10731,32 +10731,59 @@ int calc_extfld_ampl(int nexpo, double kin, double bra, double dbra, double run


double B1_smg, B2_smg, B3_smg, B3num_smg, B3denom_smg;
double den1, den2, den3, den4, reg_rescaled;

B1_smg = (bra/Dd)*(bra/(2.*(-2 + bra)*(kin + l1)))*((-6 + kin)*l1 + 3*l4);
B1_smg += (3*pow(bra,3))*(l1/Dd)/(2.*(-2 + bra)*(kin + l1));
reg_rescaled = ic_regulator_smg*(fabs(bra)+fabs(kin)+fabs(l1)); //rescale the regulator to be proportional to the alphas


den1 = (2.*(-2 + bra)*(kin + l1));
if(reg_rescaled>0 && (fabs(den1)<reg_rescaled)){
den1 = copysign(reg_rescaled,den1);
}

B1_smg = (bra/Dd)*(bra/den1)*((-6 + kin)*l1 + 3*l4);
B1_smg += (3*pow(bra,3))*(l1/Dd)/den1;
B1_smg += 2*(cs2num/Dd)*(3*bra*kin + pow(kin,2) - 3*l4)/(2.*(-2. + bra)*(kin + l1));
B1_smg += 2*(3*l2*l4/Dd + (kin/Dd)*(l1*l2 - 8*l7) - 8*l1/Dd*l7)/(2.*(-2 + bra)*(kin + l1));
B1_smg += 2*(3*l2*l4/Dd + (kin/Dd)*(l1*l2 - 8*l7) - 8*l1/Dd*l7)/den1;
B1_smg -= 2*(bra/Dd)*((kin*l1/(kin + l1) - 3*l1*l2/(kin + l1) + 3*l4/(kin + l1))/(2.*(-2 + bra)));

den2 = (4.*(-2 + bra)*(1 + DelM2)*(kin + l1));
if(reg_rescaled>0 && (fabs(den2)<reg_rescaled)){
den2 = copysign(reg_rescaled,den2);
}

B2_smg = 8*(1 + DelM2)*(3*l2*l6/Dd + 4*kin*l8/Dd);
B2_smg += 4*(l1/Dd)*(8*(1 + DelM2)*l8 + l2*(12 - 12*Omx + (1 + DelM2)*(-12 + kin + Omx*(3 - 9*wx))));
B2_smg += 2*(bra/Dd)*bra*(6*(1 + DelM2)*l6 + l1*(12 - 12*Omx + (1 + DelM2)*(-30 + kin + 6*Omx*(1 - 3*wx))));
B2_smg += 3*pow(bra,3)*(1 + DelM2)*(l1/Dd)*(6 + Omx*(-1 + 3*wx));
B2_smg += 2*(cs2num/Dd)*(2*(1 + DelM2)*pow(kin,2) - 12*(1 + DelM2)*l6 + 3*kin*(8 - 8*Omx + (1 + DelM2)*(-8 + Omx*(2 - 6*wx) + bra*(6 + Omx*(-1 + 3*wx)))));
B2_smg -= 2*(bra/Dd)*(12*(1 + DelM2)*l6 + l1*(24 - 24*Omx + (1 + DelM2)*(2*kin - 3*(8 + 2*Omx*(-1 + 3*wx) + l2*(6 + Omx*(-1 + 3*wx))))));
B2_smg /= (4.*(-2 + bra)*(1 + DelM2)*(kin + l1));
B2_smg /= den2;


den3 = ((2. * Omx)*(kin + l1));
reg_rescaled *=Omx;
if(reg_rescaled>0 && (fabs(den3)<reg_rescaled)){
den3 = copysign(reg_rescaled,den3);
}
B3num_smg = ((-(((-2. + bra) * bra + 2 * l2) *
((-2. + bra) * l1 - 4 * l3 + 2 * Dd * (-1. + nexpo))) +
cs2num * (-2 * (-2. + bra) * kin - 8 * l3 + 4 * Dd * (-1. + nexpo))) *
nexpo) /
((2. * Omx)*(kin + l1));
nexpo) / den3;

B3denom_smg = 4*(Dd/Omx)*(-2 + bra);

B3_smg = B3num_smg/B3denom_smg;

*amplitude = -(B3_smg/(B1_smg + B2_smg + nexpo + B1_smg*nexpo + pow(nexpo,2)));
reg_rescaled = ic_regulator_smg*(fabs(B1_smg)+fabs(B2_smg));
den4 = B1_smg + B2_smg + nexpo + B1_smg*nexpo + pow(nexpo,2);

if(reg_rescaled>0 && (fabs(den4)<reg_rescaled)){
den4 = copysign(reg_rescaled,den4);
}


*amplitude = -B3_smg/den4;

return _SUCCESS_;
}
}

0 comments on commit 502b835

Please sign in to comment.