diff --git a/theory/cosmo2D_fourier.c b/theory/cosmo2D_fourier.c
index 2d50676d..205e6316 100644
--- a/theory/cosmo2D_fourier.c
+++ b/theory/cosmo2D_fourier.c
@@ -17,25 +17,29 @@ 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);
-  if(cosmology.MGSigma != 0.){
-    wkappa *= (1.+MG_Sigma(a));
-  }
+  #ifdef MODGRAV
+  //if you need modgrav functions put -DMODGRAV in your make file before the include statements.
+    double omegam=cosmology.Omega_m/(a*a*a);
+    double omegav=omv_vareos(a);
+    double hub = hoverh0(a);
+    hub = hub*hub;
+    if(cosmology.MGSigma != 0.){
+      wkappa *= (1.+MG_sigma(a, omegav, hub));
+    }
+    else if (cosmology.MGalpha_B != 0 || cosmology.MGalpha_T != 0 || cosmology.MGalpha_M != 0){
+      wkappa*= (1.+horndeski_sigma(a, omegav, hub, omegam));
+    }
+  #endif
+  #ifndef MODGRAV
+    if (cosmology.MGSigma != 0. || cosmology.MGalpha_K != 0 || cosmology.MGalpha_B != 0 || cosmology.MGalpha_T != 0 || cosmology.MGalpha_M != 0){
+      error(" Need to put -DMODGRAV in compiler directions.\n");
+    }
+  #endif
   return wkappa;
 }
 double W_gal(double a, double nz){
diff --git a/theory/cosmo3D.c b/theory/cosmo3D.c
index c82c12e2..ed92b58e 100644
--- a/theory/cosmo3D.c
+++ b/theory/cosmo3D.c
@@ -98,9 +98,20 @@ int func_for_growfac(double a,const double y[],double f[],void *params)
   double one_plus_mg_mu = 1.;
   hub = hub*hub;
   f[0]=y[1];
-  if(cosmology.MGmu != 0){
-    one_plus_mg_mu += cosmology.MGmu*omegav/hub/cosmology.Omega_v;
-  }
+  //if you need modgrav functions put -DMODGRAV in your make file before the include statements.
+  #ifdef MODGRAV
+    if(cosmology.MGmu != 0){
+      one_plus_mg_mu += MG_mu(omegav, hub);
+    }
+    else if (cosmology.MGalpha_B != 0 || cosmology.MGalpha_T != 0 || cosmology.MGalpha_M != 0){
+      one_plus_mg_mu += horndeski_mu(a, omegav, hub, omegam);
+    }
+  #endif
+  #ifndef MODGRAV
+    if (cosmology.MGmu != 0. || cosmology.MGalpha_K != 0|| cosmology.MGalpha_B != 0 || cosmology.MGalpha_T != 0 || cosmology.MGalpha_M != 0){
+      error(" Need to put -DMODGRAV in compiler directions.\n");
+    }
+  #endif
   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;
 }
diff --git a/theory/modgrav.c b/theory/modgrav.c
new file mode 100644
index 00000000..04ef9e07
--- /dev/null
+++ b/theory/modgrav.c
@@ -0,0 +1,83 @@
+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;
+  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);
+  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
+  //i.e. scale with Omega_DE(a)/Omega_DE(a=1)
+  //then, the alphas are equal to the constant values at today.
+  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 - H*H*((alpha_M -alpha_T)+alpha_B/(2.*(1.+alpha_T))));
+  //!!!!!!!!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 inline horndeski_alpha(double a, double fz){
+	//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;
+}
+
diff --git a/theory/parameters.c b/theory/parameters.c
index 6e249b52..30ada720 100644
--- a/theory/parameters.c
+++ b/theory/parameters.c
@@ -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.;
diff --git a/theory/recompute.c b/theory/recompute.c
index 2d7c9301..f4f7bd5b 100644
--- a/theory/recompute.c
+++ b/theory/recompute.c
@@ -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){
@@ -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;}
 }
 
@@ -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;}
   }
diff --git a/theory/structs.c b/theory/structs.c
index 1fc7be5e..71bafcb7 100644
--- a/theory/structs.c
+++ b/theory/structs.c
@@ -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
@@ -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 {
@@ -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 {