solved some todo

gravity_smg/background_smg.c

@@ -1,6 +1,6 @@
 /** @file background_smg.c Documented background_smg module
  *
- * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx
+ * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024
  *
  * Additional functions for the background module.
  * It contains all the hi_class related functions (_smg)
@@ -30,7 +30,6 @@
 * and phi is dimensionless (we can think of phi as given in units of the Planck mass
 * - A default module to numerically compute the derivatives when no analytic functions are given should be added.
 * Numerical derivatives may further serve as a consistency check.
-* TODO: Add a background_write_alpha_primes
 *
 * @param pba           Input: pointer to background structure
 * @param a             Input: scale factor
@@ -57,7 +56,6 @@ int background_gravity_functions_smg(
 	pvecback[pba->index_bg_p_tot_wo_smg] = *ptr_p_tot;
 
 	// scalar field + curvature not yet implemented
-	// TODO_EB: look for a better place to put this test
 	class_test(pba->K !=0 ,
 	     pba->error_message,
 	     "has_smg with curvature K = %e not yet implemented",pba->K);
@@ -333,7 +331,7 @@ int background_gravity_functions_smg(
 	*ptr_p_tot += pvecback[pba->index_bg_p_smg];
 	//divide relativistic & nonrelativistic (not very meaningful for oscillatory models)
 
-	//TODO: need to define menaingfully -> separate early universe (IC, BBN...) from late (Halofit...)
+	//TODO: need to define meaningfully -> separate early universe (IC, BBN...) from late (Halofit...)
 	//BUG: causes problem with halofit!, if not, causes bug with Brans-Dicke
 	*ptr_rho_de += pvecback[pba->index_bg_rho_smg];
 
@@ -731,7 +729,7 @@ int background_solve_smg(
 		pba->error_message,
 		pba->error_message);
 
-		// TODO_EB: note that the derivative is now calculated w.r.t. loga, while our _prime are w.r.t. tau
+		// NOTE: note that the derivative is now calculated w.r.t. loga, while our _prime are w.r.t. tau
 		double factor = pvecback[pba->index_bg_a]*pvecback[pba->index_bg_H];
 		double d_over_dtau;
 
@@ -967,7 +965,6 @@ int background_store_columntitles_smg(
     class_store_columntitle(titles,"slip_eff_smg",_TRUE_);
   }
 
-  //TODO: add in output background trigger
   if (pba->output_background_smg >= 3){
     class_store_columntitle(titles,"kineticity_prime_smg",_TRUE_);
     class_store_columntitle(titles,"braiding_prime_smg",_TRUE_);
@@ -1286,7 +1283,7 @@ int derivatives_alphas_smg(
 	/* necessary for calling array_interpolate(), but never used */
 	int last_index;
 
-	// TODO_EB: note that the derivative is now calculated w.r.t. loga, while our _prime are w.r.t. tau
+	// NOTE: note that the derivative is now calculated w.r.t. log(a), while our _prime are w.r.t. tau
 	double factor = pvecback[pba->index_bg_a]*pvecback[pba->index_bg_H];
 	double d_over_dtau;
 

gravity_smg/fourier_smg.c

@@ -1,6 +1,6 @@
 /** @file fourier_smg.c Documented fourier_smg module
  *
- * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx
+ * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024
  *
  * Additional functions for the fourier module.
  * It contains all the hi_class related functions (_smg)

gravity_smg/gravity_functions_smg.c

@@ -1,6 +1,6 @@
 /** @file gravity_functions_smg.c Documented gravity_functions_smg module
  *
- * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx
+ * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024
  *
  * This module contains all the complicated expressions
  * used in hi_class. In particular, they are casted in
@@ -21,7 +21,6 @@
 #include "gravity_functions_smg.h"
 
 
-// TODO_EB: reread all the documentation and check in particular Input/Output
 /**
 * Get gravity functions Es from Gs. These are the functions
 * entering the Friedmann time-time constraint, which is
@@ -56,7 +55,6 @@ int gravity_functions_Es_from_Gs_smg(
   double rho_tot = pvecback[pba->index_bg_rho_tot_wo_smg];
   double p_tot = pvecback[pba->index_bg_p_tot_wo_smg];
 
-  // TODO_EB: decide if it is better to keep it like this or add the pointers to the equations
   double G2=pgf->G2;
   double G2_X=pgf->G2_X, G2_XX=pgf->G2_XX, G2_XXX=pgf->G2_XXX;
   double G2_phi=pgf->G2_phi, G2_Xphi=pgf->G2_Xphi, G2_XXphi=pgf->G2_XXphi;
@@ -658,7 +656,7 @@ int gravity_functions_As_from_alphas_smg(
 	)*pow(H,-2);
 
 
-  // TODO_EB: remove below this when IC are recalculated
+  // NOTE: this could be written in terms of the A's and B's variables.
 
   pvecback[pba->index_bg_lambda_1_smg] = (run + (-1.)*ten)*(-3.)*bra + (1. + ten)*kin;
 
@@ -690,15 +688,17 @@ int gravity_functions_As_from_alphas_smg(
     + (1. + beh)*bra_p/a/H
     + (2. - bra)*beh_p/a/H;
 
-  // TODO_EB: check if there is a better alternative to regularizing these quantities
+  // NOTE: this is to regularize cs2 when both the numerator and denominator are
+  // below numerical precision
 	if (pvecback[pba->index_bg_cs2num_smg] == pvecback[pba->index_bg_kinetic_D_smg]) {
 		pvecback[pba->index_bg_cs2_smg] = 1.;
 	}
 	else {
 		pvecback[pba->index_bg_cs2_smg] = pvecback[pba->index_bg_cs2num_smg]/pvecback[pba->index_bg_kinetic_D_smg];
 	}
 
-  // TODO_EB: rewrite Geff and slip for beyond Horndeski (calculate them in the hi_class.nb Mathematica notebook)
+  // NOTE: Geff and slip are calculated for Horndeski. Consider extending them to
+  // beyond Horndeski
 	double beta_1 = (run + (-1.)*ten)*2. + (1. + ten)*bra;
 	double beta_2 = 2.*beta_1 + (2. + (-2.)*M2 + bra*M2)*(rho_tot + p_tot)*(-3.)*pow(H,-2)*pow(M2,-1) + ((-2.) + bra)*(rho_smg + p_smg)*(-3.)*pow(H,-2) + 2.*pow(H,-1)*bra_p*pow(a,-1);
 

gravity_smg/gravity_models_smg.c

@@ -641,7 +641,6 @@ int gravity_models_get_Gs_smg(
   }
 
   else if(pba->gravity_model_smg == galileon){
-    //TODO: change name with ccc_galileon
 
     double M3 = pow(pba->H0,2);
 
@@ -784,14 +783,14 @@ int gravity_models_get_back_par_smg(
 
     //Doran-Robbers model astro-ph/0601544
     //as implemented in Pettorino et al. 1301.5279
-    //TODO: check these expressions, they probably assume the standard evolution/friedmann eqs, etc...
-    //TODO: rewrite the expressions integrating the equation of state
+    //NOTE: rewrite the expressions integrating the equation of state
 
     double Om0 = pba->parameters_smg[0];
     double w0 = pba->parameters_smg[1];
+    // TODO_EB: add a precision parameter for 1e-10
+    //NOTE: I've regularized the expression adding a tiny Omega_e
     double Ome = pba->parameters_smg[2] + 1e-10;
 
-    //NOTE: I've regularized the expression adding a tiny Omega_e
     double Om = ((Om0 - Ome*(1.-pow(a,-3.*w0)))/(Om0 + (1.-Om0)*pow(a,3*w0)) + Ome*(1.-pow(a,-3*w0)));
     double dOm_da = (3*pow(a,-1 - 3*w0)*(-1 + Om0)*(-2*pow(a,3*w0)*(-1 + Om0)*Ome + Om0*Ome + pow(a,6*w0)*(Om0 - 2*Ome + Om0*Ome))*w0)/pow(-(pow(a,3*w0)*(-1 + Om0)) + Om0,2); //from Mathematica
     //I took a_eq = a*rho_r/rho_m, with rho_r = 3*p_tot_wo_smg
@@ -1128,7 +1127,7 @@ int gravity_models_initial_conditions_smg(
 		    double n = pba->parameters_smg[1];
 		    double Rshift0 = pba->parameters_smg[2];
 
-	      double H = sqrt(rho_rad); //TODO: for low n -> need to solve tracker + Constraint simultaneously. Increasing H (e.g. double H = 10*sqrt(rho_rad);) works for n~0.65.
+	      double H = sqrt(rho_rad); //NOTE: for low n -> need to solve tracker + Constraint simultaneously. Increasing H (e.g. double H = 10*sqrt(rho_rad);) works for n~0.65.
 	      double H0 = pba->H0;
 
 	      double signg = copysign(1.,g);
@@ -1187,8 +1186,7 @@ int gravity_models_initial_conditions_smg(
   /* expansion_smg when w is parametrized and rho obtained with integration */
 	if (pba->rho_evolution_smg == _TRUE_){
 
-		//DT: moved here from the definition of wowa_w model, to avoid pvecback[pba->index_bg_rho_smg] mix up
-    // TODO_EB: rethink this
+		//NOTE: moved here from the definition of wowa_w model, to avoid pvecback[pba->index_bg_rho_smg] mix up
 		if (pba->expansion_model_smg == wowa_w) {
 	    double Omega_const_smg = pba->parameters_smg[0];
 	    double w0 = pba->parameters_smg[1];

gravity_smg/include/hi_class.h

@@ -1,9 +1,7 @@
-// TODO_EB: write some overall description
-
 #ifndef __HI_CLASS__
 #define __HI_CLASS__
 
-#define _HI_CLASS_VERSION_ "TODO_EB:write version"
+#define _HI_CLASS_VERSION_ "3.0"
 
 #include "gravity_functions_smg.h"
 #include "gravity_models_smg.h"

gravity_smg/include/perturbations_smg.h

@@ -59,8 +59,6 @@ int perturbations_get_h_prime_ic_from_00_smg(struct background * pba,
                                              double eta,
                                              double delta_rho_tot);
 
-int perturbations_get_x_x_prime_newtonian_smg(struct perturbations_workspace * ppw);
-
 int perturbations_einstein_scalar_smg(struct precision * ppr,
                                       struct background * pba,
                                       struct thermodynamics * pth,

gravity_smg/input_smg.c

@@ -1,6 +1,6 @@
 /** @file input_smg.c Documented input_smg module
  *
- * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx
+ * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024
  *
  * Additional functions for the input module.
  * It contains all the hi_class related functions (_smg)
@@ -35,9 +35,6 @@ int input_warnings_smg(
                        int input_verbose
                        ) {
 
-  // TODO_EB: rethink this test. The default should have no warnings, and
-  // we have to be sure that get_h_from_trace == _TRUE_ is working properly
-
   /* Here we put a warning as we want to encourage hi_class users to get
   h_prime from the trace of the Einstein ij equation than from the Einstein 00
   equation. This is because the Einstein 00 equation has a gauge dependent

gravity_smg/perturbations_smg.c

@@ -1,6 +1,6 @@
 /** @file perturbations_smg.c Documented perturbations_smg module
  *
- * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx
+ * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024
  *
  * Additional functions for the perturbation module.
  * It contains all the hi_class related functions (_smg)
@@ -364,14 +364,14 @@ int perturbations_prepare_k_output_smg(
                                        ) {
 
   if (ppt->use_pert_var_deltaphi_smg==_TRUE_) {
-    class_store_columntitle(ppt->scalar_titles, "delta_phi_smg", _TRUE_);
-    class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_smg", _TRUE_);
-    class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_prime_smg", _TRUE_);
+    class_store_columntitle(ppt->scalar_titles, "delta_phi_smg (sync)", _TRUE_);
+    class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_smg (sync)", _TRUE_);
+    class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_prime_smg (sync)", _TRUE_);
   }
   else {
-    class_store_columntitle(ppt->scalar_titles, "V_x_smg", _TRUE_);
-    class_store_columntitle(ppt->scalar_titles, "V_x_prime_smg", _TRUE_);
-    class_store_columntitle(ppt->scalar_titles, "V_x_prime_prime_smg", _TRUE_);
+    class_store_columntitle(ppt->scalar_titles, "V_x_smg (sync)", _TRUE_);
+    class_store_columntitle(ppt->scalar_titles, "V_x_prime_smg (sync)", _TRUE_);
+    class_store_columntitle(ppt->scalar_titles, "V_x_prime_prime_smg (sync)", _TRUE_);
   }
 
   /* Quasi-static functions smg*/
@@ -486,35 +486,13 @@ int perturbations_get_h_prime_ic_from_00_smg(
   double p_smg = ppw->pvecback[pba->index_bg_p_smg];
   double kin = ppw->pvecback[pba->index_bg_kineticity_smg];
   double bra = ppw->pvecback[pba->index_bg_braiding_smg];
-  /* TODO_EB: rewrite this equation with new variables */
+  /* NOTE: this could be rewritten this equation with A and B variables */
   ppw->pv->y[ppw->pv->index_pt_h_prime_from_trace] = (-4. * pow(H, -1) * pow(k, 2) * eta / a - 6. * pow(H, -1) * pow(M2, -1) * delta_rho_tot * a + 2. * H * (3. * bra + kin) * ppw->pv->y[ppw->pv->index_pt_x_prime_smg] * a + (2. * bra * pow(k, 2) + (-18. + 15. * bra + 2. * kin) * rho_smg * pow(a, 2) + (-18. * DelM2 + 15. * bra * M2 + 2. * kin * M2) * rho_tot * pow(M2, -1) * pow(a, 2) + (-2. * DelM2 + bra * M2) * 9. * pow(M2, -1) * p_tot * pow(a, 2) + 9. * (-2. + bra) * p_smg * pow(a, 2)) * ppw->pv->y[ppw->pv->index_pt_x_smg]) * pow(-2. + bra, -1);
 
   return _SUCCESS_;
 }
 
 
-/**
- * Transform synchronous gauge scalar field to newtonian.
- *
- * @param ppw              Input/Output: pointer to perturbation workspace structure
- * @return the error status
- */
-int perturbations_get_x_x_prime_newtonian_smg(
-                                              struct perturbations_workspace * ppw
-                                              ) {
-
-  int qs_array_smg[] = _VALUES_QS_SMG_FLAGS_;
-
-  /* scalar field: TODO_EB: add gauge transformations (when we add Newtonian gauge) */
-  if (qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0) {
-    ppw->pv->y[ppw->pv->index_pt_x_smg] += 0.;
-    ppw->pv->y[ppw->pv->index_pt_x_prime_smg] += 0.;
-  }
-
-  return _SUCCESS_;
-}
-
-
 /**
  * Scalar Einstein equations.
  *
@@ -2027,7 +2005,7 @@ int perturbations_adiabatic_ic_smg(
 
   a_prime_over_a = ppw->pvecback[pba->index_bg_H]*a;
 
-  H = ppw->pvecback[pba->index_bg_H];//TODO_EB
+  H = ppw->pvecback[pba->index_bg_H];
   Hprime = ppw->pvecback[pba->index_bg_H_prime];
   a = ppw->pvecback[pba->index_bg_a];
   rho_tot = ppw->pvecback[pba->index_bg_rho_tot_wo_smg];
@@ -2068,7 +2046,7 @@ int perturbations_adiabatic_ic_smg(
   DelM2 = ppw->pvecback[pba->index_bg_delta_M2_smg];//M2-1
 
 
-  /* TODO_EB: revisit initial conditions for beyond horndeski and oscillations */
+  /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */
 
   if (qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0) {
 
@@ -2412,7 +2390,7 @@ int perturbations_adiabatic_ic_smg(
         shear_ur =  A_sigma_nu_smg/A1_eta_smg* ktau_two * ppr->curvature_ini;
         // /(45.+12.*fracnu) * (3.*s2_squared-1.) * (1.+(4.*fracnu-5.)/4./(2.*fracnu+15.)*tau*om) * ppr->curvature_ini;//TBC /s2_squared; /* shear of ultra-relativistic neutrinos/relics */  //TBC:0
 
-        //TODO: needs to be modified?
+        //TODO: needs to be modified? (ask ILS)
         l3_ur = ktau_three/A1_eta_smg*2./7./(12.*fracnu+45.)* ppr->curvature_ini;//ILS
 
         if(ppt->perturbations_verbose > 8)
@@ -2633,7 +2611,7 @@ int perturbations_isocurvature_cdm_ic_smg(
 
   //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static
   if((ppt->pert_initial_conditions_smg==ext_field_attr) && (qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) {
-    /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */
+    /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */
 
     double coeff_isocurv_smg;
 
@@ -2718,7 +2696,7 @@ int perturbations_isocurvature_b_ic_smg(
 
   //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static
   if((ppt->pert_initial_conditions_smg==ext_field_attr)&&(qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) {
-    /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */
+    /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */
     double coeff_isocurv_smg;
 
     int nexpo= 1;
@@ -2808,7 +2786,7 @@ int perturbations_isocurvature_urd_ic_smg(
   //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static
   if((ppt->pert_initial_conditions_smg==ext_field_attr)&&(qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) {
 
-    /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */
+    /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */
     double coeff_isocurv_smg;
 
     double a = ppw->pvecback[pba->index_bg_a];
@@ -2913,7 +2891,7 @@ int perturbations_isocurvature_urv_ic_smg(
 
   //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static
   if((pba->has_smg == _TRUE_)&&(ppt->pert_initial_conditions_smg==ext_field_attr)&&(qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) {
-    /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */
+    /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */
 
     double coeff_isocurv_smg;
 

source/background.c

@@ -2403,7 +2403,7 @@ int background_initial_conditions(
 
   /* Just checking that our initial time indeed is deep enough in the radiation
      dominated regime (_smg) */
-  // TODO_EB: rethink this test (r+(de?)). There was a class_test_except with free(pvecback);free(pvecback_integration);background_free(pba)
+  // TODO_EB: There was a class_test_except with free(pvecback);free(pvecback_integration);background_free(pba)
   class_test(fabs(pvecback[pba->index_bg_Omega_r]+pvecback[pba->index_bg_Omega_de]-1.) > ppr->tol_initial_Omega_r,
              pba->error_message,
              "Omega_r = %e, Omega_de = %e, not close enough to 1. Decrease a_ini_over_a_today_default in order to start from radiation domination.",

source/input.c

@@ -555,7 +555,6 @@ int input_shooting(struct file_content * pfc,
   /* for each target, module up to which we need to run CLASS in order
      to compute the targetted quantities (not running the whole code
      each time to saves a lot of time) */
-  // TODO_EB: here we didn't add any computation stage in the previous version of hi_class, I think we should. Check it
   enum computation_stage target_cs[] = {cs_thermodynamics, /* computation stage for target '100*theta_s' */
                                         cs_thermodynamics, /* computation stage for target 'theta_s_100' */
                                         cs_background,     /* computation stage for target 'Omega_dcdmdr' */

source/perturbations.c

@@ -6060,14 +6060,6 @@ int perturbations_initial_conditions(struct precision * ppr,
            +ppw->pvecback[pba->index_bg_phi_prime_scf]*alpha_prime);
       }
 
-      if (pba->has_smg == _TRUE_) {
-        class_call(
-          perturbations_get_x_x_prime_newtonian_smg(ppw),
-          ppt->error_message,
-          ppt->error_message
-        );
-      }
-
       if ((pba->has_ur == _TRUE_) || (pba->has_ncdm == _TRUE_) || (pba->has_dr == _TRUE_)  || (pba->has_idr == _TRUE_)) {
 
         delta_ur -= 4.*a_prime_over_a*alpha;