Skip to content

New horndeski #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
kpardo wants to merge 5 commits into
base: master
Choose a base branch
from
Open

New horndeski #3

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1cc5abc
Initial commit of new_horndeski branch
Feb 14, 2018
7eddf13
Created separate module for modgrav routines
Feb 15, 2018
cbec3bf
Made changes according to Phil's comments; Also added an error if mod…
Feb 22, 2018
b90846a
Changed horndeski condition to ignore alpha_K. alpha_k = 1 is needed …
Feb 23, 2018
21cf815
Fixed small typo related to alpha_T
Mar 5, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 9 additions & 12 deletions theory/cosmo2D_fourier.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,24 +17,21 @@ double C_shear_tomo(double l, int ni, int nj); //shear tomography power spectra
double C_shear_tomo_nointerp(double l, int ni, int nj);
/**********************************************************/

double MG_Sigma(double a)
{
// double aa=a*a;
// double omegam=cosmology.Omega_m/(aa*a);
double omegav=omv_vareos(a);
double hub=hoverh0(a);
hub = hub*hub;

return cosmology.MGSigma*omegav/hub/cosmology.Omega_v;
}

double dchi_da(double a){
return 1./(a*a*hoverh0(a));
}
double W_kappa(double a, double fK, double nz){
double wkappa = 1.5*cosmology.Omega_m*fK/a*g_tomo(a,(int)nz);
double aa=a*a;
Copy link
Copy Markdown

@philbull philbull Feb 16, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably don't need to define this separately, as aa is only used once.

double omegam=cosmology.Omega_m/(aa*a);
double omegav=omv_vareos(a);
double hub = hoverh0(a);
hub = hub*hub;
if(cosmology.MGSigma != 0.){
Copy link
Copy Markdown

@philbull philbull Feb 16, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be enclosed in a preprocessor macro.

wkappa *= (1.+MG_Sigma(a));
wkappa *= (1.+MG_sigma(a, omegav, hub));
}
else if (cosmology.MGalpha_K != 0 || cosmology.MGalpha_B != 0 || cosmology.MGalpha_T != 0 || cosmology.MGalpha_M != 0){
wkappa*= (1.+horndeski_sigma(a, omegav, hub, omegam));
}
return wkappa;
}
Expand Down
5 changes: 4 additions & 1 deletion theory/cosmo3D.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,10 @@ int func_for_growfac(double a,const double y[],double f[],void *params)
hub = hub*hub;
f[0]=y[1];
if(cosmology.MGmu != 0){
Copy link
Copy Markdown

@philbull philbull Feb 16, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can also be enclosed in preprocessor macros (IFDEF etc).

Copy link
Copy Markdown
Author

@kpardo kpardo Feb 17, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why would you use ifdef instead of define here?

one_plus_mg_mu += cosmology.MGmu*omegav/hub/cosmology.Omega_v;
one_plus_mg_mu += MG_mu(omegav, hub);
}
else if (cosmology.MGalpha_K != 0 || cosmology.MGalpha_B != 0 || cosmology.MGalpha_T != 0 || cosmology.MGalpha_M != 0){
one_plus_mg_mu += horndeski_mu(a, omegav, hub, omegam);
}
f[1]=y[0]*3.*cosmology.Omega_m/(2.*hub*aa*aa*a)*one_plus_mg_mu-y[1]/a*(2.-(omegam+(3.*(cosmology.w0+cosmology.wa*(1.-a))+1)*omegav)/(2.*hub));
return GSL_SUCCESS;
Expand Down
81 changes: 81 additions & 0 deletions theory/modgrav.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
double MG_mu(double omegav, double hub);
double MG_sigma(double a, double omegav, double hub);
double horndeski_mu(double a, double omegav, double hub, double omegam);
double horndeski_sigma(double a, double omegav, double hub, double omegam);
double horndeski_cs_sqrd(double a, double fz, double fzprime, double H, double Hprime, double omegam, double hub);
double horndeski_alpha(double a, double fz);
/**********************************************************/

double MG_mu(double omegav, double hub){
return cosmology.MGmu*omegav/hub/cosmology.Omega_v;
}

double MG_sigma(double a, double omegav, double hub){
return cosmology.MGSigma*omegav/hub/cosmology.Omega_v;
}
double horndeski_mu(double a, double omegav, double hub, double omegam){
//this is all in the quasistatic limit!!
//given by Eqn 43 in Tessa Baker notes
//!!!!currently assuming Mpl = Mstar!!!!
double fz = omegav/hub*1./cosmology.Omega_v;
Copy link
Copy Markdown

@philbull philbull Feb 16, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add a comment saying how this is defined, please? We might want to use a different definition at some point, so good to be clear what this is.

double wde = cosmology.w0 + cosmology.wa*(1. - a);
double H = (100. * cosmology.h0 / constants.lightspeed) * sqrt(hub);
double Hprime = -0.5*H*H*(omegam/hub+omegav/hub*(1.+3.*wde))-H*H;
double fzprime = -3.*wde * a*H * omegam/hub * fz;
double alpha = horndeski_alpha(a, fz);
double cs_sqrd = horndeski_cs_sqrd(a, fz, fzprime, H, Hprime,omegam, hub);
/////////////
double alpha_K = cosmology.MGalpha_K*fz;
double alpha_B = cosmology.MGalpha_B*fz;
double alpha_T = cosmology.MGalpha_T*fz;
double alpha_M = cosmology.MGalpha_M*fz;
////////////
double mu_term1 = alpha*cs_sqrd*(1+alpha_T)*(1+alpha_T);
double mu_term2 = (-alpha_B*(1+alpha_T)+ 2*(alpha_T - alpha_M))*(-alpha_B*(1+alpha_T)+ 2*(alpha_T - alpha_M));
double mu = (mu_term1+mu_term2)/(alpha*cs_sqrd) - 1.;
return mu;
}

double horndeski_sigma(double a, double omegav, double hub, double omegam){
//this is all in the quasistatic limit!!
//also assumes Mpl = Mstar!!!!
//alpha function redshift scaling
double fz = omegav/hub*1./cosmology.Omega_v;
double wde = cosmology.w0 + cosmology.wa*(1. - a);
double H = (100. * cosmology.h0 / constants.lightspeed) * sqrt(hub);
double Hprime = -0.5*H*H*(omegam/hub+omegav/hub*(1.+3.*wde))-H*H;
double fzprime = -3.*wde * a*H * omegam/hub * fz;
//define actual alpha functions: fz*alpha_X,0
double alpha_K = cosmology.MGalpha_K*fz;
double alpha_B = cosmology.MGalpha_B*fz;
double alpha_T = cosmology.MGalpha_T*fz;
double alpha_M = cosmology.MGalpha_M*fz;
//helpful functions (see Tessa Baker's notes)
double alpha = horndeski_alpha(a, fz);
double cs_sqrd = horndeski_cs_sqrd(a, fz, fzprime, H, Hprime, omegam, hub);
//given by gamma Eqn. 44 in Tessa Baker notes
double num = alpha*cs_sqrd-alpha_B*(-alpha_B/(2.*(1.+alpha_T))+alpha_T-alpha_M);
double denom = alpha*cs_sqrd*(1.+alpha_T)+(-alpha_B*(1.+alpha_T)+2.*(alpha_T-alpha_M))*(-alpha_B*(1.+alpha_T)+2.*(alpha_T-alpha_M));
double mu = horndeski_mu(a, omegav, hub, omegam);
return 0.5*(1+mu)*(num/denom+1.);
}
double horndeski_cs_sqrd(double a, double fz, double fzprime, double H, double Hprime, double omegam, double hub){
double alpha_K = cosmology.MGalpha_K*fz;
double alpha_B = cosmology.MGalpha_B*fz;
double alpha_T = cosmology.MGalpha_T*fz;
double alpha_M = cosmology.MGalpha_M*fz;

double alpha = horndeski_alpha(a, fz);
double cs_term1 = (2.-alpha_B)*(Hprime - (alpha_M -alpha_T)*H*H-H*H*alpha_B/(2.*(1.+alpha_T)));
Copy link
Copy Markdown

@philbull philbull Feb 16, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There seems to be a common factor of H*H for the last two terms. Is there any reason why they're broken apart?
(Also, why the weird indent?)

//!!!!!!!!Assuming Mpl = Mstar!!!!!!!!!!!!!
double cs_term2 = -H*cosmology.MGalpha_B*fzprime + 1.5*H*H*omegam/hub;
double cs_sqrd = -(cs_term1 + cs_term2)/(H*H*alpha);
return cs_sqrd;
}
double horndeski_alpha(double a, double fz){
Copy link
Copy Markdown

@philbull philbull Feb 16, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function could be declared inline, which may help speed things up very slightly.

//Helpful function -- see Tessa Baker's notes
double alpha_K = cosmology.MGalpha_K*fz;
double alpha_B = cosmology.MGalpha_B*fz;
return alpha_K + 1.5*alpha_B*alpha_B;
}

1 change: 0 additions & 1 deletion theory/parameters.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -737,7 +737,6 @@ void set_cosmological_parameters_to_MILLENIUM()
printf("Cosmology set to Sato\n");
}


void set_survey_parameters_to_DES_Y1()
{
survey.area = 1000.;
Expand Down
8 changes: 6 additions & 2 deletions theory/recompute.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,10 @@ void update_cosmopara (cosmopara *C){
C->MGSigma = cosmology.MGSigma;
C->MGmu = cosmology.MGmu;
C->M_nu = cosmology.M_nu;
C->MGalpha_K = cosmology.MGalpha_K;
C->MGalpha_B = cosmology.MGalpha_B;
C->MGalpha_T = cosmology.MGalpha_T;
C->MGalpha_M = cosmology.MGalpha_M;
}

void update_galpara (galpara *G){
Expand Down Expand Up @@ -87,7 +91,7 @@ void update_nuisance (nuisancepara *N){
N->cluster_centering_sigma = nuisance.cluster_centering_sigma;
}
int recompute_expansion(cosmopara C){ //rules for recomputing growth factor & comoving distance
if (C.Omega_m != cosmology.Omega_m || C.Omega_v != cosmology.Omega_v || C.w0 != cosmology.w0 || C.wa != cosmology.wa || C.MGmu != cosmology.MGmu || C.M_nu != cosmology.M_nu){return 1;}
if (C.Omega_m != cosmology.Omega_m || C.Omega_v != cosmology.Omega_v || C.w0 != cosmology.w0 || C.wa != cosmology.wa || C.MGmu != cosmology.MGmu || C.M_nu != cosmology.M_nu || C.MGalpha_K != cosmology.MGalpha_K ||C.MGalpha_B != cosmology.MGalpha_B ||C.MGalpha_T != cosmology.MGalpha_T ||C.MGalpha_M != cosmology.MGalpha_M){return 1;}
else{return 0;}
}

Expand All @@ -103,7 +107,7 @@ int recompute_Delta(cosmopara C){ //rules for recomputing early time power spect
}

int recompute_cosmo3D(cosmopara C){
if (C.Omega_m != cosmology.Omega_m || C.Omega_v != cosmology.Omega_v || C.Omega_nu != cosmology.Omega_nu || C.M_nu != cosmology.M_nu || C.h0 != cosmology.h0 || C.omb != cosmology.omb || C.n_spec != cosmology.n_spec|| C.alpha_s != cosmology.alpha_s || C.w0 != cosmology.w0 || C.wa != cosmology.wa || C.MGSigma != cosmology.MGSigma || C.MGmu != cosmology.MGmu || C.M_nu != cosmology.M_nu){return 1;}
if (C.Omega_m != cosmology.Omega_m || C.Omega_v != cosmology.Omega_v || C.Omega_nu != cosmology.Omega_nu || C.M_nu != cosmology.M_nu || C.h0 != cosmology.h0 || C.omb != cosmology.omb || C.n_spec != cosmology.n_spec|| C.alpha_s != cosmology.alpha_s || C.w0 != cosmology.w0 || C.wa != cosmology.wa || C.MGSigma != cosmology.MGSigma || C.MGmu != cosmology.MGmu || C.M_nu != cosmology.M_nu || C.MGalpha_K != cosmology.MGalpha_K ||C.MGalpha_B != cosmology.MGalpha_B ||C.MGalpha_T != cosmology.MGalpha_T ||C.MGalpha_M != cosmology.MGalpha_M){return 1;}
if (cosmology.A_s){
if(C.A_s != cosmology.A_s){return 1;}
}
Expand Down
15 changes: 14 additions & 1 deletion theory/structs.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,8 +69,13 @@ typedef struct {
double f_NL;
double MGSigma;
double MGmu;
//Horndeski params
double MGalpha_K;
double MGalpha_B;
double MGalpha_T;
double MGalpha_M;
}cosmopara;
cosmopara cosmology = {.A_s = 0., .alpha_s =0.0, .M_nu =0., .Omega_nu =0.,.coverH0= 2997.92458, .rho_crit = 7.4775e+21,.MGSigma=0.0,.MGmu=0.0};
cosmopara cosmology = {.A_s = 0., .alpha_s =0.0, .M_nu =0., .Omega_nu =0.,.coverH0= 2997.92458, .rho_crit = 7.4775e+21,.MGSigma=0.0,.MGmu=0.0,.MGalpha_K=0.0,.MGalpha_B=0.0,.MGalpha_T=0.0,.MGalpha_M=0.0};

typedef struct {
int shear_Nbin; // number of tomography bins
Expand Down Expand Up @@ -317,6 +322,10 @@ typedef struct input_cosmo_params_mpp {
double h0;
double MGSigma;
double MGmu;
double MGalpha_K;
double MGalpha_B;
double MGalpha_T;
double MGalpha_M;
} input_cosmo_params_mpp;

typedef struct input_cosmo_params {
Expand All @@ -329,6 +338,10 @@ typedef struct input_cosmo_params {
double h0;
double MGSigma;
double MGmu;
double MGalpha_K;
double MGalpha_B;
double MGalpha_T;
double MGalpha_M;
} input_cosmo_params;

typedef struct input_nuisance_params_mpp {
Expand Down