Also compute the Hubble time and lookback time in the cosmology module.

examples/CosmologicalBox/cosmo.yml

@@ -2,7 +2,7 @@
 InternalUnitSystem:
   UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
   UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
-  UnitVelocity_in_cgs: 1e0           # km/s in centimeters per second
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 

examples/main.c

@@ -488,8 +488,8 @@ int main(int argc, char *argv[]) {
 
     cosmology_update(&cosmo, &ee);
 
-    message("z=%e H(z)=%e w=%f t=%e [yrs]", cosmo.z, cosmo.H, cosmo.w,
-            cosmo.time / prog_const.const_year);
+    message("z=%e H(z)=%e w=%f t=%e [yrs] t_l=%e [yrs]", cosmo.z, cosmo.H, cosmo.w,
+            cosmo.time / prog_const.const_year, cosmo.lookback_time / prog_const.const_year);
   }
 
   return 0;

src/cosmology.c

@@ -139,8 +139,12 @@ void cosmology_update(struct cosmology *c, const struct engine *e) {
   /* H(z) */
   c->H = c->H0 * E_z;
 
-  /* Time since Big Bang */
+  /* Expansion rate */
+  c->a_dot = c->H * c->a;
+
+  /* Time */
   c->time = cosmology_get_time_since_big_bang(c, a);
+  c->lookback_time = c->universe_age_at_present_day - c->time;
 }
 
 /**
@@ -323,6 +327,11 @@ void cosmology_init_tables(struct cosmology *c) {
                       GSL_INTEG_GAUSS61, space, &result, &abserr);
   c->time_interp_table_offset = result;
 
+  /* Integrate the time \int_{0}^{1} dt */
+  gsl_integration_qag(&F, 0., 1, 0, 1.0e-13, GSL_workspace_size,
+                      GSL_INTEG_GAUSS61, space, &result, &abserr);
+  c->universe_age_at_present_day = result;
+
   /* Free the workspace and temp array */
   gsl_integration_workspace_free(space);
   free(a_table);
@@ -377,13 +386,16 @@ void cosmology_init(const struct swift_params *params,
   const double H0_cgs =
       100. * c->h * (km / (1.e6 * phys_const->const_parsec)); /* s^-1 */
   c->H0 = H0_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+  c->Hubble_time = 1. / c->H0;
 
   /* Set remaining variables to invalid values */
   c->H = -1.;
   c->a = -1.;
   c->a_inv = -1;
   c->z = -1.;
+  c->a_dot = -1.;
   c->time = -1.;
+  c->universe_age_at_present_day = -1.;
 
   /* Initialise the interpolation tables */
   c->drift_fac_interp_table = NULL;
@@ -403,8 +415,9 @@ void cosmology_print(const struct cosmology *c) {
       "Density parameters: [O_m, O_l, O_b, O_k, O_r] = [%f, %f, %f, %f, %f]",
       c->Omega_m, c->Omega_lambda, c->Omega_b, c->Omega_k, c->Omega_r);
   message("Dark energy equation of state: w_0=%f w_a=%f", c->w_0, c->w_a);
-  message("Hubble constant: h = %f, H_0 = %e U_t^(-1) (internal units)", c->h,
-          c->H0);
+  message("Hubble constant: h = %f, H_0 = %e U_t^(-1)", c->h, c->H0);
+  message("Hubble time: 1/H0 = %e U_t", c->Hubble_time);
+  message("Universe age at present day: %e U_t", c->universe_age_at_present_day);
 }
 
 /**

src/cosmology.h

@@ -46,9 +46,15 @@ struct cosmology {
   /*! Hubble constant at the current redshift (in internal units) */
   double H;
 
+  /*! Expansion rate at the current redshift (in internal units) */
+  double a_dot;
+
   /*! Time (in internal units) since the Big Bang */
   double time;
 
+  /*! Lookback time (in internal units) */
+  double lookback_time;
+
   /*! Dark-energy equation of state at the current time */
   double w;
 
@@ -66,6 +72,9 @@ struct cosmology {
   /*! Hubble constant at z = 0 (in internal units) */
   double H0;
 
+  /*! Hubble time 1/H0 */
+  double Hubble_time;
+
   /*! Matter density parameter */
   double Omega_m;
 
@@ -107,6 +116,9 @@ struct cosmology {
 
   /*! Time between Big Bang and first entry in the table */
   double time_interp_table_offset;
+
+  /*! Time at the present-day (a=1) */
+  double universe_age_at_present_day;
 };
 
 void cosmology_update(struct cosmology *c, const struct engine *e);

theory/Cosmology/coordinates.tex

 \subsection{Choice of co-moving coordinates}
 \label{ssec:ccordinates}
+
+Note that unlike the \gadget convention we do not express quantities
+with ``little h'' included. We hence use $\rm{Mpc}$ and not
+$\rm{Mpc}/h$ for instance. Similarly, the time integration operators
+(see below) also include an $h$-factor via the explicit appearance of
+the Hubble constant.\\
+
+
+MORE INFO HERE SOON