From abd736f596b09da7dc3a13777abc069d673fd52c Mon Sep 17 00:00:00 2001 From: "feroxity@yahoo.co.uk" Date: Wed, 27 Nov 2019 16:59:18 +0100 Subject: [PATCH] Regulate amplitude in ext_fld_atr is denom zero --- include/perturbations.h | 2 +- source/perturbations.c | 61 +++++++++++++++++++++++++++++------------ 2 files changed, 45 insertions(+), 18 deletions(-) diff --git a/include/perturbations.h b/include/perturbations.h index decdc1e4..8b50a9d2 100755 --- a/include/perturbations.h +++ b/include/perturbations.h @@ -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 } diff --git a/source/perturbations.c b/source/perturbations.c index 54fcc57e..a71f195f 100755 --- a/source/perturbations.c +++ b/source/perturbations.c @@ -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, &litude); 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) @@ -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, &litude), ppt->error_message,ppt->error_message); amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2 @@ -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, &litude), ppt->error_message,ppt->error_message); amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2. @@ -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, &litude), ppt->error_message,ppt->error_message); amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2 @@ -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, &litude), ppt->error_message,ppt->error_message); amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2 @@ -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, &litude), ppt->error_message,ppt->error_message); amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2 @@ -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, &litude), ppt->error_message,ppt->error_message); amplitude *=2; //calc_extfld_ampl assumes h normalised to 1/2 @@ -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!) @@ -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)0 && (fabs(den2)0 && (fabs(den3)0 && (fabs(den4)