Skip to content

Commit

Permalink
Merge pull request #15 from wigster/hi_class
Browse files Browse the repository at this point in the history 
Regulate amplitude in ext_fld_attr if denominator zero
  • Loading branch information
@miguelzuma
miguelzuma authored Dec 3, 2019
2 parents 42f41ed + 486b9e8 commit e9110e5
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 18 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
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 e9110e5

Please sign in to comment.