diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
index 16c2fc55777eafb367896b4454f2c716538e6058..23387bcf91b294230ef4ce8b01d117e706e0f7af 100644
--- a/examples/CoolingBox/analytical_test.py
+++ b/examples/CoolingBox/analytical_test.py
@@ -136,15 +136,15 @@ log_nh = float(sys.argv[2])
 redshift = float(sys.argv[1])
 p = plt.figure(figsize = (6,6))
 p1, = plt.loglog(t_sol,u_sol,linewidth = 0.5,color = 'k', marker = '*', markersize = 1.5,label = 'explicit ODE')
-p2, = plt.loglog(t_snapshots_cgs,mean_u,linewidth = 0.5,color = 'r', marker = '*', markersize = 1.5,label = 'swift mean alexei')
+p2, = plt.loglog(t_snapshots_cgs,mean_u,linewidth = 0.5,color = 'r', marker = '*', markersize = 1.5,label = 'swift mean')
 l = legend(handles = [p1,p2])
 xlabel("${\\rm{Time~[s]}}$", labelpad=0, fontsize = 14)
 ylabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
 if int(sys.argv[4]) == 1:
-	title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, solar metallicity,\n relative error alexei: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
+	title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, solar metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
 	name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_solar.png"
 elif int(sys.argv[4]) == 0:
-	title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, zero metallicity,\n relative error alexei: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
+	title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, zero metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
 	name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_zero_metal.png"
 
 savefig(name, dpi=200)
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index 03f5626023bf1570caa65df52673c8b5565ec4a3..eb76f3d7b5b7a100fdcaf9b38b5d4c82f2f20104 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/examples/CoolingBox/makeIC.py
@@ -27,8 +27,8 @@ from numpy import *
 # Parameters
 periodic= 1           # 1 For periodic box
 boxSize = 1           # 1 kiloparsec    
-rho = 34504.4797588         # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
-P = 44554455.4455             # Pressure in code units (at 10^6K)
+rho = 3450447.97588         # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
+P = 4455445544.55             # Pressure in code units (at 10^6K)
 gamma = 5./3.         # Gas adiabatic index
 eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "coolingBox.hdf5" 
diff --git a/examples/CoolingRates/testCooling.c b/examples/CoolingRates/testCooling.c
index 4f141097a99d1b73e4eac4c6ffb3f3098cebdd9f..22ca23423ebc87cf8ca5a9ba0b08762641df13ab 100644
--- a/examples/CoolingRates/testCooling.c
+++ b/examples/CoolingRates/testCooling.c
@@ -69,299 +69,6 @@ void set_quantities(struct part *restrict p,
   if(cosmo->z >= 1) hydro_set_init_internal_energy(p,(u*pow(scale_factor,2))/cooling->internal_energy_scale);
 }
 
-/*
- * @brief Construct 1d tables from 4d EAGLE tables by 
- * interpolating over redshift, hydrogen number density
- * and helium fraction. 
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param internal_const Physical constants data structure
- * @param temp_table Pointer to 1d interpolated table of temperature values
- * @param H_plus_He_heat_table Pointer to 1d interpolated table of cooling rates
- * due to hydrogen and helium
- * @param H_plus_He_electron_abundance_table Pointer to 1d interpolated table
- * of electron abundances due to hydrogen and helium
- * @param element_electron_abundance_table Pointer to 1d interpolated table 
- * of electron abundances due to metals
- * @param element_print_cooling_table Pointer to 1d interpolated table of
- * cooling rates due to each of the metals
- * @param element_cooling_table Pointer to 1d interpolated table of cooling
- * rates due to the contribution of all the metals
- * @param abundance_ratio Pointer to array of ratios of metal abundances to solar
- * @param ub Upper bound in temperature on table construction 
- * @param lb Lower bound in temperature on table construction
- */
-void construct_1d_tables_test(const struct part *restrict p,
-			 const struct cooling_function_data *restrict cooling,
-			 const struct cosmology *restrict cosmo,
-			 const struct phys_const *restrict internal_const,
-			 double *temp_table,
-			 double *H_plus_He_heat_table,
-			 double *H_plus_He_electron_abundance_table,
-			 double *element_electron_abundance_table,
-			 double *element_print_cooling_table,
-			 double *element_cooling_table,
-			 float *abundance_ratio,
-			 float *ub, float *lb){
-
-  // obtain mass fractions and number density for particle
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                             internal_const) *
-                cooling->number_density_scale;
-
-  // find redshift, helium fraction, hydrogen number density indices and offsets
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-  if (cosmo->z > cooling->reionisation_redshift) { 
-    // Photodissociation table 
-    printf("Eagle testCooling.c photodissociation table redshift reionisation redshift %.5e %.5e\n", cosmo->z, cooling->reionisation_redshift);
-    construct_1d_table_from_3d(p, cooling, cosmo, internal_const, 
-                cooling->table.photodissociation_cooling.temperature, 
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb); 
-    construct_1d_table_from_3d( 
-                p, cooling, cosmo, internal_const, 
-                cooling->table.photodissociation_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,   
-                He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb); 
-    construct_1d_table_from_3d( 
-                p, cooling, cosmo, internal_const, 
-                cooling->table.photodissociation_cooling.H_plus_He_electron_abundance, 
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, 
-                cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-    construct_1d_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.photodissociation_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_print_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.photodissociation_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_table_from_2d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.photodissociation_cooling.electron_abundance,
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-  } else if (cosmo->z > cooling->Redshifts[cooling->N_Redshifts - 1]) {
-    printf("Eagle testCooling.c no compton table redshift max redshift %.5e %.5e\n", cosmo->z, cooling->Redshifts[cooling->N_Redshifts - 1]);
-    // High redshift table
-    construct_1d_table_from_3d(p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.temperature,
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
-    construct_1d_table_from_3d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,  
-                He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
-    construct_1d_table_from_3d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.H_plus_He_electron_abundance,
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He,
-                cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-    construct_1d_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_print_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_table_from_2d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.electron_abundance,
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-  } else {
-    // Normal tables 
-    printf("Eagle testCooling.c normal table redshift %.5e\n", cosmo->z);
-    construct_1d_table_from_4d(p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.temperature,
-                z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-                cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
-    construct_1d_table_from_4d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-                cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
-    construct_1d_table_from_4d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
-                dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
-                cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-    construct_1d_table_from_4d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_print_table_from_4d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_table_from_3d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-  } 
-    
-}
-
-/*
- * @brief Compare calculation of dlambda/du (gradient of cooling rate, 
- * needed for Newton's method) between interpolating 1d tables and 
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_dlambda_du(
-  const struct unit_system *restrict us,
-  struct part *restrict p,
-  const struct xpart *restrict xp,
-  const struct phys_const *restrict internal_const,
-  const struct cooling_function_data *restrict cooling,
-  const struct cosmology *restrict cosmo){
-
-  // allocate tables
-  double H_plus_He_heat_table[cooling->N_Temp];              
-  double H_plus_He_electron_abundance_table[cooling->N_Temp];
-  double temp_table[cooling->N_Temp];  
-  double element_cooling_table[cooling->N_Temp];         
-  double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];  
-  double element_electron_abundance_table[cooling->N_Temp];
-  double rate_element_table[cooling->N_Elements+2];
-  float *abundance_ratio, cooling_du_dt1, cooling_du_dt2;
-  double dlambda_du1, dlambda_du2;
-
-  // calculate ratio of particle metal abundances to solar 
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(p, cooling, abundance_ratio);
-  
-  // set hydrogen number density and internal energy
-  float nh = 1.0e-1, u = 1.0e9;
-  set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-  
-  // extract hydrogen and helium mass fractions
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                           internal_const) *
-              cooling->number_density_scale;
-  
-  // find redshift, hydrogen number density and helium fraction indices and offsets
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-
-  // construct tables
-  float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
-  float lower_bound = cooling->Temp[0]/eagle_log_10_e;
-  construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table, 
-                           H_plus_He_heat_table, H_plus_He_electron_abundance_table, 
-			   element_electron_abundance_table, element_print_cooling_table, 
-			   element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-  
-  // calculate dlambda/du for different values of internal energy
-  int nt = 10;
-  for(int i = 0; i < nt; i++){
-    u = pow(10.0,9 + i);
-    set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-    cooling_du_dt1 = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du1,H_plus_He_heat_table,
-                      H_plus_He_electron_abundance_table,
-       	              element_print_cooling_table,
-       	              element_electron_abundance_table,
-       	              temp_table,
-       	              p,cooling,cosmo,internal_const,rate_element_table);
-    cooling_du_dt2 = eagle_metal_cooling_rate(log10(u),
-       	              &dlambda_du2,z_index, dz, n_h_i, d_n_h, He_i, d_He,
-       	              p,cooling,cosmo,internal_const,NULL,
-       	              abundance_ratio);
-    printf("u du_dt_1d du_dt_4d dlambda_du_1d dlambda_du_4d %.5e %.5e %.5e %.5e %.5e\n",u,cooling_du_dt1,cooling_du_dt2,dlambda_du1,dlambda_du2);
-  }
-}
-
-/*
- * @brief Compare calculating temperature between interpolating 1d tables and 
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_temp(
-  const struct unit_system *restrict us,
-  struct part *restrict p,
-  const struct xpart *restrict xp,
-  const struct phys_const *restrict internal_const,
-  const struct cooling_function_data *restrict cooling,
-  const struct cosmology *restrict cosmo){
-
-  // allocate tables
-  double H_plus_He_heat_table[cooling->N_Temp];              
-  double H_plus_He_electron_abundance_table[cooling->N_Temp];
-  double temp_table[cooling->N_Temp];  
-  double element_cooling_table[cooling->N_Temp];         
-  double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];  
-  double element_electron_abundance_table[cooling->N_Temp];
-  float *abundance_ratio;
-
-  // calculate ratio of particle metal abundances to solar 
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(p, cooling, abundance_ratio);
-  
-  // set hydrogen number density and internal energy
-  float nh = 1.0e-1, u = 1.0e9;
-  set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
-  // extract hydrogen and helium mass fractions
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                           internal_const) *
-              cooling->number_density_scale;
-  
-  // find redshift, hydrogen number density and helium fraction indices and offsets
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-  // construct tables
-  float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
-  float lower_bound = cooling->Temp[0]/eagle_log_10_e;
-  construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table, 
-                           H_plus_He_heat_table, H_plus_He_electron_abundance_table, 
-			   element_electron_abundance_table, element_print_cooling_table, 
-			   element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-  
-  // calculate temperature for different values of internal energy
-  float T1d, T4d, delta_u;
-  int nt = 10;
-  for(int i = 0; i < nt; i++){
-    u = pow(10.0,9 + i);
-    set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-    T1d = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                   cooling);
-    T4d = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-        		 dz, d_n_h, d_He,cooling, cosmo);
-    printf("u T1d T4d %.5e %.5e %.5e\n",u,pow(10.0,T1d),pow(10.0,T4d));
-  }
-}
-
 /*
  * @brief Produces contributions to cooling rates for different 
  * hydrogen number densities, from different metals, 
@@ -447,14 +154,6 @@ int main(int argc, char **argv) {
   abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
   abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
   
-  // Declare 1D tables 
-  double H_plus_He_heat_table[cooling.N_Temp];              
-  double H_plus_He_electron_abundance_table[cooling.N_Temp];
-  double temp_table[cooling.N_Temp];  
-  double element_cooling_table[cooling.N_Temp];         
-  double element_print_cooling_table[cooling.N_Elements * cooling.N_Temp];  
-  double element_electron_abundance_table[cooling.N_Temp];
-  
   // extract mass fractions, calculate table indices and offsets 
   float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
   float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
@@ -464,139 +163,48 @@ int main(int argc, char **argv) {
   get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
   get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
 
-  float upper_bound = cooling.Temp[cooling.N_Temp-1]/eagle_log_10_e;
-  float lower_bound = cooling.Temp[0]/eagle_log_10_e;
-
-  int nt = 250;//, n_nh = 6;
+  int nt = 250;
   double u = pow(10.0,10);
-  //if (argc == 1 || strcmp(argv[1], "nh") == 0){
-  // Calculate cooling rates at different densities 
-    //for(int i = 0; i < n_nh; i++){
-    //  // Open files 
-    //  char output_filename[21];
-    //  sprintf(output_filename, "%s%d%s", "cooling_output_", i, ".dat");
-    //  FILE *output_file = fopen(output_filename, "w");
-    //  if (output_file == NULL) {
-    //    printf("Error opening file!\n");
-    //    exit(1);
-    //  }
-
-    //  // set hydrogen number density, construct tables
-    //  nh = pow(10.0,-i);
-    //  set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
-    //  construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const, 
-    //    			temp_table, H_plus_He_heat_table, 
-    //      		H_plus_He_electron_abundance_table, 
-    //      		element_electron_abundance_table, 
-    //      		element_print_cooling_table, 
-    //      		element_cooling_table, 
-    //      		abundance_ratio, &upper_bound, &lower_bound);
-    //  get_index_1d(cooling.nH, cooling.N_nH, log10(nh), &n_h_i, &d_n_h);
-
-    //  for(int j = 0; j < nt; j++){
-    //    // set internal energy
-    //    set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,11.0 + j*8.0/nt));
-    //    u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
-    //    double dlambda_du;
-    //    float delta_u, cooling_du_dt, logT;
-
-    //    // calculate cooling rates using 1d tables
-    //    if (tables == 1) {
-    //      cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du,H_plus_He_heat_table,
-    //                              H_plus_He_electron_abundance_table,
-    //         	              element_print_cooling_table,
-    //         	              element_electron_abundance_table,
-    //         	              temp_table,
-    //         	              z_index, dz, n_h_i, d_n_h, He_i, d_He,
-    //         	              &p,&cooling,&cosmo,&internal_const,NULL,
-    //         	              abundance_ratio);
-    //      logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-    //                     &p,&cooling, &cosmo, &internal_const);
-    //    } else {
-    //    // calculate cooling rates using 4d tables
-    //      cooling_du_dt = eagle_metal_cooling_rate(log10(u),
-    //         	              &dlambda_du,z_index, dz, n_h_i, d_n_h, He_i, d_He,
-    //         	              &p,&cooling,&cosmo,&internal_const,NULL,
-    //         	              abundance_ratio);
-    //      logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-    //          		 dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
-    //    }
-    //    float temperature_swift = pow(10.0,logT);
 
-    //    fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
-    //  }
-    //  fclose(output_file);
-    //}
-  //}
   // Calculate contributions from metals to cooling rate
-  //if (argc >= 1 || strcmp(argv[1],"metals") == 0) {
-    // open file
-    char output_filename[21];
-    sprintf(output_filename, "%s", "cooling_output.dat");
-    FILE *output_file = fopen(output_filename, "w");
-    if (output_file == NULL) {
-      printf("Error opening file!\n");
-      exit(1);
-    }
-
-    // set hydrogen number density, construct 1d tables
-    if (log_10_nh == 100) {
-      nh = 1.0e-1;
-    } else {
-      nh = pow(10.0,log_10_nh);
-    }
-    //if (argc > 2) nh = pow(10.0,strtod(argv[2],NULL));
-    printf("Eagle testcooling.c nh %.5e\n", nh);
-    u = pow(10.0,14.0);
-    set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
-    construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const, 
-      			temp_table, H_plus_He_heat_table, 
-        		H_plus_He_electron_abundance_table, 
-        		element_electron_abundance_table, 
-        		element_print_cooling_table, 
-        		element_cooling_table, 
-        		abundance_ratio, &upper_bound, &lower_bound);
-    float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
-                                             &internal_const) *
-                cooling.number_density_scale;
-    printf("Eagle testcooling.c inn_h %.5e\n", inn_h);
-    get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-    // Loop over internal energy
-    for(int j = 0; j < nt; j++){
-      set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
-      u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
-      float delta_u, cooling_du_dt, logT;
-
-      // calculate cooling rates using 1d tables
-      if (tables == 1) {
-        cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
-                              H_plus_He_electron_abundance_table,
-           	              element_print_cooling_table,
-           	              element_electron_abundance_table,
-           	              temp_table,
-           	              &p,&cooling,&cosmo,&internal_const);
-        logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                       &cooling);
-      } else {
-      // calculate cooling rates using 4d tables
-        cooling_du_dt = eagle_print_metal_cooling_rate(
-			      z_index, dz, n_h_i, d_n_h, He_i, d_He,
-           	              &p,&cooling,&cosmo,&internal_const,
-           	              abundance_ratio);
-          logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-	      		 dz, d_n_h, d_He,&cooling, &cosmo);
-      }
-      //float temperature_swift = pow(10.0,logT);
-      //fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
-      fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
-    }
-    fclose(output_file);
-  //}
+  // open file
+  char output_filename[21];
+  sprintf(output_filename, "%s", "cooling_output.dat");
+  FILE *output_file = fopen(output_filename, "w");
+  if (output_file == NULL) {
+    printf("Error opening file!\n");
+    exit(1);
+  }
 
-  // compare temperatures and dlambda/du calculated from 1d and 4d tables
-  //compare_temp(&us, &p, &xp, &internal_const, &cooling, &cosmo);
-  //compare_dlambda_du(&us, &p, &xp, &internal_const, &cooling, &cosmo);
+  // set hydrogen number density, construct 1d tables
+  if (log_10_nh == 100) {
+    nh = 1.0e-1;
+  } else {
+    nh = pow(10.0,log_10_nh);
+  }
+  //if (argc > 2) nh = pow(10.0,strtod(argv[2],NULL));
+  u = pow(10.0,14.0);
+  set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
+  float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
+                                           &internal_const) *
+              cooling.number_density_scale;
+  get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
+
+  // Loop over internal energy
+  for(int j = 0; j < nt; j++){
+    set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
+    u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
+    float cooling_du_dt;
+
+    // calculate cooling rates using 4d tables
+    cooling_du_dt = eagle_print_metal_cooling_rate(
+    		      z_index, dz, n_h_i, d_n_h, He_i, d_He,
+       	              &p,&cooling,&cosmo,&internal_const,
+       	              abundance_ratio);
+    fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
+  }
+  fclose(output_file);
+  printf("done\n");
 
   free(params);
   return 0;
diff --git a/examples/CoolingRates/testCooling.yml b/examples/CoolingRates/testCooling.yml
index db3a7cc00a529dcc403d4c24c11b30186baa7896..3010636433d35421009539d1858722ea89820bc9 100644
--- a/examples/CoolingRates/testCooling.yml
+++ b/examples/CoolingRates/testCooling.yml
@@ -59,16 +59,16 @@ InitialConditions:
   cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
 
 EAGLEChemistry:
-  InitMetallicity: 0.0
-  InitAbundance_Hydrogen: 0.752
-  InitAbundance_Helium: 0.248
-  InitAbundance_Carbon: 0.000
-  InitAbundance_Nitrogen: 0.000
-  InitAbundance_Oxygen: 0.000
-  InitAbundance_Neon: 0.000
-  InitAbundance_Magnesium: 0.000
-  InitAbundance_Silicon: 0.000
-  InitAbundance_Iron: 0.000
+  InitMetallicity: 0.014
+  InitAbundance_Hydrogen: 0.70649785
+  InitAbundance_Helium: 0.28055534
+  InitAbundance_Carbon: 2.0665436e-3
+  InitAbundance_Nitrogen: 8.3562563e-4
+  InitAbundance_Oxygen: 5.4926244e-3
+  InitAbundance_Neon: 1.4144605e-3
+  InitAbundance_Magnesium: 5.907064e-4
+  InitAbundance_Silicon: 6.825874e-4
+  InitAbundance_Iron: 1.1032152e-3
   CalciumOverSilicon: 0.0941736
   SulphurOverSilicon: 0.6054160
 
diff --git a/examples/CoolingTest_Alexei/Makefile.am b/examples/CoolingTest_Alexei/Makefile.am
deleted file mode 100644
index e878a3792a84f7ace5f3b60fcc5875169145731b..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/Makefile.am
+++ /dev/null
@@ -1,32 +0,0 @@
-# tHIS FIle is part of SWIFT.
-# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-# Add the source directory and the non-standard paths to the included library headers to CFLAGS
-AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS)
-
-AM_LDFLAGS = $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS) $(PROFILER_LIBS)
-
-# Extra libraries.
-EXTRA_LIBS = $(HDF5_LIBS)
-
-# Programs.
-bin_PROGRAMS = testCooling
-
-# Sources
-testCooling_SOURCES = testCooling.c
-testCooling_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
-testCooling_LDADD =  ../../src/.libs/libswiftsim.a $(EXTRA_LIBS)
-
diff --git a/examples/CoolingTest_Alexei/plots.py b/examples/CoolingTest_Alexei/plots.py
deleted file mode 100644
index 969e6de713a4c421cc14d5f88ef68bafaab59a41..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/plots.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import matplotlib.pyplot as plt
-
-elements = 11
-temperature = []
-cooling_rate = [[] for i in range(elements+1)]
-u = []
-fu = []
-length = 0
-file_in = open('cooling_output.dat', 'r')
-for line in file_in:
-	data = line.split()
-	temperature.append(float(data[0]))
-	cooling_rate[0].append(-float(data[1]))
-
-file_in.close()
-
-for elem in range(elements):
-	file_in = open('cooling_element_'+str(elem)+'.dat','r')
-	for line in file_in:
-        	data = line.split()
-        	cooling_rate[elem+1].append(-float(data[0]))
-	file_in.close()
-
-#file_in = open('newton_output.dat', 'r')
-#for line in file_in:
-#	data = line.split()
-#	u.append(float(data[0]))
-#	fu.append(float(data[1]))
-#
-#file_in.close()
-
-#p0, = plt.plot(u,fu, color = 'k')
-#p1, = plt.plot(u,[0 for x in u], color = 'r')
-#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
-#plt.xlim([1e13,2.0e14])
-#plt.ylim([0e13,2e14])
-#plt.xlabel('u')
-#plt.ylabel('f(u)')
-#plt.show()
-
-p0, = plt.loglog(temperature, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
-p2, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
-p3, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
-p4, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
-p5, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
-p6, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
-p7, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
-p8, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
-p9, = plt.loglog(temperature, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
-p10, = plt.loglog(temperature, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
-#plt.ylim([1e-24,1e-21])
-#plt.xlim([1e2,1e8])
-plt.xlabel('Temperature (K)')
-plt.ylabel('Cooling rate (eagle units...)')
-plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
-#plt.legend(handles = [p0,p2,p3,p4,p5,p6,p7,p8,p9])
-plt.show()
diff --git a/examples/CoolingTest_Alexei/plots_nh.py b/examples/CoolingTest_Alexei/plots_nh.py
deleted file mode 100644
index ab3300750e2497ec8026dbe2d5c6c7d358f87831..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/plots_nh.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import matplotlib.pyplot as plt
-
-elements = 6
-temperature = [[] for i in range(elements+1)]
-cooling_rate = [[] for i in range(elements+1)]
-u = []
-fu = []
-length = 0
-
-for elem in range(elements):
-	file_in = open('cooling_output_'+str(elem)+'.dat','r')
-	for line in file_in:
-		data = line.split()
-		temperature[elem+1].append(float(data[0]))
-		cooling_rate[elem+1].append(-float(data[1]))
-	file_in.close()
-
-#file_in = open('newton_output.dat', 'r')
-#for line in file_in:
-#	data = line.split()
-#	u.append(float(data[0]))
-#	fu.append(float(data[1]))
-#
-#file_in.close()
-
-#p0, = plt.plot(u,fu, color = 'k')
-#p1, = plt.plot(u,[0 for x in u], color = 'r')
-#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
-#plt.xlim([1e13,2.0e14])
-#plt.ylim([0e13,2e14])
-#plt.xlabel('u')
-#plt.ylabel('f(u)')
-#plt.show()
-
-#p0, = plt.loglog(temperature[0], cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1, = plt.loglog(temperature[1], cooling_rate[1], linewidth = 0.5, color = 'k', label = 'nh = 10^0')
-p2, = plt.loglog(temperature[2], cooling_rate[2], linewidth = 0.5, color = 'b', label = 'nh = 10^-1')
-p3, = plt.loglog(temperature[3], cooling_rate[3], linewidth = 0.5, color = 'g', label = 'nh = 10^-2')
-p4, = plt.loglog(temperature[4], cooling_rate[4], linewidth = 0.5, color = 'r', label = 'nh = 10^-3')
-p5, = plt.loglog(temperature[5], cooling_rate[5], linewidth = 0.5, color = 'c', label = 'nh = 10^-4')
-p6, = plt.loglog(temperature[6], cooling_rate[6], linewidth = 0.5, color = 'm', label = 'nh = 10^-5')
-plt.ylim([1e-24,1e-21])
-plt.xlim([1e4,1e8])
-plt.legend(handles = [p1,p2,p3,p4,p5,p6])
-plt.xlabel('Temperature (K)')
-plt.ylabel('Cooling rate (eagle units...)')
-plt.show()
diff --git a/examples/CoolingTest_Alexei/testCooling.c b/examples/CoolingTest_Alexei/testCooling.c
deleted file mode 100644
index 3a8cc301812097baf1d9e136887de86000097198..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling.c
+++ /dev/null
@@ -1,420 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-  
-enum table_index {
-  EAGLE_Hydrogen = 0,
-  EAGLE_Helium,
-  EAGLE_Carbon,
-  EAGLE_Nitrogen,
-  EAGLE_Oxygen,
-  EAGLE_Neon,
-  EAGLE_Magnesium,
-  EAGLE_Silicon,
-  EAGLE_Iron
-};
-
-void set_quantities(struct part *restrict p, 
-                    const struct unit_system *restrict us,
-		    const struct cooling_function_data *restrict cooling,
-		    const struct cosmology *restrict cosmo,
-		    const struct phys_const *restrict internal_const,
-		    float nh,
-		    double u){
-     
-  const float gamma = 5.0/3.0;
-  double hydrogen_number_density = nh*pow(units_cgs_conversion_factor(us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo->z,3)));
-  p->rho = hydrogen_number_density*internal_const->const_proton_mass*
-          (1.0+p->chemistry_data.metal_mass_fraction[EAGLE_Helium]/p->chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-
-  float scale_factor = 1.0/(1.0+cosmo->z);
-  float pressure = (u*pow(scale_factor,3))/cooling->internal_energy_scale*p->rho*(gamma -1.0);
-  p->entropy = pressure/(pow(p->rho,gamma));
-}
-
-void construct_1d_tables_test(const struct part *restrict p,
-			 const struct cooling_function_data *restrict cooling,
-			 const struct cosmology *restrict cosmo,
-			 const struct phys_const *restrict internal_const,
-			 double *temp_table,
-			 double *H_plus_He_heat_table,
-			 double *H_plus_He_electron_abundance_table,
-			 double *element_electron_abundance_table,
-			 double *element_print_cooling_table,
-			 double *element_cooling_table,
-			 float *abundance_ratio,
-			 float *ub, float *lb){
-
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                             internal_const) *
-                cooling->number_density_scale;
-
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-  printf("testCooling 1d z_index dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e\n", z_index,dz, n_h_i, d_n_h, He_i, d_He);
-
-  construct_1d_table_from_4d(
-      p, cooling, cosmo, internal_const,
-      cooling->table.element_cooling.temperature,
-      z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-      cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
-  construct_1d_table_from_4d(
-      p, cooling, cosmo, internal_const,
-      cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-      cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
-  construct_1d_table_from_4d(
-      p, cooling, cosmo, internal_const,
-      cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
-      dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
-      cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-  construct_1d_table_from_4d_elements(
-      p, cooling, cosmo, internal_const,
-      cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
-      n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table,abundance_ratio, ub, lb);
-  construct_1d_print_table_from_4d_elements(
-      p, cooling, cosmo, internal_const,
-      cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
-      n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table,abundance_ratio, ub, lb);
-  construct_1d_table_from_3d(
-      p, cooling, cosmo, internal_const,
-      cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
-      n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-    
-}
-
-void compare_dlambda_du(
-  const struct unit_system *restrict us,
-  struct part *restrict p,
-  const struct xpart *restrict xp,
-  const struct phys_const *restrict internal_const,
-  const struct cooling_function_data *restrict cooling,
-  const struct cosmology *restrict cosmo){
-
-  double H_plus_He_heat_table[176];              
-  double H_plus_He_electron_abundance_table[176];
-  double temp_table[176];  
-  double element_cooling_table[176];         
-  double element_print_cooling_table[9 * 176];  
-  double element_electron_abundance_table[176];
-  double rate_element_table[11];
-  float *abundance_ratio, cooling_du_dt1, cooling_du_dt2;
-  double dlambda_du1, dlambda_du2;
-
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(p, cooling, abundance_ratio);
-  
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-
-  float nh = 1.0e-1, u = 1.0e9;
-  set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
-  const float log_10_e = 0.43429448190325182765;
-  float upper_bound = cooling->Temp[cooling->N_Temp-1]/log_10_e;
-  float lower_bound = cooling->Temp[0]/log_10_e;
-  construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table, 
-                           H_plus_He_heat_table, H_plus_He_electron_abundance_table, 
-			   element_electron_abundance_table, element_print_cooling_table, 
-			   element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-  
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                           internal_const) *
-              cooling->number_density_scale;
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-  printf("testCooling.c 4d z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e \n", z_index, dz, n_h_i, d_n_h, He_i, d_He);
-  int nt = 10;
-  for(int i = 0; i < nt; i++){
-    u = pow(10.0,9 + i);
-    set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-    cooling_du_dt1 = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du1,H_plus_He_heat_table,
-                      H_plus_He_electron_abundance_table,
-       	              element_print_cooling_table,
-       	              element_electron_abundance_table,
-       	              temp_table,
-       	              z_index, dz, n_h_i, d_n_h, He_i, d_He,
-       	              p,cooling,cosmo,internal_const,rate_element_table,
-       	              abundance_ratio);
-    cooling_du_dt2 = eagle_metal_cooling_rate(log10(u),
-       	              &dlambda_du2,z_index, dz, n_h_i, d_n_h, He_i, d_He,
-       	              p,cooling,cosmo,internal_const,NULL,
-       	              abundance_ratio);
-    printf("u du_dt_1d du_dt_4d dlambda_du_1d dlambda_du_4d %.5e %.5e %.5e %.5e %.5e\n",u,cooling_du_dt1,cooling_du_dt2,dlambda_du1,dlambda_du2);
-  }
-}
-
-void compare_temp(
-  const struct unit_system *restrict us,
-  struct part *restrict p,
-  const struct xpart *restrict xp,
-  const struct phys_const *restrict internal_const,
-  const struct cooling_function_data *restrict cooling,
-  const struct cosmology *restrict cosmo){
-
-  double H_plus_He_heat_table[176];              
-  double H_plus_He_electron_abundance_table[176];
-  double temp_table[176];  
-  double element_cooling_table[176];         
-  double element_print_cooling_table[9 * 176];  
-  double element_electron_abundance_table[176];
-  float *abundance_ratio;
-
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(p, cooling, abundance_ratio);
-  
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-
-  float nh = 1.0e-1, u = 1.0e9;
-  set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-  const float log_10_e = 0.43429448190325182765;
-  float upper_bound = cooling->Temp[cooling->N_Temp-1]/log_10_e;
-  float lower_bound = cooling->Temp[0]/log_10_e;
-
-  construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table, 
-                           H_plus_He_heat_table, H_plus_He_electron_abundance_table, 
-			   element_electron_abundance_table, element_print_cooling_table, 
-			   element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-  float T1d, T4d, delta_u;
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                           internal_const) *
-              cooling->number_density_scale;
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-  printf("testCooling.c z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e \n", z_index, dz, n_h_i, d_n_h, He_i, d_He);
-  int nt = 10;
-  for(int i = 0; i < nt; i++){
-    u = pow(10.0,9 + i);
-    set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-    T1d = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                   p,cooling, cosmo, internal_const);
-    T4d = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-        		 dz, d_n_h, d_He, p,cooling, cosmo, internal_const);
-    printf("u T1d T4d %.5e %.5e %.5e\n",u,pow(10.0,T1d),pow(10.0,T4d));
-  }
-}
-
-int main(int argc, char **argv) {
-  /* Declare relevant structs */
-  struct swift_params *params = malloc(sizeof(struct swift_params));
-  struct unit_system us;
-  struct chemistry_global_data chem_data;
-  struct part p;
-  struct xpart xp;
-  struct phys_const internal_const;
-  struct cooling_function_data cooling;
-  struct cosmology cosmo;
-  char *parametersFileName = "./testCooling.yml";
-
-  int tables = 0;
-
-  /* Read the parameter file */
-  if (params == NULL) error("Error allocating memory for the parameter file.");
-  message("Reading runtime parameters from file '%s'", parametersFileName);
-  parser_read_file(parametersFileName, params);
-
-  /* And dump the parameters as used. */
-  parser_write_params_to_file(params, "used_parameters.yml");
-
-  /* Init units */
-  units_init_from_params(&us, params, "InternalUnitSystem");
-  phys_const_init(&us, params, &internal_const);
-
-  /* Init chemistry */
-  chemistry_init(params, &us, &internal_const, &chem_data);
-  chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
-  chemistry_print(&chem_data);
-
-  /* Init cosmology */
-  cosmology_init(params, &us, &internal_const, &cosmo);
-  cosmology_print(&cosmo);
-  cosmo.z = 0.9090909;
-
-  /* Init cooling */
-  cooling_init(params, &us, &internal_const, &cooling);
-  cooling_print(&cooling);
-
-  /* Calculate abundance ratios */
-  float *abundance_ratio;
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-  
-  /* Declare 1D tables */
-  double H_plus_He_heat_table[176];              
-  double H_plus_He_electron_abundance_table[176];
-  double temp_table[176];  
-  double element_cooling_table[176];         
-  double element_print_cooling_table[9 * 176];  
-  double element_electron_abundance_table[176];
-  
-  float nh;
-
-  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
-  get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
-  const float log_10_e = 0.43429448190325182765;
-  float upper_bound = cooling.Temp[cooling.N_Temp-1]/log_10_e;
-  float lower_bound = cooling.Temp[0]/log_10_e;
-
-  /* Loop over densities */
-  int nt = 250, n_nh = 6;
-  double u = pow(10.0,10);
-  if (argc == 1 || strcmp(argv[1], "nh") == 0){
-    for(int i = 0; i < n_nh; i++){
-      /* Open files */
-      char output_filename[21];
-      sprintf(output_filename, "%s%d%s", "cooling_output_", i, ".dat");
-      FILE *output_file = fopen(output_filename, "w");
-      if (output_file == NULL) {
-        printf("Error opening file!\n");
-        exit(1);
-      }
-
-      nh = pow(10.0,-i);
-      set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
-      construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const, 
-        			temp_table, H_plus_He_heat_table, 
-          		H_plus_He_electron_abundance_table, 
-          		element_electron_abundance_table, 
-          		element_print_cooling_table, 
-          		element_cooling_table, 
-          		abundance_ratio, &upper_bound, &lower_bound);
-      get_index_1d(cooling.nH, cooling.N_nH, log10(nh), &n_h_i, &d_n_h);
-
-      for(int j = 0; j < nt; j++){
-        set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,11.0 + j*8.0/nt));
-	u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
-        double dlambda_du;
-        float delta_u, cooling_du_dt, logT;
-        if (tables == 1) {
-	  cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du,H_plus_He_heat_table,
-                                  H_plus_He_electron_abundance_table,
-             	              element_print_cooling_table,
-             	              element_electron_abundance_table,
-             	              temp_table,
-             	              z_index, dz, n_h_i, d_n_h, He_i, d_He,
-             	              &p,&cooling,&cosmo,&internal_const,NULL,
-             	              abundance_ratio);
-          logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                         &p,&cooling, &cosmo, &internal_const);
-	} else {
-          cooling_du_dt = eagle_metal_cooling_rate(log10(u),
-             	              &dlambda_du,z_index, dz, n_h_i, d_n_h, He_i, d_He,
-             	              &p,&cooling,&cosmo,&internal_const,NULL,
-             	              abundance_ratio);
-          logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-	      		 dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
-	}
-        float temperature_swift = pow(10.0,logT);
-
-        fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
-      }
-      fclose(output_file);
-    }
-  }
-  if (argc == 1 || strcmp(argv[1],"metals") == 0) {
-    printf("here \n");
-    char output_filename[21];
-    sprintf(output_filename, "%s", "cooling_output.dat");
-    FILE *output_file = fopen(output_filename, "w");
-    if (output_file == NULL) {
-      printf("Error opening file!\n");
-      exit(1);
-    }
-
-    // set hydrogen number density, construct 1d tables
-    //nh = pow(10.0,strtod(argv[2],NULL)); 
-    nh = 1.0e-1;
-    u = pow(10.0,14.0);
-    set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
-    float nh_swift = chemistry_get_number_density(&p, &cosmo, chemistry_element_H, &internal_const)*cooling.number_density_scale;
-    printf("hydrogen number density %.5e\n", nh_swift);
-    construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const, 
-      			temp_table, H_plus_He_heat_table, 
-        		H_plus_He_electron_abundance_table, 
-        		element_electron_abundance_table, 
-        		element_print_cooling_table, 
-        		element_cooling_table, 
-        		abundance_ratio, &upper_bound, &lower_bound);
-    float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
-                                             &internal_const) *
-                cooling.number_density_scale;
-    get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-    printf("testcooling.c z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e\n",z_index, dz, n_h_i, d_n_h, He_i, d_He);
-    for(int j = 0; j < nt; j++){
-      set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
-      u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
-      float delta_u, cooling_du_dt, logT;
-      if (tables == 1) {
-        cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
-                              H_plus_He_electron_abundance_table,
-           	              element_print_cooling_table,
-           	              element_electron_abundance_table,
-           	              temp_table,
-           	              z_index, dz, n_h_i, d_n_h, He_i, d_He,
-           	              &p,&cooling,&cosmo,&internal_const,
-           	              abundance_ratio);
-        logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                       &p,&cooling, &cosmo, &internal_const);
-      } else {
-        cooling_du_dt = eagle_print_metal_cooling_rate(
-			      z_index, dz, n_h_i, d_n_h, He_i, d_He,
-           	              &p,&cooling,&cosmo,&internal_const,
-           	              abundance_ratio);
-      }
-      fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
-    }
-    fclose(output_file);
-  }
-
-  if (argc == 1 || strcmp(argv[1],"compare_temp") == 0) compare_temp(&us, &p, &xp, &internal_const, &cooling, &cosmo);
-  if (argc == 1 || strcmp(argv[1],"compare_dlambda_du") == 0) compare_dlambda_du(&us, &p, &xp, &internal_const, &cooling, &cosmo);
-
-  free(params);
-  return 0;
-}
-
-
diff --git a/examples/CoolingTest_Alexei/testCooling.yml b/examples/CoolingTest_Alexei/testCooling.yml
deleted file mode 100644
index f795ed6836786ae1ad72234fda8184c946e4ebe9..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling.yml
+++ /dev/null
@@ -1,94 +0,0 @@
-# Define the system of units to use internally. 
-InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
-  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
-  UnitCurrent_in_cgs:  1             # Amperes
-  UnitTemp_in_cgs:     1             # Kelvin
-
-# Cosmological parameters
-Cosmology:
-  h:              0.6777        # Reduced Hubble constant
-  #a_begin:        0.9090909     # Initial scale-factor of the simulation
-  a_begin:        1.0           # Initial scale-factor of the simulation
-  a_end:          1.0           # Final scale factor of the simulation
-  Omega_m:        0.307         # Matter density parameter
-  Omega_lambda:   0.693         # Dark-energy density parameter
-  Omega_b:        0.0455        # Baryon density parameter
-
-# Parameters governing the time integration
-TimeIntegration:
-  time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1e-2  # The end time of the simulation (in internal units).
-  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
-  dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
-  
-Scheduler:
-  max_top_level_cells:    15
-  
-# Parameters governing the snapshots
-Snapshots:
-  basename:            eagle # Common part of the name of output files
-  scale_factor_first:  0.92  # Scale-factor of the first snaphot (cosmological run)
-  time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
-  delta_time:          1.10  # Time difference between consecutive outputs (in internal units)
-  compression:         1
-
-# Parameters governing the conserved quantities statistics
-Statistics:
-  scale_factor_first:  0.92 # Scale-factor of the first stat dump (cosmological run)
-  time_first:          0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
-  delta_time:          1.05 # Time between statistics output
-
-# Parameters for the self-gravity scheme
-Gravity:
-  eta:                    0.025     # Constant dimensionless multiplier for time integration.
-  theta:                  0.85      # Opening angle (Multipole acceptance criterion)
-  comoving_softening:     0.0026994 # Comoving softening length (in internal units).
-  max_physical_softening: 0.0007    # Physical softening length (in internal units).
-  
-# Parameters for the hydrodynamics scheme
-SPH:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
-  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
-  minimal_temperature:   100      # (internal units)
-
-# Parameters related to the initial conditions
-InitialConditions:
-  file_name:  ./EAGLE_ICs_12.hdf5     # The file to read
-  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
-
-EAGLEChemistry:
-  InitMetallicity:         0.014
-  InitAbundance_Hydrogen:  0.70649785
-  InitAbundance_Helium:    0.28055534
-  InitAbundance_Carbon:    2.0665436e-3
-  InitAbundance_Nitrogen:  8.3562563e-4
-  InitAbundance_Oxygen:    5.4926244e-3
-  InitAbundance_Neon:      1.4144605e-3
-  InitAbundance_Magnesium: 5.907064e-4
-  InitAbundance_Silicon:   6.825874e-4
-  InitAbundance_Iron:      1.1032152e-3
-  CalciumOverSilicon:      0.0941736
-  SulphurOverSilicon:      0.6054160
-  #InitMetallicity:         0.
-  #InitAbundance_Hydrogen:  0.752
-  #InitAbundance_Helium:    0.248
-  #InitAbundance_Carbon:    0.000
-  #InitAbundance_Nitrogen:  0.000
-  #InitAbundance_Oxygen:    0.000
-  #InitAbundance_Neon:      0.000
-  #InitAbundance_Magnesium: 0.000
-  #InitAbundance_Silicon:   0.000
-  #InitAbundance_Iron:      0.000
-  #CalciumOverSilicon:      0.0941736
-  #SulphurOverSilicon:      0.6054160
-
-EagleCooling:
-  filename:                /cosma5/data/Eagle/BG_Tables/CoolingTables/
-  reionisation_redshift:   20
-  he_reion:                0
-  he_reion_z_center:       3.5
-  he_reion_z_sigma:        0.5
-  he_reion_ev_pH:          2.0
-
diff --git a/examples/CoolingTest_Alexei/testCooling_messy.c b/examples/CoolingTest_Alexei/testCooling_messy.c
deleted file mode 100644
index b5c00e5625c43e57e46c5f8426640b9f511bc58a..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling_messy.c
+++ /dev/null
@@ -1,354 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-int main(int argc, char **argv) {
-
-  /* Read the parameter file */
-  struct swift_params *params = malloc(sizeof(struct swift_params));
-  struct unit_system us;
-  struct chemistry_global_data chem_data;
-  struct part p;
-  struct xpart xp;
-  struct phys_const internal_const;
-  struct cooling_function_data cooling;
-  struct cosmology cosmo;
-  char *parametersFileName = "./testCooling.yml";
-  enum table_index {
-    EAGLE_Hydrogen = 0,
-    EAGLE_Helium,
-    EAGLE_Carbon,
-    EAGLE_Nitrogen,
-    EAGLE_Oxygen,
-    EAGLE_Neon,
-    EAGLE_Magnesium,
-    EAGLE_Silicon,
-    EAGLE_Iron
-  };
-
-  if (params == NULL) error("Error allocating memory for the parameter file.");
-  message("Reading runtime parameters from file '%s'", parametersFileName);
-  parser_read_file(parametersFileName, params);
-
-  /* And dump the parameters as used. */
-  parser_write_params_to_file(params, "used_parameters.yml");
-
-  double UnitMass_in_cgs = 1.989e43;
-  double UnitLength_in_cgs = 3.085678e24;
-  double UnitVelocity_in_cgs = 1e5;
-  double UnitCurrent_in_cgs = 1;
-  double UnitTemp_in_cgs = 1;
-  units_init(&us, UnitMass_in_cgs, UnitLength_in_cgs, UnitVelocity_in_cgs, UnitCurrent_in_cgs, UnitTemp_in_cgs);
-  phys_const_init(&us, params, &internal_const);
-
-  double hydrogen_number_density_cgs = 1e-4;
-  double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
-      temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
-  u_cgs = 0;
-  cooling_du_dt = 0;
-  temperature_cgs = 0;
-  newton_func = 0;
-  // int n_t_i = 2000;
-  dt_cgs = 4.0e-8 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
-
-  char *output_filename = "cooling_output.dat";
-  FILE *output_file = fopen(output_filename, "w");
-  if (output_file == NULL) {
-    printf("Error opening file!\n");
-    exit(1);
-  }
-  char *output_filename2 = "temperature_output.dat";
-  FILE *output_file2 = fopen(output_filename2, "w");
-  if (output_file2 == NULL) {
-    printf("Error opening file!\n");
-    exit(1);
-  }
-  char *output_filename3 = "newton_output.dat";
-  FILE *output_file3 = fopen(output_filename3, "w");
-  if (output_file3 == NULL) {
-    printf("Error opening file!\n");
-    exit(1);
-  }
-
-  gamma = 5.0 / 3.0;
-
-  chemistry_init(params, &us, &internal_const, &chem_data);
-  chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
-  chemistry_print(&chem_data);
-
-  cosmology_init(params, &us, &internal_const, &cosmo);
-  cosmology_print(&cosmo);
-
-  u = 1.0 * pow(10.0, 11) /
-      (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
-       units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
-  pressure = u * p.rho * (gamma - 1.0);
-  hydrogen_number_density =
-      hydrogen_number_density_cgs *
-      pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
-  p.rho = hydrogen_number_density * internal_const.const_proton_mass *
-          (1.0 +
-           p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
-               p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-  p.entropy = pressure / (pow(p.rho, gamma));
-
-  cooling_init(params, &us, &internal_const, &cooling);
-  cooling_print(&cooling);
-
-  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
-         (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
-                                             &internal_const) *
-                cooling.number_density_scale;
-  ratefact = inn_h * (XH / eagle_proton_mass_cgs);
-  
-  float *abundance_ratio;
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
-  // construct 1d table of cooling rates wrt temperature
-  double H_plus_He_heat_table[176];              
-  double H_plus_He_electron_abundance_table[176];
-  double temp_table[176];  
-  double element_cooling_table[176];         
-  double element_print_cooling_table[9 * 176];  
-  double element_electron_abundance_table[176];
-
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
-  get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-  construct_1d_table_from_4d(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.temperature,
-      z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
-      cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, temp_table);
-  construct_1d_table_from_4d(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.H_plus_He_heating, z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
-      cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, H_plus_He_heat_table);
-  construct_1d_table_from_4d(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.H_plus_He_electron_abundance, z_index,
-      dz, cooling.N_Redshifts, n_h_i, d_n_h, cooling.N_nH, He_i, d_He,
-      cooling.N_He, cooling.N_Temp, H_plus_He_electron_abundance_table);
-  construct_1d_table_from_4d_elements(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
-      n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_cooling_table,abundance_ratio);
-  construct_1d_print_table_from_4d_elements(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
-      n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_print_cooling_table,abundance_ratio);
-  construct_1d_table_from_3d(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.electron_abundance, z_index, dz, cooling.N_Redshifts,
-      n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_electron_abundance_table);
-
-
-  //
-  // printf("cooling table values \n");
-  // for( int j=0; j < 176; j++) {
-  //    printf("   %.5e",  H_plus_He_heat_table[j]+element_cooling_table[j]
-  u_ini_cgs = pow(10.0, 14);
-  float delta_u;
-
-
-   //Compute contributions to cooling rate from different metals
-  // int n_t_i = 100;
-  // for(int t_i = 0; t_i < n_t_i; t_i++){
-  //   u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
-  //   //u = u_cgs/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
-  //   u = u_cgs/cooling.internal_energy_scale;
-  //   hydrogen_number_density =
-  //   hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
-  //   p.rho =
-  //   hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-  //   pressure = u*p.rho*(gamma -1.0);
-  //   p.entropy = pressure/(pow(p.rho,gamma));
-  //   double u_swift = hydro_get_physical_internal_energy(&p, &cosmo) *
-  //                cooling.internal_energy_scale;
-
-  //   cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
-  //                     H_plus_He_electron_abundance_table,
-  //      	       element_print_cooling_table,
-  //      	       element_electron_abundance_table,
-  //      	       temp_table,
-  //      	       z_index, dz, n_h_i, d_n_h, He_i, d_He,
-  //      	       &p,&cooling,&cosmo,&internal_const,
-  //      	       abundance_ratio);
-  //   float logT = eagle_convert_u_to_temp_1d_table(log10(u_swift), &delta_u, temp_table,
-  //                  &p,&cooling, &cosmo, &internal_const);
-  //   float temperature_swift = pow(10.0,logT);
-
-  //   fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
-  //   fprintf(output_file2,"%.5e %.5e\n",u_swift*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
-  //   temperature_cgs);
-
-  //   newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
-  //   fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
-  //   //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e
-  //   // \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
-  //}
-   
-   int n_t_i = 100;
-   for(int nh_i = 1; nh_i < 7; nh_i++){
-     char output_filename4[21];
-     sprintf(output_filename4, "%s%d%s", "cooling_output_", nh_i, ".dat");
-     FILE *output_file4 = fopen(output_filename4, "w");
-     if (output_file4 == NULL) {
-       printf("Error opening file!\n");
-       exit(1);
-     }
-     hydrogen_number_density = pow(10.0,-nh_i)*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo.z,3)));
-     p.rho =
-     hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-
-     XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
-     HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
-      (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
-     inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
-                                                &internal_const) *
-                   cooling.number_density_scale;
-     ratefact = inn_h * (XH / eagle_proton_mass_cgs);
-     printf("test cooling hydrogen num dens inn_h %.5e %.5e\n", hydrogen_number_density*cooling.number_density_scale, inn_h);
-     
-     abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-     abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
-     get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
-     get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-     get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-     construct_1d_table_from_4d(
-         &p, &cooling, &cosmo, &internal_const,
-         cooling.table.element_cooling.temperature,
-         z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
-         cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, temp_table);
-     construct_1d_table_from_4d(
-         &p, &cooling, &cosmo, &internal_const,
-         cooling.table.element_cooling.H_plus_He_heating, z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
-         cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, H_plus_He_heat_table);
-     construct_1d_table_from_4d(
-         &p, &cooling, &cosmo, &internal_const,
-         cooling.table.element_cooling.H_plus_He_electron_abundance, z_index,
-         dz, cooling.N_Redshifts, n_h_i, d_n_h, cooling.N_nH, He_i, d_He,
-         cooling.N_He, cooling.N_Temp, H_plus_He_electron_abundance_table);
-     construct_1d_table_from_4d_elements(
-         &p, &cooling, &cosmo, &internal_const,
-         cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
-         n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_cooling_table,abundance_ratio);
-     construct_1d_print_table_from_4d_elements(
-         &p, &cooling, &cosmo, &internal_const,
-         cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
-         n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_print_cooling_table,abundance_ratio);
-     construct_1d_table_from_3d(
-         &p, &cooling, &cosmo, &internal_const,
-         cooling.table.element_cooling.electron_abundance, z_index, dz, cooling.N_Redshifts,
-         n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_electron_abundance_table);
-     for(int t_i = 0; t_i < n_t_i; t_i++){
-       u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
-       u = u_cgs/cooling.internal_energy_scale;
-       pressure = u*p.rho*(gamma -1.0);
-       p.entropy = pressure/(pow(p.rho,gamma));
-       double u_swift = hydro_get_physical_internal_energy(&p, &cosmo) *
-                    cooling.internal_energy_scale;
-
-       double dlambda_du;
-       cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u_swift),&dlambda_du,H_plus_He_heat_table,
-                         H_plus_He_electron_abundance_table,
-          	       element_print_cooling_table,
-          	       element_electron_abundance_table,
-          	       temp_table,
-          	       z_index, dz, n_h_i, d_n_h, He_i, d_He,
-          	       &p,&cooling,&cosmo,&internal_const,NULL,
-          	       abundance_ratio);
-       float logT = eagle_convert_u_to_temp_1d_table(log10(u_swift), &delta_u, temp_table,
-                      &p,&cooling, &cosmo, &internal_const);
-       float temperature_swift = pow(10.0,logT);
-
-       fprintf(output_file4,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
-       fprintf(output_file2,"%.5e %.5e\n",u_swift*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
-       temperature_cgs);
-
-       newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
-       fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
-     }
-     fclose(output_file4);
-   }
-   fclose(output_file);
-   fclose(output_file2);
-   fclose(output_file3);
-
-  //double dLambdaNet_du, LambdaNet;  //, LambdaNext;
-  //float x_init, u_eq = 2.0e12;
-  //for (int j = 5; j < 6; j++) {
-  //  float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0, 0.5 * (j + 5)),
-  //                                                 temp_table, &p, &cooling,
-  //                                                 &cosmo, &internal_const),
-  //        x, du;
-  //  float dt = 2.0e-4 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
-  //  LambdaNet = eagle_cooling_rate_1d_table(
-  //      u_ini, &dLambdaNet_du, H_plus_He_heat_table,
-  //      H_plus_He_electron_abundance_table, element_cooling_table,
-  //      element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
-  //      &internal_const);
-  //  float u_temp = u_ini + LambdaNet * ratefact * dt;
-  //  /* RGB removed this **
-  //  if (u_temp > 0) LambdaNext = eagle_cooling_rate_1d_table(u_temp,
-  //  &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
-  //  element_cooling_table, element_electron_abundance_table, temp_table, &p,
-  //  &cooling, &cosmo, &internal_const);
-  //  if (fabs(LambdaNet - LambdaNext)/LambdaNet < 0.5) {
-  //    u_temp = u_ini;
-  //  } else {
-  //    u_temp =
-  //  eagle_convert_temp_to_u_1d_table(1.0e4,temp_table,&p,&cooling,&cosmo,&internal_const);
-  //              } */
-  //  if (u_temp > u_eq) {
-  //    x_init = log(u_temp);
-  //  } else {
-  //    x_init = log(u_eq);
-  //  }
-  //  x = newton_iter(x_init, u_ini, H_plus_He_heat_table,
-  //                  H_plus_He_electron_abundance_table, element_cooling_table,
-  //                  element_electron_abundance_table, temp_table, &p, &cosmo,
-  //                  &cooling, &internal_const, dt);
-  //  printf(
-  //      "testing newton integration, u_ini, u %.5e %.5e, temperature initial, "
-  //      "final %.5e %.5e\n",
-  //      u_ini, exp(x),
-  //      eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
-  //                                       &cosmo, &internal_const),
-  //      eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
-  //                                       &cosmo, &internal_const));
-  //}
-
-  free(params);
-
-  return 0;
-}
diff --git a/examples/CoolingTest_Alexei/testCooling_script.c b/examples/CoolingTest_Alexei/testCooling_script.c
deleted file mode 100644
index ab4200b67be67e1d45913e56652f9e64913fe3c1..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling_script.c
+++ /dev/null
@@ -1,637 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include <unistd.h>
-#include "cooling.h"
-#include "cooling_struct.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-  
-enum table_index {
-  EAGLE_Hydrogen = 0,
-  EAGLE_Helium,
-  EAGLE_Carbon,
-  EAGLE_Nitrogen,
-  EAGLE_Oxygen,
-  EAGLE_Neon,
-  EAGLE_Magnesium,
-  EAGLE_Silicon,
-  EAGLE_Iron
-};
-
-/*
- * @brief Assign particle density and entropy corresponding to the 
- * hydrogen number density and internal energy specified.
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param internal_const Physical constants data structure
- * @param nh Hydrogen number density (cgs units)
- * @param u Internal energy (cgs units)
- */
-void set_quantities(struct part *restrict p, 
-                    const struct unit_system *restrict us,
-		    const struct cooling_function_data *restrict cooling,
-		    const struct cosmology *restrict cosmo,
-		    const struct phys_const *restrict internal_const,
-		    float nh,
-		    double u){
-     
-  const float gamma = 5.0/3.0;
-  double hydrogen_number_density = nh*pow(units_cgs_conversion_factor(us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo->z,3)));
-  p->rho = hydrogen_number_density*internal_const->const_proton_mass*
-          (1.0+p->chemistry_data.metal_mass_fraction[EAGLE_Helium]/p->chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-
-  float scale_factor = 1.0/(1.0+cosmo->z);
-  float pressure = (u*pow(scale_factor,3))/cooling->internal_energy_scale*p->rho*(gamma -1.0);
-  p->entropy = pressure/(pow(p->rho,gamma));
-}
-
-/*
- * @brief Construct 1d tables from 4d EAGLE tables by 
- * interpolating over redshift, hydrogen number density
- * and helium fraction. 
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param internal_const Physical constants data structure
- * @param temp_table Pointer to 1d interpolated table of temperature values
- * @param H_plus_He_heat_table Pointer to 1d interpolated table of cooling rates
- * due to hydrogen and helium
- * @param H_plus_He_electron_abundance_table Pointer to 1d interpolated table
- * of electron abundances due to hydrogen and helium
- * @param element_electron_abundance_table Pointer to 1d interpolated table 
- * of electron abundances due to metals
- * @param element_print_cooling_table Pointer to 1d interpolated table of
- * cooling rates due to each of the metals
- * @param element_cooling_table Pointer to 1d interpolated table of cooling
- * rates due to the contribution of all the metals
- * @param abundance_ratio Pointer to array of ratios of metal abundances to solar
- * @param ub Upper bound in temperature on table construction 
- * @param lb Lower bound in temperature on table construction
- */
-void construct_1d_tables_test(const struct part *restrict p,
-			 const struct cooling_function_data *restrict cooling,
-			 const struct cosmology *restrict cosmo,
-			 const struct phys_const *restrict internal_const,
-			 double *temp_table,
-			 double *H_plus_He_heat_table,
-			 double *H_plus_He_electron_abundance_table,
-			 double *element_electron_abundance_table,
-			 double *element_print_cooling_table,
-			 double *element_cooling_table,
-			 float *abundance_ratio,
-			 float *ub, float *lb){
-
-  // obtain mass fractions and number density for particle
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                             internal_const) *
-                cooling->number_density_scale;
-
-  // find redshift, helium fraction, hydrogen number density indices and offsets
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-  // construct tables
-  //construct_1d_table_from_4d(
-  //    p, cooling, cosmo, internal_const,
-  //    cooling->table.element_cooling.temperature,
-  //    z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-  //    cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
-  //construct_1d_table_from_4d(
-  //    p, cooling, cosmo, internal_const,
-  //    cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-  //    cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
-  //construct_1d_table_from_4d(
-  //    p, cooling, cosmo, internal_const,
-  //    cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
-  //    dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
-  //    cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-  //construct_1d_table_from_4d_elements(
-  //    p, cooling, cosmo, internal_const,
-  //    cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
-  //    n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table,abundance_ratio, ub, lb);
-  //construct_1d_print_table_from_4d_elements(
-  //    p, cooling, cosmo, internal_const,
-  //    cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
-  //    n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table,abundance_ratio, ub, lb);
-  //construct_1d_table_from_3d(
-  //    p, cooling, cosmo, internal_const,
-  //    cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
-  //    n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-
-  if (cosmo->z > cooling->reionisation_redshift) { 
-    // Photodissociation table 
-    printf("Eagle testCooling.c photodissociation table redshift reionisation redshift %.5e %.5e\n", cosmo->z, cooling->reionisation_redshift);
-    construct_1d_table_from_3d(p, cooling, cosmo, internal_const, 
-                cooling->table.photodissociation_cooling.temperature, 
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb); 
-    construct_1d_table_from_3d( 
-                p, cooling, cosmo, internal_const, 
-                cooling->table.photodissociation_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,   
-                He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb); 
-    construct_1d_table_from_3d( 
-                p, cooling, cosmo, internal_const, 
-                cooling->table.photodissociation_cooling.H_plus_He_electron_abundance, 
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, 
-                cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-    construct_1d_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.photodissociation_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_print_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.photodissociation_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_table_from_2d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.photodissociation_cooling.electron_abundance,
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-  } else if (cosmo->z > cooling->Redshifts[cooling->N_Redshifts - 1]) {
-    printf("Eagle testCooling.c no compton table redshift max redshift %.5e %.5e\n", cosmo->z, cooling->Redshifts[cooling->N_Redshifts - 1]);
-    // High redshift table
-    construct_1d_table_from_3d(p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.temperature,
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
-    construct_1d_table_from_3d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,  
-                He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
-    construct_1d_table_from_3d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.H_plus_He_electron_abundance,
-                n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He,
-                cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-    construct_1d_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_print_table_from_3d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.metal_heating, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_table_from_2d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.no_compton_cooling.electron_abundance,
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-  } else {
-    // Normal tables 
-    printf("Eagle testCooling.c normal table redshift %.5e\n", cosmo->z);
-    construct_1d_table_from_4d(p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.temperature,
-                z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-                cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
-    construct_1d_table_from_4d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
-                cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
-    construct_1d_table_from_4d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
-                dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
-                cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
-    construct_1d_table_from_4d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_print_table_from_4d_elements(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts, 
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
-    construct_1d_table_from_3d(
-                p, cooling, cosmo, internal_const,
-                cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
-                n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-  } 
-    
-}
-
-/*
- * @brief Compare calculation of dlambda/du (gradient of cooling rate, 
- * needed for Newton's method) between interpolating 1d tables and 
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_dlambda_du(
-  const struct unit_system *restrict us,
-  struct part *restrict p,
-  const struct xpart *restrict xp,
-  const struct phys_const *restrict internal_const,
-  const struct cooling_function_data *restrict cooling,
-  const struct cosmology *restrict cosmo){
-
-  // allocate tables
-  double H_plus_He_heat_table[cooling->N_Temp];              
-  double H_plus_He_electron_abundance_table[cooling->N_Temp];
-  double temp_table[cooling->N_Temp];  
-  double element_cooling_table[cooling->N_Temp];         
-  double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];  
-  double element_electron_abundance_table[cooling->N_Temp];
-  double rate_element_table[cooling->N_Elements+2];
-  float *abundance_ratio, cooling_du_dt1, cooling_du_dt2;
-  double dlambda_du1, dlambda_du2;
-
-  // calculate ratio of particle metal abundances to solar 
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(p, cooling, abundance_ratio);
-  
-  // set hydrogen number density and internal energy
-  float nh = 1.0e-1, u = 1.0e9;
-  set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-  
-  // extract hydrogen and helium mass fractions
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                           internal_const) *
-              cooling->number_density_scale;
-  
-  // find redshift, hydrogen number density and helium fraction indices and offsets
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-
-  // construct tables
-  float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
-  float lower_bound = cooling->Temp[0]/eagle_log_10_e;
-  construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table, 
-                           H_plus_He_heat_table, H_plus_He_electron_abundance_table, 
-			   element_electron_abundance_table, element_print_cooling_table, 
-			   element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-  
-  // calculate dlambda/du for different values of internal energy
-  int nt = 10;
-  for(int i = 0; i < nt; i++){
-    u = pow(10.0,9 + i);
-    set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-    cooling_du_dt1 = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du1,H_plus_He_heat_table,
-                      H_plus_He_electron_abundance_table,
-       	              element_print_cooling_table,
-       	              element_electron_abundance_table,
-       	              temp_table,
-       	              z_index, dz, n_h_i, d_n_h, He_i, d_He,
-       	              p,cooling,cosmo,internal_const,rate_element_table,
-       	              abundance_ratio);
-    cooling_du_dt2 = eagle_metal_cooling_rate(log10(u),
-       	              &dlambda_du2,z_index, dz, n_h_i, d_n_h, He_i, d_He,
-       	              p,cooling,cosmo,internal_const,NULL,
-       	              abundance_ratio);
-    printf("u du_dt_1d du_dt_4d dlambda_du_1d dlambda_du_4d %.5e %.5e %.5e %.5e %.5e\n",u,cooling_du_dt1,cooling_du_dt2,dlambda_du1,dlambda_du2);
-  }
-}
-
-/*
- * @brief Compare calculating temperature between interpolating 1d tables and 
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_temp(
-  const struct unit_system *restrict us,
-  struct part *restrict p,
-  const struct xpart *restrict xp,
-  const struct phys_const *restrict internal_const,
-  const struct cooling_function_data *restrict cooling,
-  const struct cosmology *restrict cosmo){
-
-  // allocate tables
-  double H_plus_He_heat_table[cooling->N_Temp];              
-  double H_plus_He_electron_abundance_table[cooling->N_Temp];
-  double temp_table[cooling->N_Temp];  
-  double element_cooling_table[cooling->N_Temp];         
-  double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];  
-  double element_electron_abundance_table[cooling->N_Temp];
-  float *abundance_ratio;
-
-  // calculate ratio of particle metal abundances to solar 
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(p, cooling, abundance_ratio);
-  
-  // set hydrogen number density and internal energy
-  float nh = 1.0e-1, u = 1.0e9;
-  set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
-  // extract hydrogen and helium mass fractions
-  float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                           internal_const) *
-              cooling->number_density_scale;
-  
-  // find redshift, hydrogen number density and helium fraction indices and offsets
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-  get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-  get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-  // construct tables
-  float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
-  float lower_bound = cooling->Temp[0]/eagle_log_10_e;
-  construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table, 
-                           H_plus_He_heat_table, H_plus_He_electron_abundance_table, 
-			   element_electron_abundance_table, element_print_cooling_table, 
-			   element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-  
-  // calculate temperature for different values of internal energy
-  float T1d, T4d, delta_u;
-  int nt = 10;
-  for(int i = 0; i < nt; i++){
-    u = pow(10.0,9 + i);
-    set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-    T1d = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                   p,cooling, cosmo, internal_const);
-    T4d = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-        		 dz, d_n_h, d_He, p,cooling, cosmo, internal_const);
-    printf("u T1d T4d %.5e %.5e %.5e\n",u,pow(10.0,T1d),pow(10.0,T4d));
-  }
-}
-
-/*
- * @brief Produces contributions to cooling rates for different 
- * hydrogen number densities, from different metals, 
- * tests 1d and 4d table interpolations produce 
- * same results for cooling rate, dlambda/du and temperature.
- *
- * @param pass string "nh" to produce cooling rates for one internal
- * energy, different redshifts
- * @param pass string "metals" to produce contribution to cooling
- * rates from different metals
- * @param pass string "compare_temp" to compare temperature interpolation
- * from 1d and 4d tables
- * @param pass string "compare_dlambda_du" to compare dlambda/du calculation
- * from 1d and 4d tables
- * @param if no string passed, all of the above are performed.
- */
-int main(int argc, char **argv) {
-  // Declare relevant structs
-  struct swift_params *params = malloc(sizeof(struct swift_params));
-  struct unit_system us;
-  struct chemistry_global_data chem_data;
-  struct part p;
-  struct xpart xp;
-  struct phys_const internal_const;
-  struct cooling_function_data cooling;
-  struct cosmology cosmo;
-  char *parametersFileName = "./testCooling.yml";
-
-  float nh;
-
-  // Read options
-  int param, tables = 0;
-  float redshift = -1.0, log_10_nh = 100;
-  while ((param = getopt(argc, argv, "z:d:t")) != -1)
-  switch(param){
-    case 'z':
-      redshift = atof(optarg);
-      break;
-    case 'd':
-      log_10_nh = atof(optarg);
-      break;
-    case 't':
-      tables = 1;
-      break;
-    case '?':
-      if (optopt == 'z')
-        printf ("Option -%c requires an argument.\n", optopt);
-      else
-        printf ("Unknown option character `\\x%x'.\n",
-                 optopt);
-      error("invalid option(s) to testCooling");
-  }
-
-  // Read the parameter file
-  if (params == NULL) error("Error allocating memory for the parameter file.");
-  message("Reading runtime parameters from file '%s'", parametersFileName);
-  parser_read_file(parametersFileName, params);
-
-  // And dump the parameters as used. 
-  parser_write_params_to_file(params, "used_parameters.yml");
-
-  // Init units 
-  units_init_from_params(&us, params, "InternalUnitSystem");
-  phys_const_init(&us, params, &internal_const);
-
-  // Init chemistry 
-  chemistry_init(params, &us, &internal_const, &chem_data);
-  chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
-  chemistry_print(&chem_data);
-
-  // Init cosmology 
-  cosmology_init(params, &us, &internal_const, &cosmo);
-  cosmology_print(&cosmo);
-  if (redshift == -1.0) {
-    cosmo.z = 3.0;
-  } else {
-    cosmo.z = redshift;
-  }
-
-  // Init cooling 
-  cooling_init(params, &us, &internal_const, &cooling);
-  cooling_print(&cooling);
-
-  // Calculate abundance ratios 
-  float *abundance_ratio;
-  abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
-  abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-  
-  // Declare 1D tables 
-  double H_plus_He_heat_table[cooling.N_Temp];              
-  double H_plus_He_electron_abundance_table[cooling.N_Temp];
-  double temp_table[cooling.N_Temp];  
-  double element_cooling_table[cooling.N_Temp];         
-  double element_print_cooling_table[cooling.N_Elements * cooling.N_Temp];  
-  double element_electron_abundance_table[cooling.N_Temp];
-  
-  // extract mass fractions, calculate table indices and offsets 
-  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
-                 (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
-  int z_index,He_i,n_h_i;
-  float dz,d_He,d_n_h;
-  get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
-  get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
-  float upper_bound = cooling.Temp[cooling.N_Temp-1]/eagle_log_10_e;
-  float lower_bound = cooling.Temp[0]/eagle_log_10_e;
-
-  int nt = 250, n_nh = 6;
-  double u = pow(10.0,10);
-  //if (argc == 1 || strcmp(argv[1], "nh") == 0){
-  // Calculate cooling rates at different densities 
-    for(int i = 0; i < n_nh; i++){
-      // Open files 
-      char output_filename[21];
-      sprintf(output_filename, "%s%d%s", "cooling_output_", i, ".dat");
-      FILE *output_file = fopen(output_filename, "w");
-      if (output_file == NULL) {
-        printf("Error opening file!\n");
-        exit(1);
-      }
-
-      // set hydrogen number density, construct tables
-      nh = pow(10.0,-i);
-      set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
-      construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const, 
-        			temp_table, H_plus_He_heat_table, 
-          		H_plus_He_electron_abundance_table, 
-          		element_electron_abundance_table, 
-          		element_print_cooling_table, 
-          		element_cooling_table, 
-          		abundance_ratio, &upper_bound, &lower_bound);
-      get_index_1d(cooling.nH, cooling.N_nH, log10(nh), &n_h_i, &d_n_h);
-
-      for(int j = 0; j < nt; j++){
-	// set internal energy
-        set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,11.0 + j*8.0/nt));
-	u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
-        double dlambda_du;
-        float delta_u, cooling_du_dt, logT;
-
-	// calculate cooling rates using 1d tables
-        if (tables == 1) {
-	  cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du,H_plus_He_heat_table,
-                                  H_plus_He_electron_abundance_table,
-             	              element_print_cooling_table,
-             	              element_electron_abundance_table,
-             	              temp_table,
-             	              z_index, dz, n_h_i, d_n_h, He_i, d_He,
-             	              &p,&cooling,&cosmo,&internal_const,NULL,
-             	              abundance_ratio);
-          logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                         &p,&cooling, &cosmo, &internal_const);
-	} else {
-	// calculate cooling rates using 4d tables
-          cooling_du_dt = eagle_metal_cooling_rate(log10(u),
-             	              &dlambda_du,z_index, dz, n_h_i, d_n_h, He_i, d_He,
-             	              &p,&cooling,&cosmo,&internal_const,NULL,
-             	              abundance_ratio);
-          logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-	      		 dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
-	}
-        float temperature_swift = pow(10.0,logT);
-
-        fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
-      }
-      fclose(output_file);
-    }
-  //}
-  // Calculate contributions from metals to cooling rate
-  //if (argc >= 1 || strcmp(argv[1],"metals") == 0) {
-    // open file
-    char output_filename[21];
-    sprintf(output_filename, "%s", "cooling_output.dat");
-    FILE *output_file = fopen(output_filename, "w");
-    if (output_file == NULL) {
-      printf("Error opening file!\n");
-      exit(1);
-    }
-
-    // set hydrogen number density, construct 1d tables
-    if (log_10_nh == 100) {
-      nh = 1.0e-1;
-    } else {
-      nh = pow(10.0,log_10_nh);
-    }
-    //if (argc > 2) nh = pow(10.0,strtod(argv[2],NULL));
-    printf("Eagle testcooling.c nh %.5e\n", nh);
-    u = pow(10.0,14.0);
-    set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
-    construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const, 
-      			temp_table, H_plus_He_heat_table, 
-        		H_plus_He_electron_abundance_table, 
-        		element_electron_abundance_table, 
-        		element_print_cooling_table, 
-        		element_cooling_table, 
-        		abundance_ratio, &upper_bound, &lower_bound);
-    float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
-                                             &internal_const) *
-                cooling.number_density_scale;
-    printf("Eagle testcooling.c inn_h %.5e\n", inn_h);
-    get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-    // Loop over internal energy
-    for(int j = 0; j < nt; j++){
-      set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
-      u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
-      float delta_u, cooling_du_dt, logT;
-
-      // calculate cooling rates using 1d tables
-      if (tables == 1) {
-        cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
-                              H_plus_He_electron_abundance_table,
-           	              element_print_cooling_table,
-           	              element_electron_abundance_table,
-           	              temp_table,
-           	              z_index, dz, n_h_i, d_n_h, He_i, d_He,
-           	              &p,&cooling,&cosmo,&internal_const,
-           	              abundance_ratio);
-        logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
-                       &p,&cooling, &cosmo, &internal_const);
-      } else {
-      // calculate cooling rates using 4d tables
-        cooling_du_dt = eagle_print_metal_cooling_rate(
-			      z_index, dz, n_h_i, d_n_h, He_i, d_He,
-           	              &p,&cooling,&cosmo,&internal_const,
-           	              abundance_ratio);
-          logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
-	      		 dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
-      }
-      //float temperature_swift = pow(10.0,logT);
-      //fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
-      fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
-    }
-    fclose(output_file);
-  //}
-
-  // compare temperatures and dlambda/du calculated from 1d and 4d tables
-  compare_temp(&us, &p, &xp, &internal_const, &cooling, &cosmo);
-  compare_dlambda_du(&us, &p, &xp, &internal_const, &cooling, &cosmo);
-
-  free(params);
-  return 0;
-}
-
-
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 6d2b6ce5f647b91c8e778df918341b0539f20046..f4fc179248b45dd79d985518fbfcae4ebeab9d73 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -44,11 +44,12 @@
 #include "physical_constants.h"
 #include "units.h"
 
-const float explicit_tolerance = 0.05;
-const float newton_tolerance =
+static const float explicit_tolerance = 0.05;
+static const float newton_tolerance =
     1.0e-2;  // lower than bisection_tol because using log(u)
-const float bisection_tolerance = 1.0e-6;
-const double bracket_factor = 1.0488088481701;  // sqrt(1.1) to match EAGLE
+static const float bisection_tolerance = 1.0e-6;
+static const double bracket_factor = 1.0488088481701;  // sqrt(1.1) to match EAGLE
+static float print_cooling_rate_contribution_flag = 0;
 
 /*
  * @brief calculates heating due to helium reionization
@@ -73,9 +74,9 @@ eagle_helium_reionization_extraheat(
 #endif
 
   /* Helium reionization */
-  if (cooling->he_reion == 1) {
+  if (cooling->he_reion_flag == 1) {
     double he_reion_erg_pG =
-        cooling->he_reion_ev_pH * eagle_ev_to_erg / eagle_proton_mass_cgs;
+        cooling->he_reion_ev_pH / cooling->proton_mass_cgs;
     extra_heating += he_reion_erg_pG *
                      (erf((z - dz - cooling->he_reion_z_center) /
                           (pow(2., 0.5) * cooling->he_reion_z_sigma)) -
@@ -196,10 +197,10 @@ __attribute__((always_inline)) INLINE double eagle_metal_cooling_rate(
   if (z > cooling->Redshifts[cooling->N_Redshifts - 1] ||
       z > cooling->reionisation_redshift) {
 
-    T_gam = eagle_cmb_temperature * (1 + z);
+    T_gam = cooling->T_CMB_0 * (1 + z);
 
-    temp_lambda = -eagle_compton_rate *
-                  (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
+    temp_lambda = -cooling->compton_rate_cgs *
+                  (temp - cooling->T_CMB_0 * (1 + z)) * pow((1 + z), 4) *
                   h_plus_he_electron_abundance / n_h;
     cooling_rate += temp_lambda;
     if (element_lambda != NULL) element_lambda[1] = temp_lambda;
@@ -284,130 +285,6 @@ __attribute__((always_inline)) INLINE double eagle_metal_cooling_rate(
   return cooling_rate;
 }
 
-/*
- * @brief Calculates cooling rate by interpolating precomputed 1d tables
- * of cooling rates which depend only on (log base 10 of) temperature.
- * May be used in bisection scheme when need to perform many table
- * interpolations.
- *
- * @param log_10_u Log base 10 of internal energy
- * @param dlambda_du Pointer to value to be set to rate
- * @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
- * to H,He
- * @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due
- * to metals
- * @param temp_table 1D table of temperatures
- * of change of cooling rate with respect to internal energy
- * @param z_index Redshift index of particle
- * @param dz Redshift offset of particle
- * @param n_h_i Particle hydrogen number density index
- * @param d_n_h Particle hydrogen number density offset
- * @param He_i Particle helium fraction index
- * @param d_He Particle helium fraction offset
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param element_lambda Pointer to array for printing contribution
- * to cooling rate from each of the metals.
- * @param solar_ratio Array of ratios of particle metal abundances
- * to solar metal abundances
- */
-__attribute__((always_inline)) INLINE double eagle_metal_cooling_rate_1d_table(
-    double log_10_u, double *dlambda_du, double *H_plus_He_heat_table,
-    double *H_plus_He_electron_abundance_table, double *element_cooling_table,
-    double *element_electron_abundance_table, double *temp_table,
-    const struct part *restrict p,
-    const struct cooling_function_data *restrict cooling,
-    const struct cosmology *restrict cosmo,
-    const struct phys_const *internal_const, double *element_lambda) {
-  double T_gam, solar_electron_abundance;
-  double n_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
-                                            internal_const) *
-               cooling->number_density_scale;  // chemistry_data
-  double z = cosmo->z;
-  double cooling_rate = 0.0, temp_lambda;
-  float du;
-  float h_plus_he_electron_abundance;
-
-  int i, temp_i;
-  double temp;
-  float d_temp;
-
-  *dlambda_du = 0.0;
-
-  // interpolate to get temperature of particles, find where we are in
-  // the temperature table.
-  temp = eagle_convert_u_to_temp_1d_table(log_10_u, &du, temp_table, cooling);
-  get_index_1d(cooling->Temp, cooling->N_Temp, temp, &temp_i, &d_temp);
-
-  /* ------------------ */
-  /* Metal-free cooling */
-  /* ------------------ */
-
-  // contribution to cooling and electron abundance from H, He.
-  temp_lambda = interpol_1d_dbl(H_plus_He_heat_table, temp_i, d_temp);
-  h_plus_he_electron_abundance =
-      interpol_1d_dbl(H_plus_He_electron_abundance_table, temp_i, d_temp);
-  cooling_rate += temp_lambda;
-  *dlambda_du += (((double)H_plus_He_heat_table[temp_i + 1]) -
-                  ((double)H_plus_He_heat_table[temp_i])) /
-                 ((double)du);
-  if (element_lambda != NULL) element_lambda[0] = temp_lambda;
-
-  /* ------------------ */
-  /* Compton cooling    */
-  /* ------------------ */
-
-  // inverse Compton cooling is not in collisional table
-  // before reionisation so add now
-
-  if (z > cooling->Redshifts[cooling->N_Redshifts - 1] ||
-      z > cooling->reionisation_redshift) {
-    /* inverse Compton cooling is not in collisional table
-       before reionisation so add now */
-
-    T_gam = eagle_cmb_temperature * (1 + z);
-
-    temp_lambda = -eagle_compton_rate *
-                  (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
-                  h_plus_he_electron_abundance / n_h;
-    cooling_rate += temp_lambda;
-    if (element_lambda != NULL) element_lambda[1] = temp_lambda;
-  }
-
-  /* ------------- */
-  /* Metal cooling */
-  /* ------------- */
-
-  // for each element the cooling rate is multiplied by the ratio of H, He
-  // electron abundance to solar electron abundance. Note: already multiplied by
-  // ratio of particle to solar metal abundances when creating 1d tables.
-
-  solar_electron_abundance =
-      interpol_1d_dbl(element_electron_abundance_table, temp_i, d_temp);
-
-  for (i = 0; i < cooling->N_Elements; i++) {
-    temp_lambda = interpol_1d_dbl(element_cooling_table + i * cooling->N_Temp,
-                                  temp_i, d_temp) *
-                  (h_plus_he_electron_abundance / solar_electron_abundance);
-    cooling_rate += temp_lambda;
-    *dlambda_du +=
-        (((double)element_cooling_table[i * cooling->N_Temp + temp_i + 1]) *
-             ((double)H_plus_He_electron_abundance_table[temp_i + 1]) /
-             ((double)element_electron_abundance_table[temp_i + 1]) -
-         ((double)element_cooling_table[i * cooling->N_Temp + temp_i]) *
-             ((double)H_plus_He_electron_abundance_table[temp_i]) /
-             ((double)element_electron_abundance_table[temp_i])) /
-        ((double)du);
-    if (element_lambda != NULL) element_lambda[i + 2] = temp_lambda;
-  }
-
-  return cooling_rate;
-}
-
 /*
  * @brief Wrapper function used to calculate cooling rate and dLambda_du.
  * Table indices and offsets for redshift, hydrogen number density and
@@ -441,7 +318,7 @@ __attribute__((always_inline)) INLINE double eagle_cooling_rate(
   double lambda_net = 0.0;
 
   lambda_net = eagle_metal_cooling_rate(
-      logu / eagle_log_10, dLambdaNet_du, z_index, dz, n_h_i, d_n_h, He_i, d_He,
+      logu / M_LN10, dLambdaNet_du, z_index, dz, n_h_i, d_n_h, He_i, d_He,
       p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
   if (isnan(lambda_net))
     printf(
@@ -452,53 +329,6 @@ __attribute__((always_inline)) INLINE double eagle_cooling_rate(
   return lambda_net;
 }
 
-/*
- * @brief Wrapper function used to calculate cooling rate and dLambda_du.
- * Table indices and offsets for redshift, hydrogen number density and
- * helium fraction are passed it so as to compute them only once per particle.
- *
- * @param logu Natural log of internal energy
- * @param dlambda_du Pointer to gradient of cooling rate (set in
- * eagle_metal_cooling_rate)
- * @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
- * to H,He
- * @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due
- * to metals
- * @param temp_table 1D table of temperatures
- * @param z_index Redshift index of particle
- * @param dz Redshift offset of particle
- * @param n_h_i Particle hydrogen number density index
- * @param d_n_h Particle hydrogen number density offset
- * @param He_i Particle helium fraction index
- * @param d_He Particle helium fraction offset
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param abundance_ratio Ratio of element abundance to solar
- */
-__attribute__((always_inline)) INLINE double eagle_cooling_rate_1d_table(
-    double logu, double *dLambdaNet_du, double *H_plus_He_heat_table,
-    double *H_plus_He_electron_abundance_table, double *element_cooling_table,
-    double *element_electron_abundance_table, double *temp_table, int z_index,
-    float dz, int n_h_i, float d_n_h, int He_i, float d_He,
-    const struct part *restrict p,
-    const struct cooling_function_data *restrict cooling,
-    const struct cosmology *restrict cosmo,
-    const struct phys_const *internal_const, float *abundance_ratio) {
-  double *element_lambda = NULL, lambda_net = 0.0;
-
-  lambda_net = eagle_metal_cooling_rate_1d_table(
-      logu / eagle_log_10, dLambdaNet_du, H_plus_He_heat_table,
-      H_plus_He_electron_abundance_table, element_cooling_table,
-      element_electron_abundance_table, temp_table, p, cooling, cosmo,
-      internal_const, element_lambda);
-
-  return lambda_net;
-}
-
 /*
  * @brief Wrapper function used to calculate cooling rate and dLambda_du.
  * Writes to file contribution from each element to cooling rate for testing
@@ -530,9 +360,18 @@ __attribute__((always_inline)) INLINE double eagle_print_metal_cooling_rate(
 
   char output_filename[25];
   FILE **output_file = malloc((cooling->N_Elements + 2) * sizeof(FILE *));
+  print_cooling_rate_contribution_flag += 1.0/(cooling->N_Elements+2);
   for (int element = 0; element < cooling->N_Elements + 2; element++) {
     sprintf(output_filename, "%s%d%s", "cooling_element_", element, ".dat");
-    output_file[element] = fopen(output_filename, "a");
+    if (print_cooling_rate_contribution_flag < 1) {
+      // If this is the first time we're running this function, overwrite the
+      // output files
+      output_file[element] = fopen(output_filename, "w");
+      print_cooling_rate_contribution_flag += 1.0/(cooling->N_Elements+2);
+    } else {
+      // append to existing files
+      output_file[element] = fopen(output_filename, "a");
+    }
     if (output_file == NULL) {
       printf("Error opening file!\n");
       exit(1);
@@ -552,66 +391,6 @@ __attribute__((always_inline)) INLINE double eagle_print_metal_cooling_rate(
   return lambda_net;
 }
 
-/*
- * @brief Wrapper function used to calculate cooling rate and dLambda_du.
- * Prints contribution from each element to cooling rate for testing purposes.
- * Table indices and offsets for redshift, hydrogen number density and
- * helium fraction are passed it so as to compute them only once per particle.
- *
- * @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
- * to H,He
- * @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due
- * to metals
- * @param temp_table 1D table of temperatures
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param abundance_ratio Ratio of element abundance to solar
- */
-__attribute__((always_inline)) INLINE double
-eagle_print_metal_cooling_rate_1d_table(
-    double *H_plus_He_heat_table, double *H_plus_He_electron_abundance_table,
-    double *element_print_cooling_table,
-    double *element_electron_abundance_table, double *temp_table,
-    const struct part *restrict p,
-    const struct cooling_function_data *restrict cooling,
-    const struct cosmology *restrict cosmo,
-    const struct phys_const *internal_const) {
-  double *element_lambda, lambda_net = 0.0;
-  element_lambda = malloc((cooling->N_Elements + 2) * sizeof(double));
-  double u = hydro_get_physical_internal_energy(p, cosmo) *
-             cooling->internal_energy_scale;
-  double dLambdaNet_du;
-
-  char output_filename[25];
-  FILE **output_file = malloc((cooling->N_Elements + 2) * sizeof(FILE *));
-  for (int element = 0; element < cooling->N_Elements + 2; element++) {
-    sprintf(output_filename, "%s%d%s", "cooling_element_", element, ".dat");
-    output_file[element] = fopen(output_filename, "a");
-    if (output_file == NULL) {
-      printf("Error opening file!\n");
-      exit(1);
-    }
-  }
-
-  for (int j = 0; j < cooling->N_Elements + 2; j++) element_lambda[j] = 0.0;
-  lambda_net = eagle_metal_cooling_rate_1d_table(
-      log10(u), &dLambdaNet_du, H_plus_He_heat_table,
-      H_plus_He_electron_abundance_table, element_print_cooling_table,
-      element_electron_abundance_table, temp_table, p, cooling, cosmo,
-      internal_const, element_lambda);
-  for (int j = 0; j < cooling->N_Elements + 2; j++) {
-    fprintf(output_file[j], "%.5e\n", element_lambda[j]);
-  }
-
-  for (int i = 0; i < cooling->N_Elements + 2; i++) fclose(output_file[i]);
-
-  return lambda_net;
-}
-
 /*
  * @brief Bisection integration scheme used when Newton Raphson method fails to
  * converge
@@ -650,7 +429,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
 
   /* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
    * equivalent expression  below */
-  double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+  double ratefact = inn_h * (XH / cooling->proton_mass_cgs);
 
   // Bracketing
   u_lower = u_init;
@@ -760,14 +539,14 @@ __attribute__((always_inline)) INLINE float newton_iter(
 
   /* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
    * equivalent expression  below */
-  double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+  double ratefact = inn_h * (XH / cooling->proton_mass_cgs);
 
   logu_old = logu_init;
   logu = logu_old;
   int i = 0;
   float log_table_bound_high =
-      (cooling->Therm[cooling->N_Temp - 1] + 1) / eagle_log_10_e;
-  float log_table_bound_low = (cooling->Therm[0] + 1) / eagle_log_10_e;
+      (cooling->Therm[cooling->N_Temp - 1] + 1) / M_LOG10E;
+  float log_table_bound_low = (cooling->Therm[0] + 1) / M_LOG10E;
 
   do /* iterate to convergence */
   {
@@ -891,10 +670,10 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
 
   /* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
    * equivalent expression  below */
-  ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+  ratefact = inn_h * (XH / cooling->proton_mass_cgs);
 
   /* set helium and hydrogen reheating term */
-  if (cooling->he_reion == 1)
+  if (cooling->he_reion_flag == 1)
     LambdaTune = eagle_helium_reionization_extraheat(z_index, dz, cooling);
 
   // compute hydrogen number density and helium fraction table indices
@@ -1011,70 +790,6 @@ __attribute__((always_inline)) INLINE float cooling_get_radiated_energy(
   return xp->cooling_data.radiated_energy;
 }
 
-/**
- * To do: incorporate into code so that not all tables are read at the beginning
- * @brief Checks the tables that are currently loaded in memory and read
- * new ones if necessary.
- *
- * @param cooling The #cooling_function_data we play with.
- * @param index_z The index along the redshift axis of the tables of the current
- * z.
- */
-inline void eagle_check_cooling_tables(struct cooling_function_data *cooling,
-                                       int index_z) {
-
-  /* Do we already have the right table in memory? */
-  if (cooling->low_z_index == index_z) return;
-
-  if (index_z >= cooling->N_Redshifts) {
-
-    error("Missing implementation");
-    // Add reading of high-z tables
-
-    /* Record the table indices */
-    cooling->low_z_index = index_z;
-    cooling->high_z_index = index_z;
-  } else {
-
-    /* Record the table indices */
-    cooling->low_z_index = index_z;
-    cooling->high_z_index = index_z + 1;
-
-    /* Load the damn thing */
-    eagle_read_two_tables(cooling);
-  }
-}
-
-/**
- * To do: incorporate into code so that not all tables are read at the beginning
- * @brief Common operations performed on the cooling function at a
- * given time-step or redshift.
- *
- * Here we load the cooling tables corresponding to our redshift.
- *
- * @param phys_const The physical constants in internal units.
- * @param us The internal system of units.
- * @param cosmo The current cosmological model.
- * @param cooling The #cooling_function_data used in the run.
- */
-INLINE void cooling_update(const struct phys_const *phys_const,
-                           const struct unit_system *us,
-                           const struct cosmology *cosmo,
-                           struct cooling_function_data *cooling) {
-  /* Current redshift */
-  const float redshift = cosmo->z;
-
-  /* Get index along the redshift index of the tables */
-  int index_z;
-  float delta_z_table;
-  get_redshift_index(redshift, &index_z, &delta_z_table, cooling);
-  cooling->index_z = index_z;
-  cooling->delta_z_table = delta_z_table;
-
-  /* Load the current table (if different from what we have) */
-  eagle_check_cooling_tables(cooling, index_z);
-}
-
 /**
  * @brief Initialises the cooling properties.
  *
@@ -1099,7 +814,7 @@ INLINE void cooling_init_backend(struct swift_params *parameter_file,
       parameter_file, "EAGLEChemistry:CalciumOverSilicon");
   cooling->sulphur_over_silicon_ratio = parser_get_param_float(
       parameter_file, "EAGLEChemistry:SulphurOverSilicon");
-  cooling->he_reion =
+  cooling->he_reion_flag =
       parser_get_param_float(parameter_file, "EagleCooling:he_reion");
   cooling->he_reion_z_center =
       parser_get_param_float(parameter_file, "EagleCooling:he_reion_z_center");
@@ -1108,15 +823,16 @@ INLINE void cooling_init_backend(struct swift_params *parameter_file,
   cooling->he_reion_ev_pH =
       parser_get_param_float(parameter_file, "EagleCooling:he_reion_ev_pH");
 
+  /* convert to cgs */
+  cooling->he_reion_ev_pH *= phys_const->const_electron_volt *
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
+
   // read in cooling tables
   GetCoolingRedshifts(cooling);
   sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
   ReadCoolingHeader(fname, cooling);
   cooling->table = eagle_readtable(cooling->cooling_table_path, cooling);
 
-  // allocate array for solar abundances
-  cooling->solar_abundances = malloc(cooling->N_Elements * sizeof(double));
-
   // compute conversion factors
   cooling->internal_energy_scale =
       units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) /
@@ -1128,6 +844,31 @@ INLINE void cooling_init_backend(struct swift_params *parameter_file,
                          units_cgs_conversion_factor(us, UNIT_CONV_MASS);
   cooling->temperature_scale =
       units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+  cooling->proton_mass_cgs = phys_const->const_proton_mass *
+      units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+  cooling->T_CMB_0 = phys_const->const_T_CMB_0 *
+      units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+  /* Compute the coefficient at the front of the Compton cooling expression */
+  const double radiation_constant =
+      4. * phys_const->const_stefan_boltzmann / phys_const->const_speed_light_c;
+  const double compton_coefficient =
+      4. * radiation_constant * phys_const->const_thomson_cross_section *
+      phys_const->const_boltzmann_k /
+      (phys_const->const_electron_mass * phys_const->const_speed_light_c);
+  const float dimension_coefficient[5] = {1, 2, -3, 0, -5};
+
+  /* This should be ~1.0178085e-37 g cm^2 s^-3 K^-5 */
+  const double compton_coefficient_cgs =
+      compton_coefficient *
+      units_general_cgs_conversion_factor(us, dimension_coefficient);
+
+  /* And now the Compton rate */
+  cooling->compton_rate_cgs =
+      compton_coefficient_cgs * cooling->T_CMB_0 * 
+      cooling->T_CMB_0 * cooling->T_CMB_0 * cooling->T_CMB_0;
+
 }
 
 /**
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index f1ad66b055861cd80c5ab9b1d99799449ac5ed79..2eb887ae51377e6922699d16c352fce11d6392aa 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -32,15 +32,6 @@
 #include "physical_constants.h"
 #include "units.h"
 
-enum hdf5_allowed_types {
-  hdf5_short,
-  hdf5_int,
-  hdf5_long,
-  hdf5_float,
-  hdf5_double,
-  hdf5_char
-};
-
 double eagle_helium_reionization_extraheat(
     double, double, const struct cooling_function_data *restrict);
 
@@ -50,35 +41,17 @@ double eagle_metal_cooling_rate(double, double *, int, float, int, float, int,
                                 const struct cosmology *restrict,
                                 const struct phys_const *, double *, float *);
 
-double eagle_metal_cooling_rate_1d_table(
-    double, double *, double *, double *, double *, double *, double *,
-    const struct part *restrict, const struct cooling_function_data *restrict,
-    const struct cosmology *restrict, const struct phys_const *, double *);
-
 double eagle_cooling_rate(double, double *, int, float, int, float, int, float,
                           const struct part *restrict,
                           const struct cooling_function_data *restrict,
                           const struct cosmology *restrict,
                           const struct phys_const *, float *);
 
-double eagle_cooling_rate_1d_table(double, double *, double *, double *,
-                                   double *, double *, double *, int, float,
-                                   int, float, int, float,
-                                   const struct part *restrict,
-                                   const struct cooling_function_data *restrict,
-                                   const struct cosmology *restrict,
-                                   const struct phys_const *, float *);
-
 double eagle_print_metal_cooling_rate(
     int, float, int, float, int, float, const struct part *restrict,
     const struct cooling_function_data *restrict,
     const struct cosmology *restrict, const struct phys_const *, float *);
 
-double eagle_print_metal_cooling_rate_1d_table(
-    double *, double *, double *, double *, double *,
-    const struct part *restrict, const struct cooling_function_data *restrict,
-    const struct cosmology *restrict, const struct phys_const *);
-
 float bisection_iter(float, double, int, float, int, float, int, float, float,
                      struct part *restrict, const struct cosmology *restrict,
                      const struct cooling_function_data *restrict,
@@ -116,11 +89,6 @@ void cooling_first_init_part(const struct phys_const *restrict,
 
 float cooling_get_radiated_energy(const struct xpart *restrict);
 
-void eagle_check_cooling_tables(struct cooling_function_data *, int);
-
-void cooling_update(const struct phys_const *, const struct unit_system *,
-                    const struct cosmology *, struct cooling_function_data *);
-
 void cooling_init_backend(struct swift_params *, const struct unit_system *,
                           const struct phys_const *,
                           struct cooling_function_data *);
diff --git a/src/cooling/EAGLE/cooling_read_table.h b/src/cooling/EAGLE/cooling_read_table.h
deleted file mode 100644
index 82bf92bdd54020bce12c9ffc8a42b19b0ea400fa..0000000000000000000000000000000000000000
--- a/src/cooling/EAGLE/cooling_read_table.h
+++ /dev/null
@@ -1,132 +0,0 @@
-#ifndef SWIFT_COOLING_EAGLE_READ_TABLES_H
-#define SWIFT_COOLING_EAGLE_READ_TABLES_H
-
-/*
- * @brief Reads cooling header data into cooling data structure
- *
- * @param fname Filepath
- * @param cooling Cooling data structure
- */
-INLINE void ReadCoolingHeader(char *fname,
-                              struct cooling_function_data *cooling) {
-  int i;
-
-  hid_t tempfile_id, dataset, datatype;
-
-  herr_t status;
-
-  /* fill the constants */
-  tempfile_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-  if (tempfile_id < 0) {
-    error("[ReadCoolingHeader()]: unable to open file %s\n", fname);
-  }
-
-  dataset =
-      H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_Temp);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Number_of_density_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_nH);
-  status = H5Dclose(dataset);
-
-  dataset =
-      H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_He);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_SolarAbundances);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_Elements);
-  status = H5Dclose(dataset);
-
-  /* allocate arrays for cooling table header */
-  // allocate_header_arrays();
-
-  cooling->Temp = malloc(cooling->N_Temp * sizeof(float));
-  cooling->nH = malloc(cooling->N_Temp * sizeof(float));
-  cooling->HeFrac = malloc(cooling->N_Temp * sizeof(float));
-  cooling->SolarAbundances = malloc(cooling->N_Temp * sizeof(float));
-
-  /* fill the arrays */
-  dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->Temp);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Solar/Hydrogen_density_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->nH);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->HeFrac);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->SolarAbundances);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->Therm);
-  status = H5Dclose(dataset);
-
-  char element_names[cooling->N_Elements][element_name_length];
-  hsize_t string_length = element_name_length;
-
-  /* names of chemical elements stored in table */
-  datatype = H5Tcopy(H5T_C_S1);
-  status = H5Tset_size(datatype, string_length);
-  dataset = H5Dopen(tempfile_id, "/Header/Metal_names", H5P_DEFAULT);
-  status =
-      H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
-  status = H5Dclose(dataset);
-
-  for (i = 0; i < cooling->N_Elements; i++)
-    cooling->ElementNames[i] = mystrdup(element_names[i]);
-
-  // char solar_abund_names[cooling_N_SolarAbundances][EL_NAME_LENGTH];
-
-  /* assumed solar abundances used in constructing the tables, and corresponding
-   * names */
-  // dataset = H5Dopen(tempfile_id, "/Header/Abundances/Abund_names",
-  // H5P_DEFAULT);
-  // status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-  //                 solar_abund_names);
-  // status = H5Dclose(dataset);
-  // H5Tclose(datatype);
-
-  // for (i = 0; i < cooling_N_SolarAbundances; i++)
-  //  cooling_SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
-
-  // status = H5Fclose(tempfile_id);
-
-  /* Convert to temperature, density and internal energy arrays to log10 */
-  for (i = 0; i < cooling->N_Temp; i++) {
-    cooling->Temp[i] = log10(cooling->Temp[i]);
-    cooling->Therm[i] = log10(cooling->Therm[i]);
-  }
-
-  for (i = 0; i < cooling->N_nH; i++) cooling->nH[i] = log10(cooling->nH[i]);
-
-  printf("Done with cooling table header.\n");
-  fflush(stdout);
-}
-
-#endif /* SWIFT_COOLING_EAGLE_READ_TABLES_H */
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index da157ad8db558a162bc70c753d59c7e91e9b303c..8b86239e03cb91621d25348067d3ffbb6fa700e7 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/src/cooling/EAGLE/cooling_struct.h
@@ -18,31 +18,11 @@
  ******************************************************************************/
 #ifndef SWIFT_COOLING_STRUCT_EAGLE_H
 #define SWIFT_COOLING_STRUCT_EAGLE_H
-#include <stdbool.h>
-#include <time.h>
 
 // EAGLE defined constants
 #define eagle_element_name_length 20
-#define eagle_cmb_temperature 2.728
-#define eagle_compton_rate 1.0178e-37 * 2.728 * 2.728 * 2.728 * 2.728
 #define eagle_metal_cooling_on 1
 #define eagle_max_iterations 15
-#define eagle_proton_mass_cgs 1.6726e-24
-#define eagle_ev_to_erg 1.60217646e-12
-#define eagle_log_10 2.30258509299404
-#define eagle_log_10_e 0.43429448190325182765
-
-/*
- * @brief struct containing radiative and heating rates
- */
-// struct radiative_rates{
-//  float Cooling[NUMBER_OF_CTYPES];
-//  float Heating[NUMBER_OF_HTYPES];
-//  float TotalCoolingRate;
-//  float TotalHeatingRate;
-//  float CoolingTime; /*Cooling time in Myr*/
-//  float HeatingTime; /*Heating time in Myr*/
-//};
 
 /*
  * @brief struct containing cooling tables independent of redshift
@@ -99,50 +79,59 @@ struct cooling_function_data {
   /* Cooling table variables */
   struct eagle_cooling_table table;
 
+  /* Size of table dimensions */
   int N_Redshifts;
   int N_nH;
   int N_Temp;
   int N_He;
+
+  /* Number of metals and solar abundances tracked in EAGLE */
   int N_Elements;
   int N_SolarAbundances;
 
+  /* Arrays of grid values in tables */
   float *Redshifts;
   float *nH;
   float *Temp;
   float *HeFrac;
   float *Therm;
+
+  /* Array of values of solar metal and electron abundance */
   float *SolarAbundances;
   float *SolarElectronAbundance;
-  float *ElementAbundance_SOLARM1;
-  double *solar_abundances;
+
+
+  /* Arrays containing names of tracked elements. Used for 
+   * reading in the contributions to the cooling rate from
+   * each element */
   char **ElementNames;
   char **SolarAbundanceNames;
-  int *ElementNamePointers;
-  int *SolarAbundanceNamePointers;
-
-  /*! Constant multiplication factor for time-step criterion */
-  float cooling_tstep_mult;
-  float Redshift;
-  float min_energy;
-  float cooling_rate;
+
+  /* Normalisation constants that are frequently used */
   double internal_energy_scale;
   double number_density_scale;
   double temperature_scale;
   double power_scale;
+
+  /* filepath to EAGLE cooling tables */
   char cooling_table_path[500];
+
+  /* Some constants read in from yml file relevant to EAGLE cooling */
   float reionisation_redshift;
   float calcium_over_silicon_ratio;
   float sulphur_over_silicon_ratio;
 
-  int he_reion;
+  /* Helium reionisation parameters */
+  int he_reion_flag;
   float he_reion_ev_pH;
   float he_reion_z_center;
   float he_reion_z_sigma;
 
-  int index_z, low_z_index, high_z_index;
-  float delta_z_table;
+  /* Proton mass in cgs */
+  double proton_mass_cgs;
+  double T_CMB_0;
+  double compton_rate_cgs;
 
-  double delta_u;
 };
 
 /**
diff --git a/src/cooling/EAGLE/eagle_cool_tables.c b/src/cooling/EAGLE/eagle_cool_tables.c
index 9b4bd6b17ba43c2112368afd30e35ff16dfee76f..4aa8d8496095ebd703d4ecf5057c7037223aca25 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.c
+++ b/src/cooling/EAGLE/eagle_cool_tables.c
@@ -36,19 +36,6 @@
 #include "error.h"
 #include "interpolate.h"
 
-/*
- * @brief String assignment function from EAGLE
- *
- * @param s String
- */
-char *mystrdup(const char *s) {
-  char *p;
-
-  p = (char *)malloc((strlen(s) + 1) * sizeof(char));
-  strcpy(p, s);
-  return p;
-}
-
 /*
  * @brief Reads in EAGLE table of redshift values
  *
@@ -134,10 +121,14 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
   cooling->HeFrac = malloc(cooling->N_He * sizeof(float));
   cooling->SolarAbundances = malloc(cooling->N_SolarAbundances * sizeof(float));
   cooling->Therm = malloc(cooling->N_Temp * sizeof(float));
-  cooling->ElementNames =
-      malloc(cooling->N_Elements * eagle_element_name_length * sizeof(char));
-  cooling->SolarAbundanceNames = malloc(
-      cooling->N_SolarAbundances * eagle_element_name_length * sizeof(char));
+
+  cooling->ElementNames = malloc(cooling->N_Elements * sizeof(char *));
+  for (i = 0; i < cooling->N_Elements; i++)
+    cooling->ElementNames[i] = malloc(eagle_element_name_length * sizeof(char));
+
+  cooling->SolarAbundanceNames = malloc(cooling->N_SolarAbundances * sizeof(char *));
+  for (i = 0; i < cooling->N_SolarAbundances; i++)
+    cooling->SolarAbundanceNames[i] = malloc(eagle_element_name_length * sizeof(char));
 
   // read in values for each of the arrays
   dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
@@ -181,7 +172,8 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
   H5Tclose(datatype);
 
   for (i = 0; i < cooling->N_Elements; i++)
-    cooling->ElementNames[i] = mystrdup(element_names[i]);
+    strcpy(cooling->ElementNames[i],element_names[i]);
+  
 
   char solar_abund_names[cooling->N_SolarAbundances][eagle_element_name_length];
 
@@ -196,7 +188,7 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
   H5Tclose(datatype);
 
   for (i = 0; i < cooling->N_SolarAbundances; i++)
-    cooling->SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
+    strcpy(cooling->SolarAbundanceNames[i],element_names[i]);
 
   status = H5Fclose(tempfile_id);
 
@@ -355,7 +347,7 @@ struct cooling_tables_redshift_invariant get_no_compt_table(
   free(he_net_cooling_rate);
   free(he_electron_abundance);
 
-  printf("eagle_cool_tables.h done reading in no compton table\n");
+  message("done reading in no compton table");
 
   return cooling_table;
 }
@@ -490,7 +482,7 @@ struct cooling_tables_redshift_invariant get_collisional_table(
   free(he_net_cooling_rate);
   free(he_electron_abundance);
 
-  printf("eagle_cool_tables.h done reading in collisional table\n");
+  message("done reading in collisional table");
   return cooling_table;
 }
 
@@ -639,7 +631,7 @@ struct cooling_tables_redshift_invariant get_photodis_table(
   free(he_net_cooling_rate);
   free(he_electron_abundance);
 
-  printf("eagle_cool_tables.h done reading in photodissociation table\n");
+  message("done reading in photodissociation table");
   return cooling_table;
 }
 
@@ -793,7 +785,7 @@ struct cooling_tables get_cooling_table(
   free(he_net_cooling_rate);
   free(he_electron_abundance);
 
-  printf("eagle_cool_tables.h done reading in general cooling table\n");
+  message("done reading in general cooling table");
 
   return cooling_table;
 }
@@ -820,179 +812,4 @@ struct eagle_cooling_table eagle_readtable(
   return table;
 }
 
-/*
- * @brief Work in progress: read in only two tables for given redshift
- *
- * @param cooling_table_path Directory containing cooling tables
- * @param cooling Cooling function data structure
- * @param filename1 filename2 Names of the two tables we need to read
- */
-struct cooling_tables get_two_cooling_tables(
-    char *cooling_table_path,
-    const struct cooling_function_data *restrict cooling, char *filename1,
-    char *filename2) {
-
-  struct cooling_tables cooling_table;
-  hid_t file_id, dataset;
-
-  herr_t status;
-
-  char fname[500], set_name[500];
-
-  int specs, i, j, k, table_index, cooling_index;
-
-  float *net_cooling_rate;
-  float *electron_abundance;
-  float *temperature;
-  float *he_net_cooling_rate;
-  float *he_electron_abundance;
-
-  // allocate temporary arrays (needed to change order of dimensions
-  // of arrays)
-  net_cooling_rate =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  electron_abundance =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                cooling->N_nH * sizeof(float));
-  he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                        cooling->N_nH * sizeof(float));
-  he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                          cooling->N_nH * sizeof(float));
-
-  // allocate arrays that the values will be assigned to
-  cooling_table.metal_heating =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_Elements *
-                      cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.electron_abundance = (float *)malloc(
-      cooling->N_Redshifts * cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.temperature =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_heating =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_electron_abundance =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-
-  // repeat for both tables we need to read
-  for (int file = 0; file < 2; file++) {
-    if (file == 0) {
-      sprintf(fname, "%s%s.hdf5", cooling_table_path, filename1);
-    } else {
-      sprintf(fname, "%s%s.hdf5", cooling_table_path, filename2);
-    }
-    file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-    if (file_id < 0) {
-      error("[GetCoolingTables()]: unable to open file %s\n", fname);
-    }
-
-    // read in cooling rates due to metals
-    for (specs = 0; specs < cooling->N_Elements; specs++) {
-      sprintf(set_name, "/%s/Net_Cooling", cooling->ElementNames[specs]);
-      dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-      status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                       net_cooling_rate);
-      status = H5Dclose(dataset);
-
-      for (i = 0; i < cooling->N_nH; i++) {
-        for (j = 0; j < cooling->N_Temp; j++) {
-          table_index =
-              row_major_index_2d(j, i, cooling->N_Temp, cooling->N_nH);
-          cooling_index = row_major_index_4d(
-              file, i, j, specs, cooling->N_Redshifts, cooling->N_nH,
-              cooling->N_Temp, cooling->N_Elements);
-          cooling_table.metal_heating[cooling_index] =
-              -net_cooling_rate[table_index];
-        }
-      }
-    }
-
-    // read in cooling rates due to hydrogen and helium, H + He electron
-    // abundances, temperatures
-    sprintf(set_name, "/Metal_free/Net_Cooling");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     he_net_cooling_rate);
-    status = H5Dclose(dataset);
-
-    sprintf(set_name, "/Metal_free/Temperature/Temperature");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     temperature);
-    status = H5Dclose(dataset);
-
-    sprintf(set_name, "/Metal_free/Electron_density_over_n_h");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     he_electron_abundance);
-    status = H5Dclose(dataset);
-
-    for (i = 0; i < cooling->N_He; i++) {
-      for (j = 0; j < cooling->N_Temp; j++) {
-        for (k = 0; k < cooling->N_nH; k++) {
-          table_index = row_major_index_3d(i, j, k, cooling->N_He,
-                                           cooling->N_Temp, cooling->N_nH);
-          cooling_index =
-              row_major_index_4d(file, k, i, j, cooling->N_Redshifts,
-                                 cooling->N_nH, cooling->N_He, cooling->N_Temp);
-          cooling_table.H_plus_He_heating[cooling_index] =
-              -he_net_cooling_rate[table_index];
-          cooling_table.H_plus_He_electron_abundance[cooling_index] =
-              he_electron_abundance[table_index];
-          cooling_table.temperature[cooling_index] =
-              log10(temperature[table_index]);
-        }
-      }
-    }
-
-    // read in electron densities due to metals
-    sprintf(set_name, "/Solar/Electron_density_over_n_h");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     electron_abundance);
-    status = H5Dclose(dataset);
-
-    for (i = 0; i < cooling->N_Temp; i++) {
-      for (j = 0; j < cooling->N_nH; j++) {
-        table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
-        cooling_index = row_major_index_3d(file, j, i, cooling->N_Redshifts,
-                                           cooling->N_nH, cooling->N_Temp);
-        cooling_table.electron_abundance[cooling_index] =
-            electron_abundance[table_index];
-      }
-    }
-
-    status = H5Fclose(file_id);
-  }
-
-  free(net_cooling_rate);
-  free(electron_abundance);
-  free(temperature);
-  free(he_net_cooling_rate);
-  free(he_electron_abundance);
-
-  printf("eagle_cool_tables.h done reading in general cooling table\n");
-
-  return cooling_table;
-}
-
-/*
- * @brief Work in progress: constructs the data structure containting cooling
- * tables bracketing the redshift of a particle.
- *
- * @param cooling_table_path Directory containing cooling table
- * @param cooling Cooling data structure
- */
-void eagle_read_two_tables(struct cooling_function_data *restrict cooling) {
-
-  struct eagle_cooling_table table;
-
-  table.element_cooling =
-      get_cooling_table(cooling->cooling_table_path, cooling);
-  cooling->table = table;
-}
-
 #endif
diff --git a/src/cooling/EAGLE/eagle_cool_tables.h b/src/cooling/EAGLE/eagle_cool_tables.h
index 760ed59487f549e026cd6a63ae533a785ba868b1..3899d5396b49efad5a1361bab41c2d156ba0d7e5 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.h
+++ b/src/cooling/EAGLE/eagle_cool_tables.h
@@ -1,1020 +1,57 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
 #ifndef SWIFT_EAGLE_COOL_TABLES_H
 #define SWIFT_EAGLE_COOL_TABLES_H
 
+/**
+ * @file src/cooling/EAGLE/cooling.h
+ * @brief EAGLE cooling function
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
 #include <hdf5.h>
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
 #include "cooling_struct.h"
+#include "interpolate.h"
 #include "error.h"
 
-/*
- * @brief String assignment function from EAGLE
- *
- * @param s String
- */
-inline char *mystrdup(const char *s) {
-  char *p;
-
-  p = (char *)malloc((strlen(s) + 1) * sizeof(char));
-  strcpy(p, s);
-  return p;
-}
-
-int row_major_index_2d(int, int, int, int);
-int row_major_index_3d(int, int, int, int, int, int);
-int row_major_index_4d(int, int, int, int, int, int, int, int);
-
-/*
- * @brief Reads in EAGLE table redshift values
- *
- * @param cooling Cooling data structure
- */
-inline void GetCoolingRedshifts(struct cooling_function_data *cooling) {
-  FILE *infile;
-
-  int i = 0;
-
-  char buffer[500], redfilename[500];
-
-  sprintf(redfilename, "%s/redshifts.dat", cooling->cooling_table_path);
-  infile = fopen(redfilename, "r");
-  if (infile == NULL) puts("GetCoolingRedshifts can't open a file");
-
-  if (fscanf(infile, "%s", buffer) != EOF) {
-    cooling->N_Redshifts = atoi(buffer);
-    cooling->Redshifts = (float *)malloc(cooling->N_Redshifts * sizeof(float));
-
-    while (fscanf(infile, "%s", buffer) != EOF) {
-      cooling->Redshifts[i] = atof(buffer);
-      i += 1;
-    }
-  }
-  fclose(infile);
-}
-
-/*
- * @brief Reads in EAGLE cooling table header
- *
- * @param fname Filepath
- * @param cooling Cooling data structure
- */
-inline void ReadCoolingHeader(char *fname,
-                              struct cooling_function_data *cooling) {
-  int i;
-
-  hid_t tempfile_id, dataset, datatype;
-
-  herr_t status;
-
-  /* fill the constants */
-  tempfile_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-  if (tempfile_id < 0) {
-    error("[ReadCoolingHeader()]: unable to open file %s\n", fname);
-  }
-
-  dataset =
-      H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_Temp);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Number_of_density_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_nH);
-  status = H5Dclose(dataset);
-
-  dataset =
-      H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_He);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_SolarAbundances);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   &cooling->N_Elements);
-  status = H5Dclose(dataset);
-
-  /* allocate arrays for cooling table header */
-  // allocate_header_arrays();
-
-  cooling->Temp = malloc(cooling->N_Temp * sizeof(float));
-  cooling->nH = malloc(cooling->N_nH * sizeof(float));
-  cooling->HeFrac = malloc(cooling->N_He * sizeof(float));
-  cooling->SolarAbundances = malloc(cooling->N_SolarAbundances * sizeof(float));
-  cooling->Therm = malloc(cooling->N_Temp * sizeof(float));
-  cooling->ElementNames =
-      malloc(cooling->N_Elements * eagle_element_name_length * sizeof(char));
-  cooling->SolarAbundanceNames = malloc(
-      cooling->N_SolarAbundances * eagle_element_name_length * sizeof(char));
-
-  /* fill the arrays */
-  dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->Temp);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Solar/Hydrogen_density_bins", H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->nH);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->HeFrac);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->SolarAbundances);
-  status = H5Dclose(dataset);
-
-  dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins",
-                    H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   cooling->Therm);
-  status = H5Dclose(dataset);
-
-  char element_names[cooling->N_Elements][eagle_element_name_length];
-  hsize_t string_length = eagle_element_name_length;
-
-  /* names of chemical elements stored in table */
-  datatype = H5Tcopy(H5T_C_S1);
-  status = H5Tset_size(datatype, string_length);
-  dataset = H5Dopen(tempfile_id, "/Header/Metal_names", H5P_DEFAULT);
-  status =
-      H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
-  status = H5Dclose(dataset);
-  H5Tclose(datatype);
-
-  for (i = 0; i < cooling->N_Elements; i++)
-    cooling->ElementNames[i] = mystrdup(element_names[i]);
-
-  char solar_abund_names[cooling->N_SolarAbundances][eagle_element_name_length];
-
-  /* assumed solar abundances used in constructing the tables, and corresponding
-   * names */
-  datatype = H5Tcopy(H5T_C_S1);
-  status = H5Tset_size(datatype, string_length);
-  dataset = H5Dopen(tempfile_id, "/Header/Abundances/Abund_names", H5P_DEFAULT);
-  status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   solar_abund_names);
-  status = H5Dclose(dataset);
-  H5Tclose(datatype);
-
-  for (i = 0; i < cooling->N_SolarAbundances; i++)
-    cooling->SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
-
-  status = H5Fclose(tempfile_id);
-
-  /* Convert to temperature, density and internal energy arrays to log10 */
-  for (i = 0; i < cooling->N_Temp; i++) {
-    cooling->Temp[i] = log10(cooling->Temp[i]);
-    cooling->Therm[i] = log10(cooling->Therm[i]);
-  }
-
-  for (i = 0; i < cooling->N_nH; i++) cooling->nH[i] = log10(cooling->nH[i]);
-
-  printf("Done with cooling table header.\n");
-  fflush(stdout);
-}
-
-/*
- * @brief Get the cooling table for photoionized cooling (before redshift ~9)
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-inline struct cooling_tables_redshift_invariant get_no_compt_table(
-    char *cooling_table_path,
-    const struct cooling_function_data *restrict cooling) {
-
-  struct cooling_tables_redshift_invariant cooling_table;
-  hid_t file_id, dataset;
-
-  herr_t status;
-
-  char fname[500], set_name[500];
-
-  int specs, i, j, k, table_index, cooling_index;
-
-  float *net_cooling_rate;
-  float *electron_abundance;
-  float *temperature;
-  float *he_net_cooling_rate;
-  float *he_electron_abundance;
-
-  net_cooling_rate =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  electron_abundance =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                cooling->N_nH * sizeof(float));
-  he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                        cooling->N_nH * sizeof(float));
-  he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                          cooling->N_nH * sizeof(float));
-
-  cooling_table.metal_heating = (float *)malloc(
-      cooling->N_Elements * cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.electron_abundance =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                              cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_heating = (float *)malloc(
-      cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_electron_abundance = (float *)malloc(
-      cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
-
-  sprintf(fname, "%sz_8.989nocompton.hdf5", cooling_table_path);
-
-  file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-  // printf("GetNoCompTable Redshift 1 %ld %s\n", (long int)file_id, fname);
-  // fflush(stdout);
-
-  /* For normal elements */
-  for (specs = 0; specs < cooling->N_Elements; specs++) {
-    sprintf(set_name, "/%s/Net_Cooling", cooling->ElementNames[specs]);
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     net_cooling_rate);
-    status = H5Dclose(dataset);
-
-    for (j = 0; j < cooling->N_Temp; j++) {
-      for (k = 0; k < cooling->N_nH; k++) {
-        table_index = row_major_index_2d(j, k, cooling->N_Temp, cooling->N_nH);
-        cooling_index = row_major_index_4d(
-            0, k, j, specs, 1, cooling->N_nH, cooling->N_Temp,
-            cooling->N_Elements);  // Redshift invariant table!!!
-        cooling_table.metal_heating[cooling_index] =
-            -net_cooling_rate[table_index];
-      }
-    }
-  }
-
-  /* Helium */
-  sprintf(set_name, "/Metal_free/Net_Cooling");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   he_net_cooling_rate);
-  status = H5Dclose(dataset);
-
-  sprintf(set_name, "/Metal_free/Temperature/Temperature");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   temperature);
-  status = H5Dclose(dataset);
-
-  sprintf(set_name, "/Metal_free/Electron_density_over_n_h");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   he_electron_abundance);
-  status = H5Dclose(dataset);
-
-  for (i = 0; i < cooling->N_He; i++) {
-    for (j = 0; j < cooling->N_Temp; j++) {
-      for (k = 0; k < cooling->N_nH; k++) {
-        table_index = row_major_index_3d(i, j, k, cooling->N_He,
-                                         cooling->N_Temp, cooling->N_nH);
-        cooling_index =
-            row_major_index_4d(0, k, i, j, 1, cooling->N_nH, cooling->N_He,
-                               cooling->N_Temp);  // Redshift invariant table!!!
-        cooling_table.H_plus_He_heating[cooling_index] =
-            -he_net_cooling_rate[table_index];
-        cooling_table.H_plus_He_electron_abundance[cooling_index] =
-            he_electron_abundance[table_index];
-        cooling_table.temperature[cooling_index] =
-            log10(temperature[table_index]);
-      }
-    }
-  }
-
-  sprintf(set_name, "/Solar/Electron_density_over_n_h");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   electron_abundance);
-  status = H5Dclose(dataset);
-
-  for (i = 0; i < cooling->N_Temp; i++) {
-    for (j = 0; j < cooling->N_nH; j++) {
-      table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
-      cooling_index =
-          row_major_index_3d(0, j, i, 1, cooling->N_nH,
-                             cooling->N_Temp);  // Redshift invariant table!!!
-      cooling_table.electron_abundance[cooling_index] =
-          electron_abundance[table_index];
-    }
-  }
-
-  status = H5Fclose(file_id);
-
-  // cooling_table.metals_heating = cooling_MetalsNetHeating;
-  // cooling_table.H_plus_He_heating = cooling_HplusHeNetHeating;
-  // cooling_table.H_plus_He_electron_abundance =
-  // cooling_HplusHeElectronAbundance;
-  // cooling_table.temperature = cooling_ThermalToTemp;
-  // cooling_table.electron_abundance = cooling_SolarElectronAbundance;
-
-  free(net_cooling_rate);
-  free(electron_abundance);
-  free(temperature);
-  free(he_net_cooling_rate);
-  free(he_electron_abundance);
-
-  printf("eagle_cool_tables.h done reading in no compton table\n");
-
-  return cooling_table;
-}
-
-/*
- * @brief Get the cooling table for collisional cooling (before reionisation)
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-inline struct cooling_tables_redshift_invariant get_collisional_table(
-    char *cooling_table_path,
-    const struct cooling_function_data *restrict cooling) {
-
-  struct cooling_tables_redshift_invariant cooling_table;
-  hid_t file_id, dataset;
-
-  herr_t status;
-
-  char fname[500], set_name[500];
-
-  int specs, i, j, table_index, cooling_index;
-
-  float *net_cooling_rate;
-  float *electron_abundance;
-  float *temperature;
-  float *he_net_cooling_rate;
-  float *he_electron_abundance;
-
-  net_cooling_rate = (float *)malloc(cooling->N_Temp * sizeof(float));
-  electron_abundance = (float *)malloc(cooling->N_Temp * sizeof(float));
-  temperature =
-      (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
-  he_net_cooling_rate =
-      (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
-  he_electron_abundance =
-      (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
-
-  cooling_table.metal_heating =
-      (float *)malloc(cooling->N_Elements * cooling->N_Temp * sizeof(float));
-  cooling_table.electron_abundance =
-      (float *)malloc(cooling->N_Temp * sizeof(float));
-  cooling_table.temperature =
-      (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
-  cooling_table.H_plus_He_heating =
-      (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
-  cooling_table.H_plus_He_electron_abundance =
-      (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
-
-  sprintf(fname, "%sz_collis.hdf5", cooling_table_path);
-
-  file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-  /* For normal elements */
-  for (specs = 0; specs < cooling->N_Elements; specs++) {
-    sprintf(set_name, "/%s/Net_Cooling", cooling->ElementNames[specs]);
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     net_cooling_rate);
-    status = H5Dclose(dataset);
-
-    for (j = 0; j < cooling->N_Temp; j++) {
-      cooling_index =
-          row_major_index_3d(0, specs, j, 1, cooling->N_Elements,
-                             cooling->N_Temp);  // Redshift invariant table!!!
-      cooling_table.metal_heating[cooling_index] = -net_cooling_rate[j];
-    }
-  }
-
-  /* Helium */
-  sprintf(set_name, "/Metal_free/Net_Cooling");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   he_net_cooling_rate);
-  status = H5Dclose(dataset);
-
-  sprintf(set_name, "/Metal_free/Temperature/Temperature");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   temperature);
-  status = H5Dclose(dataset);
-
-  sprintf(set_name, "/Metal_free/Electron_density_over_n_h");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   he_electron_abundance);
-  status = H5Dclose(dataset);
-
-  for (i = 0; i < cooling->N_He; i++) {
-    for (j = 0; j < cooling->N_Temp; j++) {
-      table_index = row_major_index_2d(i, j, cooling->N_He, cooling->N_Temp);
-      cooling_index =
-          row_major_index_3d(0, i, j, 1, cooling->N_He,
-                             cooling->N_Temp);  // Redshift invariant table!!!
-      cooling_table.H_plus_He_heating[cooling_index] =
-          -he_net_cooling_rate[table_index];
-      cooling_table.H_plus_He_electron_abundance[cooling_index] =
-          he_electron_abundance[table_index];
-      cooling_table.temperature[cooling_index] =
-          log10(temperature[table_index]);
-    }
-  }
-
-  sprintf(set_name, "/Solar/Electron_density_over_n_h");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   electron_abundance);
-  status = H5Dclose(dataset);
-
-  for (i = 0; i < cooling->N_Temp; i++) {
-    cooling_table.electron_abundance[i] = electron_abundance[i];
-  }
-
-  status = H5Fclose(file_id);
-
-  free(net_cooling_rate);
-  free(electron_abundance);
-  free(temperature);
-  free(he_net_cooling_rate);
-  free(he_electron_abundance);
-
-  printf("eagle_cool_tables.h done reading in collisional table\n");
-  return cooling_table;
-}
-
-/*
- * @brief Get the cooling table for photodissociation
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-inline struct cooling_tables_redshift_invariant get_photodis_table(
-    char *cooling_table_path,
-    const struct cooling_function_data *restrict cooling) {
-
-  struct cooling_tables_redshift_invariant cooling_table;
-  hid_t file_id, dataset;
-
-  herr_t status;
-
-  char fname[500], set_name[500];
-
-  int specs, i, j, k, table_index, cooling_index;
-
-  float *net_cooling_rate;
-  float *electron_abundance;
-  float *temperature;
-  float *he_net_cooling_rate;
-  float *he_electron_abundance;
-
-  net_cooling_rate =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  electron_abundance =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                cooling->N_nH * sizeof(float));
-  he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                        cooling->N_nH * sizeof(float));
-  he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                          cooling->N_nH * sizeof(float));
-
-  cooling_table.metal_heating = (float *)malloc(
-      cooling->N_Elements * cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.electron_abundance =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                              cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_heating = (float *)malloc(
-      cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_electron_abundance = (float *)malloc(
-      cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
-
-  sprintf(fname, "%sz_photodis.hdf5", cooling_table_path);
-
-  file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-  /* For normal elements */
-  for (specs = 0; specs < cooling->N_Elements; specs++) {
-    sprintf(set_name, "/%s/Net_Cooling", cooling->ElementNames[specs]);
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     net_cooling_rate);
-    status = H5Dclose(dataset);
-
-    for (j = 0; j < cooling->N_Temp; j++) {
-      for (k = 0; k < cooling->N_nH; k++) {
-        table_index = row_major_index_2d(j, k, cooling->N_Temp, cooling->N_nH);
-        cooling_index = row_major_index_3d(
-            k, j, specs, cooling->N_nH, cooling->N_Temp,
-            cooling->N_Elements);  // Redshift invariant table!!!
-        cooling_table.metal_heating[cooling_index] =
-            -net_cooling_rate[table_index];
-      }
-    }
-  }
-
-  /* Helium */
-  sprintf(set_name, "/Metal_free/Net_Cooling");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   he_net_cooling_rate);
-  status = H5Dclose(dataset);
-
-  sprintf(set_name, "/Metal_free/Temperature/Temperature");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   temperature);
-  status = H5Dclose(dataset);
-
-  sprintf(set_name, "/Metal_free/Electron_density_over_n_h");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   he_electron_abundance);
-  status = H5Dclose(dataset);
-
-  for (i = 0; i < cooling->N_He; i++) {
-    for (j = 0; j < cooling->N_Temp; j++) {
-      for (k = 0; k < cooling->N_nH; k++) {
-        table_index = row_major_index_3d(i, j, k, cooling->N_He,
-                                         cooling->N_Temp, cooling->N_nH);
-        cooling_index =
-            row_major_index_3d(k, i, j, cooling->N_nH, cooling->N_He,
-                               cooling->N_Temp);  // Redshift invariant table!!!
-        cooling_table.H_plus_He_heating[cooling_index] =
-            -he_net_cooling_rate[table_index];
-        cooling_table.H_plus_He_electron_abundance[cooling_index] =
-            he_electron_abundance[table_index];
-        cooling_table.temperature[cooling_index] =
-            log10(temperature[table_index]);
-      }
-    }
-  }
-
-  sprintf(set_name, "/Solar/Electron_density_over_n_h");
-  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                   electron_abundance);
-  status = H5Dclose(dataset);
-
-  for (i = 0; i < cooling->N_Temp; i++) {
-    for (j = 0; j < cooling->N_nH; j++) {
-      table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
-      cooling_index = row_major_index_2d(
-          j, i, cooling->N_nH, cooling->N_Temp);  // Redshift invariant table!!!
-      cooling_table.electron_abundance[cooling_index] =
-          electron_abundance[table_index];
-    }
-  }
-
-  status = H5Fclose(file_id);
-
-  free(net_cooling_rate);
-  free(electron_abundance);
-  free(temperature);
-  free(he_net_cooling_rate);
-  free(he_electron_abundance);
-
-  printf("eagle_cool_tables.h done reading in photodissociation table\n");
-  return cooling_table;
-}
-
-inline struct cooling_tables get_two_cooling_tables(
-    char *cooling_table_path,
-    const struct cooling_function_data *restrict cooling, char *filename1,
-    char *filename2) {
-
-  struct cooling_tables cooling_table;
-  hid_t file_id, dataset;
-
-  herr_t status;
-
-  char fname[500], set_name[500];
-
-  int specs, i, j, k, table_index, cooling_index;
-
-  float *net_cooling_rate;
-  float *electron_abundance;
-  float *temperature;
-  float *he_net_cooling_rate;
-  float *he_electron_abundance;
-
-  net_cooling_rate =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  electron_abundance =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                cooling->N_nH * sizeof(float));
-  he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                        cooling->N_nH * sizeof(float));
-  he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                          cooling->N_nH * sizeof(float));
-
-  cooling_table.metal_heating =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_Elements *
-                      cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.electron_abundance = (float *)malloc(
-      cooling->N_Redshifts * cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.temperature =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_heating =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_electron_abundance =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-
-  /* For normal elements */
-  for (int file = 0; file < 2; file++) {
-    if (file == 0) {
-      sprintf(fname, "%s%s.hdf5", cooling_table_path, filename1);
-    } else {
-      sprintf(fname, "%s%s.hdf5", cooling_table_path, filename2);
-    }
-    file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-    if (file_id < 0) {
-      error("[GetCoolingTables()]: unable to open file %s\n", fname);
-    }
-
-    for (specs = 0; specs < cooling->N_Elements; specs++) {
-      sprintf(set_name, "/%s/Net_Cooling", cooling->ElementNames[specs]);
-      dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-      status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                       net_cooling_rate);
-      status = H5Dclose(dataset);
-
-      for (i = 0; i < cooling->N_nH; i++) {
-        for (j = 0; j < cooling->N_Temp; j++) {
-          table_index =
-              row_major_index_2d(j, i, cooling->N_Temp, cooling->N_nH);
-          cooling_index = row_major_index_4d(
-              file, i, j, specs, cooling->N_Redshifts, cooling->N_nH,
-              cooling->N_Temp, cooling->N_Elements);
-          cooling_table.metal_heating[cooling_index] =
-              -net_cooling_rate[table_index];
-        }
-      }
-    }
-
-    /* Helium */
-    sprintf(set_name, "/Metal_free/Net_Cooling");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     he_net_cooling_rate);
-    status = H5Dclose(dataset);
-
-    sprintf(set_name, "/Metal_free/Temperature/Temperature");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     temperature);
-    status = H5Dclose(dataset);
-
-    sprintf(set_name, "/Metal_free/Electron_density_over_n_h");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     he_electron_abundance);
-    status = H5Dclose(dataset);
-
-    for (i = 0; i < cooling->N_He; i++) {
-      for (j = 0; j < cooling->N_Temp; j++) {
-        for (k = 0; k < cooling->N_nH; k++) {
-          table_index = row_major_index_3d(i, j, k, cooling->N_He,
-                                           cooling->N_Temp, cooling->N_nH);
-          cooling_index =
-              row_major_index_4d(file, k, i, j, cooling->N_Redshifts,
-                                 cooling->N_nH, cooling->N_He, cooling->N_Temp);
-          cooling_table.H_plus_He_heating[cooling_index] =
-              -he_net_cooling_rate[table_index];
-          cooling_table.H_plus_He_electron_abundance[cooling_index] =
-              he_electron_abundance[table_index];
-          cooling_table.temperature[cooling_index] =
-              log10(temperature[table_index]);
-        }
-      }
-    }
-
-    sprintf(set_name, "/Solar/Electron_density_over_n_h");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     electron_abundance);
-    status = H5Dclose(dataset);
-
-    for (i = 0; i < cooling->N_Temp; i++) {
-      for (j = 0; j < cooling->N_nH; j++) {
-        table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
-        cooling_index = row_major_index_3d(file, j, i, cooling->N_Redshifts,
-                                           cooling->N_nH, cooling->N_Temp);
-        cooling_table.electron_abundance[cooling_index] =
-            electron_abundance[table_index];
-      }
-    }
-
-    status = H5Fclose(file_id);
-  }
-
-  free(net_cooling_rate);
-  free(electron_abundance);
-  free(temperature);
-  free(he_net_cooling_rate);
-  free(he_electron_abundance);
-
-  printf("eagle_cool_tables.h done reading in general cooling table\n");
-
-  return cooling_table;
-}
-
-/*
- * @brief Get the cooling tables dependent on redshift
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-inline struct cooling_tables get_cooling_table(
-    char *cooling_table_path,
-    const struct cooling_function_data *restrict cooling) {
-
-  struct cooling_tables cooling_table;
-  hid_t file_id, dataset;
-
-  herr_t status;
-
-  char fname[500], set_name[500];
-
-  int specs, i, j, k, table_index, cooling_index;
-
-  float *net_cooling_rate;
-  float *electron_abundance;
-  float *temperature;
-  float *he_net_cooling_rate;
-  float *he_electron_abundance;
-
-  net_cooling_rate =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  electron_abundance =
-      (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
-  temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                cooling->N_nH * sizeof(float));
-  he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                        cooling->N_nH * sizeof(float));
-  he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
-                                          cooling->N_nH * sizeof(float));
-
-  cooling_table.metal_heating =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_Elements *
-                      cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.electron_abundance = (float *)malloc(
-      cooling->N_Redshifts * cooling->N_Temp * cooling->N_nH * sizeof(float));
-  cooling_table.temperature =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_heating =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-  cooling_table.H_plus_He_electron_abundance =
-      (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
-                      cooling->N_nH * sizeof(float));
-
-  /* For normal elements */
-  for (int z_index = 0; z_index < cooling->N_Redshifts; z_index++) {
-    sprintf(fname, "%sz_%1.3f.hdf5", cooling_table_path,
-            cooling->Redshifts[z_index]);
-    file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
-    if (file_id < 0) {
-      error("[GetCoolingTables()]: unable to open file %s\n", fname);
-    }
-
-    for (specs = 0; specs < cooling->N_Elements; specs++) {
-      sprintf(set_name, "/%s/Net_Cooling", cooling->ElementNames[specs]);
-      dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-      status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                       net_cooling_rate);
-      status = H5Dclose(dataset);
-
-      for (i = 0; i < cooling->N_nH; i++) {
-        for (j = 0; j < cooling->N_Temp; j++) {
-          table_index =
-              row_major_index_2d(j, i, cooling->N_Temp, cooling->N_nH);
-          cooling_index = row_major_index_4d(
-              z_index, i, j, specs, cooling->N_Redshifts, cooling->N_nH,
-              cooling->N_Temp, cooling->N_Elements);
-          cooling_table.metal_heating[cooling_index] =
-              -net_cooling_rate[table_index];
-        }
-      }
-    }
-
-    /* Helium */
-    sprintf(set_name, "/Metal_free/Net_Cooling");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     he_net_cooling_rate);
-    status = H5Dclose(dataset);
-
-    sprintf(set_name, "/Metal_free/Temperature/Temperature");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     temperature);
-    status = H5Dclose(dataset);
-
-    sprintf(set_name, "/Metal_free/Electron_density_over_n_h");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     he_electron_abundance);
-    status = H5Dclose(dataset);
-
-    for (i = 0; i < cooling->N_He; i++) {
-      for (j = 0; j < cooling->N_Temp; j++) {
-        for (k = 0; k < cooling->N_nH; k++) {
-          table_index = row_major_index_3d(i, j, k, cooling->N_He,
-                                           cooling->N_Temp, cooling->N_nH);
-          cooling_index =
-              row_major_index_4d(z_index, k, i, j, cooling->N_Redshifts,
-                                 cooling->N_nH, cooling->N_He, cooling->N_Temp);
-          cooling_table.H_plus_He_heating[cooling_index] =
-              -he_net_cooling_rate[table_index];
-          cooling_table.H_plus_He_electron_abundance[cooling_index] =
-              he_electron_abundance[table_index];
-          cooling_table.temperature[cooling_index] =
-              log10(temperature[table_index]);
-        }
-      }
-    }
-
-    sprintf(set_name, "/Solar/Electron_density_over_n_h");
-    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
-    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                     electron_abundance);
-    status = H5Dclose(dataset);
-
-    for (i = 0; i < cooling->N_Temp; i++) {
-      for (j = 0; j < cooling->N_nH; j++) {
-        table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
-        cooling_index = row_major_index_3d(z_index, j, i, cooling->N_Redshifts,
-                                           cooling->N_nH, cooling->N_Temp);
-        cooling_table.electron_abundance[cooling_index] =
-            electron_abundance[table_index];
-      }
-    }
-
-    status = H5Fclose(file_id);
-  }
-
-  free(net_cooling_rate);
-  free(electron_abundance);
-  free(temperature);
-  free(he_net_cooling_rate);
-  free(he_electron_abundance);
-
-  printf("eagle_cool_tables.h done reading in general cooling table\n");
-
-  return cooling_table;
-}
-
-/*
- * @brief Constructs the data structure containting all the cooling tables
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-inline void eagle_read_two_tables(
-    struct cooling_function_data *restrict cooling) {
-
-  struct eagle_cooling_table table;
-
-  table.element_cooling =
-      get_cooling_table(cooling->cooling_table_path, cooling);
-  cooling->table = table;
-}
-
-/*
- * @brief Constructs the data structure containting all the cooling tables
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-inline struct eagle_cooling_table eagle_readtable(
-    char *cooling_table_path,
-    const struct cooling_function_data *restrict cooling) {
-
-  struct eagle_cooling_table table;
-
-  table.no_compton_cooling = get_no_compt_table(cooling_table_path, cooling);
-  table.photodissociation_cooling =
-      get_photodis_table(cooling_table_path, cooling);
-  table.collisional_cooling =
-      get_collisional_table(cooling_table_path, cooling);
-  table.element_cooling = get_cooling_table(cooling_table_path, cooling);
-
-  return table;
-}
-
-/*
- * @brief Finds the element index for the corresponding element in swift
- *
- * @param element_name Element string we want to match to element index
- * @param cooling Cooling data structure
- */
-inline int element_index(char *element_name,
-                         const struct cooling_function_data *restrict cooling) {
-  int i;
-
-  for (i = 0; i < cooling->N_Elements; i++)
-    if (strcmp(cooling->ElementNames[i], element_name) == 0) return i;
-
-  /* element not found */
-  return -1;
-}
-
-/*
- * @brief Finds the element index for the corresponding element in EAGLE
- *
- * @param element_name Element string we want to match to element index
- * @param size Number of elements tracked in EAGLE
- * @param cooling Cooling data structure
- */
-inline int get_element_index(char *table[20], int size, char *element_name) {
-  int i;
-
-  for (i = 0; i < size; i++)
-    if (strcmp(table[i], element_name) == 0) return i;
-
-  /* element not found */
-  return -1;
-}
-
-/*
- * @brief Makes an array of element names which are tracked in the EAGLE tables
- *
- * @param cooling Cooling data structure
- */
-inline void MakeNamePointers(struct cooling_function_data *cooling) {
-  int i, j, sili_index = 0;
-  char ElementNames[cooling->N_Elements][eagle_element_name_length];
+void GetCoolingRedshifts(struct cooling_function_data *);
 
-  /* This is ridiculous, way too many element name arrays. Needs to be changed
-   */
-  // ElementNames =
-  // malloc(cooling->N_Elements*eagle_element_name_length*sizeof(char));
-  strcpy(ElementNames[0], "Hydrogen");
-  strcpy(ElementNames[1], "Helium");
-  strcpy(ElementNames[2], "Carbon");
-  strcpy(ElementNames[3], "Nitrogen");
-  strcpy(ElementNames[4], "Oxygen");
-  strcpy(ElementNames[5], "Neon");
-  strcpy(ElementNames[6], "Magnesium");
-  strcpy(ElementNames[7], "Silicon");
-  strcpy(ElementNames[8], "Iron");
+void ReadCoolingHeader(char *, struct cooling_function_data *);
 
-  cooling->ElementNamePointers = malloc(cooling->N_Elements * sizeof(int));
-  cooling->SolarAbundanceNamePointers =
-      malloc(cooling->N_Elements * sizeof(int));
+struct cooling_tables_redshift_invariant get_no_compt_table(
+  char *, const struct cooling_function_data *restrict);
 
-  for (i = 0; i < cooling->N_Elements; i++) {
-    if (strcmp(ElementNames[i], "Silicon") == 0) sili_index = i;
-  }
+struct cooling_tables_redshift_invariant get_collisional_table(
+    char *, const struct cooling_function_data *restrict);
 
-  for (i = 0; i < cooling->N_Elements; i++) {
-    cooling->SolarAbundanceNamePointers[i] = -999;
-    cooling->ElementNamePointers[i] = -999;
+struct cooling_tables_redshift_invariant get_photodis_table(
+    char *, const struct cooling_function_data *restrict);
 
-    for (j = 0; j < cooling->N_SolarAbundances; j++) {
-      if (strcmp(cooling->ElementNames[i], cooling->SolarAbundanceNames[j]) ==
-          0)
-        cooling->SolarAbundanceNamePointers[i] = j;
-    }
+struct cooling_tables get_cooling_table(
+    char *, const struct cooling_function_data *restrict);
 
-    if (strcmp(cooling->ElementNames[i], "Sulphur") == 0 ||
-        strcmp(cooling->ElementNames[i], "Calcium") ==
-            0) /* These elements are tracked! */
-      cooling->ElementNamePointers[i] = -1 * sili_index;
-    else {
-      for (j = 0; j < cooling->N_Elements; j++) {
-        if (strcmp(cooling->ElementNames[i], ElementNames[j]) == 0)
-          cooling->ElementNamePointers[i] = j;
-      }
-    }
-  }
-}
+struct eagle_cooling_table eagle_readtable(
+    char *, const struct cooling_function_data *restrict);
 
 #endif
diff --git a/src/cooling/EAGLE/interpolate.h b/src/cooling/EAGLE/interpolate.h
index 151bc735067d55481befa5ee5f62815b93174c27..ed82ee26e5a02e9663b9dfe1524b506a1871eb75 100644
--- a/src/cooling/EAGLE/interpolate.h
+++ b/src/cooling/EAGLE/interpolate.h
@@ -21,7 +21,7 @@
 
 /**
  * @file src/cooling/EAGLE/interpolate.h
- * @brief EAGLE interpolation function declarations
+ * @brief Interpolation functions for EAGLE tables
  */
 
 /* Config parameters. */
@@ -36,6 +36,7 @@
 /* Local includes. */
 #include "chemistry.h"
 #include "cooling_struct.h"
+#include "eagle_cool_tables.h"
 #include "error.h"
 #include "hydro.h"
 #include "io_properties.h"
@@ -44,90 +45,415 @@
 #include "physical_constants.h"
 #include "units.h"
 
-int row_major_index_2d(int, int, int, int);
+static int get_redshift_index_first_call = 0;
+static int get_redshift_index_previous = -1;
 
-int row_major_index_3d(int, int, int, int, int, int);
+static const float EPS = 1.e-4;
 
-int row_major_index_4d(int, int, int, int, int, int, int, int);
+/**
+ * @brief Returns the 1d index of element with 2d indices i,j
+ * from a flattened 2d array in row major order
+ *
+ * @param i, j Indices of element of interest
+ * @param nx, ny Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_2d(int i, int j,
+                                                             int nx, int ny) {
+  return i * ny + j;
+}
 
-void get_index_1d(float *, int, double, int *, float *);
+/**
+ * @brief Returns the 1d index of element with 3d indices i,j,k
+ * from a flattened 3d array in row major order
+ *
+ * @param i, j, k Indices of element of interest
+ * @param nx, ny, nz Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_3d(int i, int j,
+                                                             int k, int nx,
+                                                             int ny, int nz) {
+  return i * ny * nz + j * nz + k;
+}
 
-void get_redshift_index(float, int *, float *,
-                        const struct cooling_function_data *restrict);
+/**
+ * @brief Returns the 1d index of element with 4d indices i,j,k,l
+ * from a flattened 4d array in row major order
+ *
+ * @param i, j, k, l Indices of element of interest
+ * @param nx, ny, nz, nw Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_4d(int i, int j,
+                                                             int k, int l,
+                                                             int nx, int ny,
+                                                             int nz, int nw) {
+  return i * ny * nz * nw + j * nz * nw + k * nw + l;
+}
 
-float interpol_1d(float *, int, float);
+/*
+ * @brief This routine returns the position i of a value x in a 1D table and the
+ * displacement dx needed for the interpolation.  The table is assumed to
+ * be evenly spaced.
+ *
+ * @param table Pointer to table of values
+ * @param ntable Size of the table
+ * @param x Value whose position within the array we are interested in
+ * @param i Pointer to the index whose corresponding table value
+ * is the greatest value in the table less than x
+ * @param dx Pointer to offset of x within table cell
+ */
+__attribute__((always_inline)) INLINE void get_index_1d(float *table,
+                                                        int ntable, double x,
+                                                        int *i, float *dx) {
 
-double interpol_1d_dbl(double *, int, float);
+  float dxm1 = (float)(ntable - 1) / (table[ntable - 1] - table[0]);
 
-float interpol_2d(float *, int, int, float, float, int, int, double *,
-                  double *);
+  if ((float)x <= table[0] + EPS) {
+    *i = 0;
+    *dx = 0;
+  } else if ((float)x >= table[ntable - 1] - EPS) {
+    *i = ntable - 2;
+    *dx = 1;
+  } else {
+    *i = (int)floor(((float)x - table[0]) * dxm1);
+    *dx = ((float)x - table[*i]) * dxm1;
+  }
+}
 
-double interpol_2d_dbl(double *, int, int, double, double, int, int, double *,
-                       double *);
+/*
+ * @brief Returns the position i of a value z in the redshift table
+ * and computes the displacement dz needed for the interpolation.
+ * Since the redshift table is not evenly spaced, compare z with each
+ * table value in decreasing order starting with the previous redshift index
+ *
+ * @param z Redshift whose position within the redshift array we are interested
+ * in
+ * @param z_index Pointer to the index whose corresponding redshift
+ * is the greatest value in the redshift table less than x
+ * @param dz Pointer to offset of z within redshift cell
+ * @param cooling Pointer to cooling structure containing redshift table
+ */
+__attribute__((always_inline)) INLINE void get_redshift_index(
+    float z, int *z_index, float *dz,
+    const struct cooling_function_data *restrict cooling) {
+  int i, iz;
 
-float interpol_3d(float *, int, int, int, float, float, float, int, int, int,
-                  double *, double *);
+  if (get_redshift_index_first_call == 0) {
+    get_redshift_index_first_call = 1;
+    get_redshift_index_previous = cooling->N_Redshifts - 2;
 
-float interpol_4d(float *, int, int, int, int, float, float, float, float, int,
-                  int, int, int, double *, double *);
+    /* this routine assumes cooling_redshifts table is in increasing order. Test
+     * this. */
+    for (i = 0; i < cooling->N_Redshifts - 2; i++)
+      if (cooling->Redshifts[i + 1] < cooling->Redshifts[i]) {
+        error("[get_redshift_index]: table should be in increasing order\n");
+      }
+  }
 
-void construct_1d_table_from_2d(const struct part *restrict,
-                                const struct cooling_function_data *restrict,
-                                const struct cosmology *restrict,
-                                const struct phys_const *, float *, int, float,
-                                int, int, double *, float *, float *);
+  /* before the earliest redshift or before hydrogen reionization, flag for
+   * collisional cooling */
+  if (z > cooling->reionisation_redshift) {
+    *z_index = cooling->N_Redshifts;
+    *dz = 0.0;
+  }
+  /* from reionization use the cooling tables */
+  else if (z > cooling->Redshifts[cooling->N_Redshifts - 1] &&
+           z <= cooling->reionisation_redshift) {
+    *z_index = cooling->N_Redshifts + 1;
+    *dz = 0.0;
+  }
+  /* at the end, just use the last value */
+  else if (z <= cooling->Redshifts[0]) {
+    *z_index = 0;
+    *dz = 0.0;
+  } else {
+    /* start at the previous index and search */
+    for (iz = get_redshift_index_previous; iz >= 0; iz--) {
+      if (z > cooling->Redshifts[iz]) {
+        *dz = (z - cooling->Redshifts[iz]) /
+              (cooling->Redshifts[iz + 1] - cooling->Redshifts[iz]);
 
-void construct_1d_table_from_3d(const struct part *restrict,
-                                const struct cooling_function_data *restrict,
-                                const struct cosmology *restrict,
-                                const struct phys_const *, float *, int, float,
-                                int, int, float, int, int, double *, float *,
-                                float *);
+        get_redshift_index_previous = *z_index = iz;
 
-void construct_1d_print_table_from_3d_elements(
-    const struct part *restrict, const struct cooling_function_data *restrict,
-    const struct cosmology *restrict, const struct phys_const *, float *, int,
-    float, int, int, double *, float *, float *, float *);
+        break;
+      }
+    }
+  }
+}
 
-void construct_1d_table_from_3d_elements(
-    const struct part *restrict, const struct cooling_function_data *restrict,
-    const struct cosmology *restrict, const struct phys_const *, float *, int,
-    float, int, int, double *, float *, float *, float *);
+/*
+ * @brief Performs 1d interpolation (double precision)
+ *
+ * @param table Pointer to 1d table of values
+ * @param i Index of cell we are interpolating
+ * @param dx Offset within cell
+ */
+__attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table,
+                                                             int i, float dx) {
+
+  return (1 - dx) * table[i] + dx * table[i + 1];
+}
+
+/*
+ * @brief Performs 2d interpolation
+ *
+ * @param table Pointer to flattened 2d table of values
+ * @param i,j Indices of cell we are interpolating
+ * @param dx,dy Offset within cell
+ * @param nx,ny Table dimensions
+ * @param upper Pointer to value set to the table value at
+ * the when dy = 1 (used for calculating derivatives)
+ * @param lower Pointer to value set to the table value at
+ * the when dy = 0 (used for calculating derivatives)
+ */
+__attribute__((always_inline)) INLINE float interpol_2d(float *table, int i,
+                                                        int j, float dx,
+                                                        float dy, int nx,
+                                                        int ny, double *upper,
+                                                        double *lower) {
+  float result;
+  int index[4];
 
-void construct_1d_table_from_4d(const struct part *restrict,
-                                const struct cooling_function_data *restrict,
-                                const struct cosmology *restrict,
-                                const struct phys_const *, float *, int, float,
-                                int, int, float, int, int, float, int, int,
-                                double *, float *, float *);
+  index[0] = row_major_index_2d(i, j, nx, ny);
+  index[1] = row_major_index_2d(i, j + 1, nx, ny);
+  index[2] = row_major_index_2d(i + 1, j, nx, ny);
+  index[3] = row_major_index_2d(i + 1, j + 1, nx, ny);
 
-void construct_1d_print_table_from_4d_elements(
-    const struct part *restrict, const struct cooling_function_data *restrict,
-    const struct cosmology *restrict, const struct phys_const *, float *, int,
-    float, int, int, float, int, int, double *, float *, float *, float *);
+  result = (1 - dx) * (1 - dy) * table[index[0]] +
+           (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
+           dx * dy * table[index[3]];
+  dy = 1.0;
+  *upper = (1 - dx) * (1 - dy) * table[index[0]] +
+           (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
+           dx * dy * table[index[3]];
+  dy = 0.0;
+  *lower = (1 - dx) * (1 - dy) * table[index[0]] +
+           (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
+           dx * dy * table[index[3]];
+
+  return result;
+}
+
+/*
+ * @brief Performs 3d interpolation
+ *
+ * @param table Pointer to flattened 3d table of values
+ * @param i,j,k Indices of cell we are interpolating
+ * @param dx,dy,dz Offset within cell
+ * @param nx,ny,nz Table dimensions
+ * @param upper Pointer to value set to the table value at
+ * the when dz = 1 (used for calculating derivatives)
+ * @param lower Pointer to value set to the table value at
+ * the when dz = 0 (used for calculating derivatives)
+ */
+__attribute__((always_inline)) INLINE float interpol_3d(
+    float *table, int i, int j, int k, float dx, float dy, float dz, int nx,
+    int ny, int nz, double *upper, double *lower) {
+  float result;
+  int index[8];
+
+  index[0] = row_major_index_3d(i, j, k, nx, ny, nz);
+  index[1] = row_major_index_3d(i, j, k + 1, nx, ny, nz);
+  index[2] = row_major_index_3d(i, j + 1, k, nx, ny, nz);
+  index[3] = row_major_index_3d(i, j + 1, k + 1, nx, ny, nz);
+  index[4] = row_major_index_3d(i + 1, j, k, nx, ny, nz);
+  index[5] = row_major_index_3d(i + 1, j, k + 1, nx, ny, nz);
+  index[6] = row_major_index_3d(i + 1, j + 1, k, nx, ny, nz);
+  index[7] = row_major_index_3d(i + 1, j + 1, k + 1, nx, ny, nz);
+
+  float table0 = table[index[0]];
+  float table1 = table[index[1]];
+  float table2 = table[index[2]];
+  float table3 = table[index[3]];
+  float table4 = table[index[4]];
+  float table5 = table[index[5]];
+  float table6 = table[index[6]];
+  float table7 = table[index[7]];
+
+  result = (1 - dx) * (1 - dy) * (1 - dz) * table0 +
+           (1 - dx) * (1 - dy) * dz * table1 +
+           (1 - dx) * dy * (1 - dz) * table2 + (1 - dx) * dy * dz * table3 +
+           dx * (1 - dy) * (1 - dz) * table4 + dx * (1 - dy) * dz * table5 +
+           dx * dy * (1 - dz) * table6 + dx * dy * dz * table7;
+  dz = 1.0;
+  *upper = (1 - dx) * (1 - dy) * (1 - dz) * table0 +
+           (1 - dx) * (1 - dy) * dz * table1 +
+           (1 - dx) * dy * (1 - dz) * table2 + (1 - dx) * dy * dz * table3 +
+           dx * (1 - dy) * (1 - dz) * table4 + dx * (1 - dy) * dz * table5 +
+           dx * dy * (1 - dz) * table6 + dx * dy * dz * table7;
+  dz = 0.0;
+  *lower = (1 - dx) * (1 - dy) * (1 - dz) * table0 +
+           (1 - dx) * (1 - dy) * dz * table1 +
+           (1 - dx) * dy * (1 - dz) * table2 + (1 - dx) * dy * dz * table3 +
+           dx * (1 - dy) * (1 - dz) * table4 + dx * (1 - dy) * dz * table5 +
+           dx * dy * (1 - dz) * table6 + dx * dy * dz * table7;
+
+  return result;
+}
+
+/*
+ * @brief Performs 4d interpolation
+ *
+ * @param table Pointer to flattened 4d table of values
+ * @param i,j,k,l Indices of cell we are interpolating
+ * @param dx,dy,dz,dw Offset within cell
+ * @param nx,ny,nz,nw Table dimensions
+ * @param upper Pointer to value set to the table value at
+ * the when dw = 1 when used for interpolating table
+ * depending on metal species, dz = 1 otherwise
+ * (used for calculating derivatives)
+ * @param upper Pointer to value set to the table value at
+ * the when dw = 0 when used for interpolating table
+ * depending on metal species, dz = 0 otherwise
+ * (used for calculating derivatives)
+ */
+__attribute__((always_inline)) INLINE float interpol_4d(
+    float *table, int i, int j, int k, int l, float dx, float dy, float dz,
+    float dw, int nx, int ny, int nz, int nw, double *upper, double *lower) {
+  float result;
+  int index[16];
+
+  index[0] = row_major_index_4d(i, j, k, l, nx, ny, nz, nw);
+  index[1] = row_major_index_4d(i, j, k, l + 1, nx, ny, nz, nw);
+  index[2] = row_major_index_4d(i, j, k + 1, l, nx, ny, nz, nw);
+  index[3] = row_major_index_4d(i, j, k + 1, l + 1, nx, ny, nz, nw);
+  index[4] = row_major_index_4d(i, j + 1, k, l, nx, ny, nz, nw);
+  index[5] = row_major_index_4d(i, j + 1, k, l + 1, nx, ny, nz, nw);
+  index[6] = row_major_index_4d(i, j + 1, k + 1, l, nx, ny, nz, nw);
+  index[7] = row_major_index_4d(i, j + 1, k + 1, l + 1, nx, ny, nz, nw);
+  index[8] = row_major_index_4d(i + 1, j, k, l, nx, ny, nz, nw);
+  index[9] = row_major_index_4d(i + 1, j, k, l + 1, nx, ny, nz, nw);
+  index[10] = row_major_index_4d(i + 1, j, k + 1, l, nx, ny, nz, nw);
+  index[11] = row_major_index_4d(i + 1, j, k + 1, l + 1, nx, ny, nz, nw);
+  index[12] = row_major_index_4d(i + 1, j + 1, k, l, nx, ny, nz, nw);
+  index[13] = row_major_index_4d(i + 1, j + 1, k, l + 1, nx, ny, nz, nw);
+  index[14] = row_major_index_4d(i + 1, j + 1, k + 1, l, nx, ny, nz, nw);
+  index[15] = row_major_index_4d(i + 1, j + 1, k + 1, l + 1, nx, ny, nz, nw);
+
+  float table0 = table[index[0]];
+  float table1 = table[index[1]];
+  float table2 = table[index[2]];
+  float table3 = table[index[3]];
+  float table4 = table[index[4]];
+  float table5 = table[index[5]];
+  float table6 = table[index[6]];
+  float table7 = table[index[7]];
+  float table8 = table[index[8]];
+  float table9 = table[index[9]];
+  float table10 = table[index[10]];
+  float table11 = table[index[11]];
+  float table12 = table[index[12]];
+  float table13 = table[index[13]];
+  float table14 = table[index[14]];
+  float table15 = table[index[15]];
+
+  result = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table0 +
+           (1 - dx) * (1 - dy) * (1 - dz) * dw * table1 +
+           (1 - dx) * (1 - dy) * dz * (1 - dw) * table2 +
+           (1 - dx) * (1 - dy) * dz * dw * table3 +
+           (1 - dx) * dy * (1 - dz) * (1 - dw) * table4 +
+           (1 - dx) * dy * (1 - dz) * dw * table5 +
+           (1 - dx) * dy * dz * (1 - dw) * table6 +
+           (1 - dx) * dy * dz * dw * table7 +
+           dx * (1 - dy) * (1 - dz) * (1 - dw) * table8 +
+           dx * (1 - dy) * (1 - dz) * dw * table9 +
+           dx * (1 - dy) * dz * (1 - dw) * table10 +
+           dx * (1 - dy) * dz * dw * table11 +
+           dx * dy * (1 - dz) * (1 - dw) * table12 +
+           dx * dy * (1 - dz) * dw * table13 +
+           dx * dy * dz * (1 - dw) * table14 + dx * dy * dz * dw * table15;
+  if (nw == 9) {
+    // interpolating metal species
+    dz = 1.0;
+  } else {
+    dw = 1.0;
+  }
+  *upper = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table0 +
+           (1 - dx) * (1 - dy) * (1 - dz) * dw * table1 +
+           (1 - dx) * (1 - dy) * dz * (1 - dw) * table2 +
+           (1 - dx) * (1 - dy) * dz * dw * table3 +
+           (1 - dx) * dy * (1 - dz) * (1 - dw) * table4 +
+           (1 - dx) * dy * (1 - dz) * dw * table5 +
+           (1 - dx) * dy * dz * (1 - dw) * table6 +
+           (1 - dx) * dy * dz * dw * table7 +
+           dx * (1 - dy) * (1 - dz) * (1 - dw) * table8 +
+           dx * (1 - dy) * (1 - dz) * dw * table9 +
+           dx * (1 - dy) * dz * (1 - dw) * table10 +
+           dx * (1 - dy) * dz * dw * table11 +
+           dx * dy * (1 - dz) * (1 - dw) * table12 +
+           dx * dy * (1 - dz) * dw * table13 +
+           dx * dy * dz * (1 - dw) * table14 + dx * dy * dz * dw * table15;
+  if (nw == 9) {
+    // interpolating metal species
+    dz = 0.0;
+  } else {
+    dw = 0.0;
+  }
+  *lower = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table0 +
+           (1 - dx) * (1 - dy) * (1 - dz) * dw * table1 +
+           (1 - dx) * (1 - dy) * dz * (1 - dw) * table2 +
+           (1 - dx) * (1 - dy) * dz * dw * table3 +
+           (1 - dx) * dy * (1 - dz) * (1 - dw) * table4 +
+           (1 - dx) * dy * (1 - dz) * dw * table5 +
+           (1 - dx) * dy * dz * (1 - dw) * table6 +
+           (1 - dx) * dy * dz * dw * table7 +
+           dx * (1 - dy) * (1 - dz) * (1 - dw) * table8 +
+           dx * (1 - dy) * (1 - dz) * dw * table9 +
+           dx * (1 - dy) * dz * (1 - dw) * table10 +
+           dx * (1 - dy) * dz * dw * table11 +
+           dx * dy * (1 - dz) * (1 - dw) * table12 +
+           dx * dy * (1 - dz) * dw * table13 +
+           dx * dy * dz * (1 - dw) * table14 + dx * dy * dz * dw * table15;
+
+  return result;
+}
+
+/*
+ * @brief Interpolates temperature from internal energy based on table and
+ * calculates the size of the internal energy cell for the specified
+ * internal energy.
+ *
+ * @param log_10_u Log base 10 of internal energy
+ * @param delta_u Pointer to size of internal energy cell
+ * @param z_i Redshift index
+ * @param n_h_i Hydrogen number density index
+ * @param He_i Helium fraction index
+ * @param d_z Redshift offset
+ * @param d_n_h Hydrogen number density offset
+ * @param d_He Helium fraction offset
+ * @param cooling Cooling data structure
+ * @param cosmo Cosmology data structure
+ */
+__attribute__((always_inline)) INLINE double eagle_convert_u_to_temp(
+    double log_10_u, float *delta_u, int z_i, int n_h_i, int He_i, float d_z,
+    float d_n_h, float d_He,
+    const struct cooling_function_data *restrict cooling,
+    const struct cosmology *restrict cosmo) {
 
-void construct_1d_table_from_4d_elements(
-    const struct part *restrict, const struct cooling_function_data *restrict,
-    const struct cosmology *restrict, const struct phys_const *, float *, int,
-    float, int, int, float, int, int, double *, float *, float *, float *);
+  int u_i;
+  float d_u, logT;
+  double upper, lower;
 
-double eagle_convert_temp_to_u_1d_table(
-    double, float *, const struct cooling_function_data *restrict);
+  get_index_1d(cooling->Therm, cooling->N_Temp, log_10_u, &u_i, &d_u);
 
-double eagle_convert_u_to_temp(double, float *, int, int, int, float, float,
-                               float,
-                               const struct cooling_function_data *restrict,
-                               const struct cosmology *restrict);
+  if (cosmo->z > cooling->reionisation_redshift) {
+    logT = interpol_3d(cooling->table.photodissociation_cooling.temperature,
+                       n_h_i, He_i, u_i, d_n_h, d_He, d_u, cooling->N_nH,
+                       cooling->N_He, cooling->N_Temp, &upper, &lower);
+  } else if (cosmo->z > cooling->Redshifts[cooling->N_Redshifts - 1]) {
+    logT = interpol_3d(cooling->table.no_compton_cooling.temperature, n_h_i,
+                       He_i, u_i, d_n_h, d_He, d_u, cooling->N_nH,
+                       cooling->N_He, cooling->N_Temp, &upper, &lower);
+  } else {
+    logT = interpol_4d(cooling->table.element_cooling.temperature, z_i, n_h_i,
+                       He_i, u_i, d_z, d_n_h, d_He, d_u, cooling->N_Redshifts,
+                       cooling->N_nH, cooling->N_He, cooling->N_Temp, &upper,
+                       &lower);
+  }
 
-double eagle_convert_u_to_temp_1d_table(
-    double, float *, double *, const struct cooling_function_data *restrict);
+  *delta_u =
+      pow(10.0, cooling->Therm[u_i + 1]) - pow(10.0, cooling->Therm[u_i]);
 
-void construct_1d_tables(int, float, int, float, int, float,
-                         const struct phys_const *restrict,
-                         const struct cosmology *restrict,
-                         const struct cooling_function_data *restrict,
-                         const struct part *restrict, float *, double *,
-                         double *, double *, double *, double *, double *,
-                         float *, float *);
+  return logT;
+}
 
-#endif /* SWIFT_INTERPOL_EAGLE_H */
+#endif
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b816be7e9fb81057079cc4acac9d41d4e4a1a242..640c837b3bd1c3bf19d3d57ef080839b170427bf 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -20,7 +20,7 @@ AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) $(GSL_INCS) $(FFTW_INCS)
 AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS) $(PROFILER_LIBS)
 
 # List of programs and scripts to run in the test suite
-TESTS = testGreetings testCooling testMaths testReading.sh testSingle testKernel testSymmetry \
+TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
         testActivePair.sh test27cells.sh test27cellsPerturbed.sh  \
         testParser.sh testSPHStep test125cells.sh test125cellsPerturbed.sh testFFT \
         testAdiabaticIndex \
@@ -32,7 +32,7 @@ TESTS = testGreetings testCooling testMaths testReading.sh testSingle testKernel
 	test27cellsStars.sh test27cellsStarsPerturbed.sh
 
 # List of test programs to compile
-check_PROGRAMS = testGreetings testCooling testReading testSingle testTimeIntegration \
+check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
 		 testSPHStep testActivePair test27cells test27cells_subset test125cells testParser \
                  testKernel testFFT testInteractions testMaths \
                  testSymmetry testThreadpool \
@@ -49,8 +49,6 @@ $(check_PROGRAMS): ../src/.libs/libswiftsim.a
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
 
-testCooling_SOURCES = testCooling.c
-
 testMaths_SOURCES = testMaths.c
 
 testReading_SOURCES = testReading.c
diff --git a/tests/plots.py b/tests/plots.py
deleted file mode 100644
index 56d90d1dc79fed9c6ca9e24b64697c45a30f9baa..0000000000000000000000000000000000000000
--- a/tests/plots.py
+++ /dev/null
@@ -1,40 +0,0 @@
-import matplotlib.pyplot as plt
-
-elements = 11
-temperature = []
-cooling_rate = [[] for i in range(elements+1)]
-length = 0
-file_in = open('cooling_output.dat', 'r')
-for line in file_in:
-	data = line.split()
-	temperature.append(float(data[0]))
-	cooling_rate[0].append(-float(data[1]))
-
-file_in.close()
-
-for elem in range(elements):
-	file_in = open('cooling_output_'+str(elem)+'.dat','r')
-	for line in file_in:
-        	data = line.split()
-        	cooling_rate[elem+1].append(-float(data[0]))
-	file_in.close()
-
-
-p0, = plt.loglog(temperature, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
-p2, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
-p3, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
-p4, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
-p5, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
-p6, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
-p7, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
-p8, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
-p9, = plt.loglog(temperature, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
-p10, = plt.loglog(temperature, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
-plt.ylim([1e-24,1e-21])
-plt.xlim([1e4,1e8])
-plt.xlabel('Temperature (K)')
-plt.ylabel('Cooling rate (eagle units...)')
-plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9])
-#plt.legend(handles = [p0,p2,p3,p4,p5,p6,p7,p8,p9])
-plt.show()
diff --git a/tests/testCooling.c b/tests/testCooling.c
deleted file mode 100644
index 40eb7cd4596da1fdec165b244745af3c6728e01a..0000000000000000000000000000000000000000
--- a/tests/testCooling.c
+++ /dev/null
@@ -1,232 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-int main() {
-
-  /* Read the parameter file */
-  struct swift_params *params = malloc(sizeof(struct swift_params));
-  struct unit_system us;
-  struct chemistry_data chemistry_data;
-  struct part p;
-  struct xpart xp;
-  struct phys_const internal_const;
-  struct cooling_function_data cooling;
-  struct cosmology cosmo;
-  char *parametersFileName = "../examples/CoolingBox/coolingBox.yml";
-  enum table_index {
-    EAGLE_Hydrogen = 0,
-    EAGLE_Helium,
-    EAGLE_Carbon,
-    EAGLE_Nitrogen,
-    EAGLE_Oxygen,
-    EAGLE_Neon,
-    EAGLE_Magnesium,
-    EAGLE_Silicon,
-    EAGLE_Iron
-  };
-
-  if (params == NULL) error("Error allocating memory for the parameter file.");
-  message("Reading runtime parameters from file '%s'", parametersFileName);
-  parser_read_file(parametersFileName, params);
-
-  /* And dump the parameters as used. */
-  parser_write_params_to_file(params, "used_parameters.yml");
-
-  units_init(&us, params, "InternalUnitSystem");
-  phys_const_init(&us, params, &internal_const);
-
-  double hydrogen_number_density_cgs = 1e-1;
-  double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
-      temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
-  u_cgs = 0;
-  cooling_du_dt = 0;
-  temperature_cgs = 0;
-  newton_func = 0;
-  // int n_t_i = 2000;
-  dt_cgs = 4.0e-8 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
-
-  char *output_filename = "cooling_output.dat";
-  FILE *output_file = fopen(output_filename, "w");
-  if (output_file == NULL) {
-    printf("Error opening file!\n");
-    exit(1);
-  }
-  char *output_filename2 = "temperature_output.dat";
-  FILE *output_file2 = fopen(output_filename2, "w");
-  if (output_file2 == NULL) {
-    printf("Error opening file!\n");
-    exit(1);
-  }
-  char *output_filename3 = "newton_output.dat";
-  FILE *output_file3 = fopen(output_filename3, "w");
-  if (output_file3 == NULL) {
-    printf("Error opening file!\n");
-    exit(1);
-  }
-
-  gamma = 5.0 / 3.0;
-
-  chemistry_init(params, &us, &internal_const, &chemistry_data);
-  chemistry_first_init_part(&p, &xp, &chemistry_data);
-  chemistry_print(&chemistry_data);
-
-  cosmology_init(params, &us, &internal_const, &cosmo);
-  cosmology_print(&cosmo);
-
-  u = 1.0 * pow(10.0, 11) /
-      (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
-       units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
-  pressure = u * p.rho * (gamma - 1.0);
-  hydrogen_number_density =
-      hydrogen_number_density_cgs *
-      pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
-  p.rho = hydrogen_number_density * internal_const.const_proton_mass *
-          (1.0 + p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
-                     p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-  p.entropy = pressure / (pow(p.rho, gamma));
-
-  cooling_init(params, &us, &internal_const, &cooling);
-  cooling_print(&cooling);
-
-  // construct 1d table of cooling rates wrt temperature
-  float H_plus_He_heat_table[176];                // WARNING sort out how it is
-                                                  // declared/allocated
-  float H_plus_He_electron_abundance_table[176];  // WARNING sort out how it is
-                                                  // declared/allocated
-  float temp_table[176];  // WARNING sort out how it is declared/allocated
-  float element_cooling_table[9 * 176];         // WARNING sort out how it is
-                                                // declared/allocated
-  float element_electron_abundance_table[176];  // WARNING sort out how it is
-                                                // declared/allocated
-  construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
-                             cooling.table.element_cooling.H_plus_He_heating,
-                             H_plus_He_heat_table);
-  construct_1d_table_from_4d(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.H_plus_He_electron_abundance,
-      H_plus_He_electron_abundance_table);
-  construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
-                             cooling.table.element_cooling.temperature,
-                             temp_table);
-  construct_1d_table_from_4d_elements(
-      &p, &cooling, &cosmo, &internal_const,
-      cooling.table.element_cooling.metal_heating, element_cooling_table);
-  construct_1d_table_from_3d(&p, &cooling, &cosmo, &internal_const,
-                             cooling.table.element_cooling.electron_abundance,
-                             element_electron_abundance_table);
-
-  //
-  // printf("cooling table values \n");
-  // for( int j=0; j < 176; j++) {
-  //    printf("   %.5e",  H_plus_He_heat_table[j]+element_cooling_table[j]
-
-  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
-                                             &internal_const) *
-                cooling.number_density_scale;
-  ratefact = inn_h * (XH / eagle_proton_mass_cgs);
-  u_ini_cgs = pow(10.0, 14);
-
-  // Compute contributions to cooling rate from different metals
-  // for(int t_i = 0; t_i < n_t_i; t_i++){
-  //
-  //  u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
-  //  u =
-  //  u_cgs/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
-  //  pressure = u*p.rho*(gamma -1.0);
-  //  hydrogen_number_density =
-  //  hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
-  //  p.rho =
-  //  hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-  //  p.entropy = pressure/(pow(p.rho,gamma));
-
-  //  //cooling_du_dt =
-  //  eagle_print_metal_cooling_rate(&p,&cooling,&cosmo,&internal_const);
-  //  cooling_du_dt =
-  //  eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cooling,&cosmo,&internal_const);
-  //  temperature_cgs =
-  //  eagle_convert_u_to_temp(u_cgs,&p,&cooling,&cosmo,&internal_const);
-  //  fprintf(output_file,"%.5e %.5e\n",temperature_cgs,cooling_du_dt);
-  //  fprintf(output_file2,"%.5e
-  //  %.5e\n",u*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
-  //  temperature_cgs);
-
-  //  newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
-  //  fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
-  //  //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e
-  //  \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
-  //}
-  // fclose(output_file);
-  // fclose(output_file2);
-  // fclose(output_file3);
-
-  double dLambdaNet_du, LambdaNet;  //, LambdaNext;
-  float x_init, u_eq = 2.0e12;
-  for (int j = 5; j < 6; j++) {
-    float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0, 0.5 * (j + 5)),
-                                                   temp_table, &p, &cooling,
-                                                   &cosmo, &internal_const),
-          x, du;
-    float dt = 2.0e-4 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
-    LambdaNet = eagle_cooling_rate_1d_table(
-        u_ini, &dLambdaNet_du, H_plus_He_heat_table,
-        H_plus_He_electron_abundance_table, element_cooling_table,
-        element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
-        &internal_const);
-    float u_temp = u_ini + LambdaNet * ratefact * dt;
-    /* RGB removed this **
-    if (u_temp > 0) LambdaNext = eagle_cooling_rate_1d_table(u_temp,
-    &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
-    element_cooling_table, element_electron_abundance_table, temp_table, &p,
-    &cooling, &cosmo, &internal_const);
-    if (fabs(LambdaNet - LambdaNext)/LambdaNet < 0.5) {
-      u_temp = u_ini;
-    } else {
-      u_temp =
-    eagle_convert_temp_to_u_1d_table(1.0e4,temp_table,&p,&cooling,&cosmo,&internal_const);
-                } */
-    if (u_temp > u_eq) {
-      x_init = log(u_temp);
-    } else {
-      x_init = log(u_eq);
-    }
-    x = newton_iter(x_init, u_ini, H_plus_He_heat_table,
-                    H_plus_He_electron_abundance_table, element_cooling_table,
-                    element_electron_abundance_table, temp_table, &p, &cosmo,
-                    &cooling, &internal_const, dt);
-    printf(
-        "testing newton integration, u_ini, u %.5e %.5e, temperature initial, "
-        "final %.5e %.5e\n",
-        u_ini, exp(x),
-        eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
-                                         &cosmo, &internal_const),
-        eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
-                                         &cosmo, &internal_const));
-  }
-
-  free(params);
-
-  return 0;
-}
diff --git a/tests/testCoolingIntegration.c b/tests/testCoolingIntegration.c
deleted file mode 100644
index 925ca8a6f7889e8f29112514a538bf7fac8d0759..0000000000000000000000000000000000000000
--- a/tests/testCoolingIntegration.c
+++ /dev/null
@@ -1,265 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-int main() {
-
-  /* Read the parameter file */
-  struct swift_params *params = malloc(sizeof(struct swift_params));
-  struct unit_system us;
-  struct chemistry_data chemistry_data;
-  struct part p;
-  struct xpart xp;
-  struct phys_const internal_const;
-  struct cooling_function_data cooling;
-  struct cosmology cosmo;
-  // char *parametersFileName = "../examples/CoolingBox/coolingBox.yml";
-  char *parametersFileName = "../examples/EAGLE_12/eagle_12.yml";
-  enum table_index {
-    EAGLE_Hydrogen = 0,
-    EAGLE_Helium,
-    EAGLE_Carbon,
-    EAGLE_Nitrogen,
-    EAGLE_Oxygen,
-    EAGLE_Neon,
-    EAGLE_Magnesium,
-    EAGLE_Silicon,
-    EAGLE_Iron
-  };
-
-  if (params == NULL) error("Error allocating memory for the parameter file.");
-  message("Reading runtime parameters from file '%s'", parametersFileName);
-  parser_read_file(parametersFileName, params);
-
-  char output_filename[40];
-  FILE **output_file = malloc(10 * sizeof(FILE *));
-
-  /* And dump the parameters as used. */
-  parser_write_params_to_file(params, "used_parameters.yml");
-
-  units_init(&us, params, "InternalUnitSystem");
-  phys_const_init(&us, params, &internal_const);
-
-  double hydrogen_number_density_cgs = 1.582e-3;
-  // double hydrogen_number_density_cgs = 1.0e-1;
-  double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
-      temperature_cgs, newton_func, ratefact;
-  u_cgs = 0;
-  cooling_du_dt = 0;
-  temperature_cgs = 0;
-  newton_func = 0;
-  // int n_t_i = 2000;
-
-  gamma = 5.0 / 3.0;
-
-  chemistry_init(params, &us, &internal_const, &chemistry_data);
-  chemistry_first_init_part(&p, &xp, &chemistry_data);
-  chemistry_print(&chemistry_data);
-
-  cosmology_init(params, &us, &internal_const, &cosmo);
-  cosmology_print(&cosmo);
-
-  float u_ini = 3.357e15;
-  u = u_ini / (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
-               units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
-  pressure = u * p.rho * (gamma - 1.0);
-  hydrogen_number_density =
-      hydrogen_number_density_cgs *
-      pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
-  p.rho = hydrogen_number_density * internal_const.const_proton_mass *
-          (1.0 + p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
-                     p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-  p.entropy = pressure / (pow(p.rho, gamma));
-  cosmo.z = 0.0999744;
-
-  cooling_init(params, &us, &internal_const, &cooling);
-  cooling_print(&cooling);
-
-  // construct 1d table of cooling rates wrt temperature
-  float H_plus_He_heat_table[176];                // WARNING sort out how it is
-                                                  // declared/allocated
-  float H_plus_He_electron_abundance_table[176];  // WARNING sort out how it is
-                                                  // declared/allocated
-  float temp_table[176];  // WARNING sort out how it is declared/allocated
-  float element_cooling_table[176];             // WARNING sort out how it is
-                                                // declared/allocated
-  float element_electron_abundance_table[176];  // WARNING sort out how it is
-                                                // declared/allocated
-  for (int k = 0; k < cooling.N_Temp; k++) H_plus_He_heat_table[k] = 0.0;
-  for (int k = 0; k < cooling.N_Temp; k++)
-    H_plus_He_electron_abundance_table[k] = 0.0;
-  for (int k = 0; k < cooling.N_Temp; k++) temp_table[k] = 0.0;
-  for (int k = 0; k < cooling.N_Temp; k++) element_cooling_table[k] = 0.0;
-  for (int k = 0; k < cooling.N_Temp; k++)
-    element_electron_abundance_table[k] = 0.0;
-  // construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_heating,H_plus_He_heat_table);
-
-  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
-  // float inn_h =
-  // chemistry_get_number_density(&p,&cosmo,chemistry_element_H,&internal_const)*cooling.number_density_scale;
-  float inn_h = hydro_get_physical_density(&p, &cosmo) *
-                units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY) * XH /
-                eagle_proton_mass_cgs;
-  ratefact = inn_h * (XH / eagle_proton_mass_cgs);
-  printf("XH inn_h ratefact %.5e %.5e %.5e\n", XH, inn_h, ratefact);
-
-  double dLambdaNet_du, LambdaNet;
-  float x_init;
-  for (int j = 0; j < 1; j++) {
-    // float u_ini =
-    // eagle_convert_temp_to_u_1d_table(pow(10.0,0.5*(j+5)),temp_table,&p,&cooling,&cosmo,&internal_const);
-    // float u_ini =  3.357e15;
-    u = u_ini / (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
-                 units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
-    pressure = u * p.rho * (gamma - 1.0);
-    hydrogen_number_density =
-        hydrogen_number_density_cgs *
-        pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
-    p.rho = hydrogen_number_density * internal_const.const_proton_mass *
-            (1.0 + p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
-                       p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
-    p.entropy = pressure / (pow(p.rho, gamma));
-
-    float x, du;
-    double u_temp;
-    // float dt = 2.0e-4*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
-    float dt = 1.73798e-3 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
-    construct_1d_table_from_4d_H_He(
-        &p, &cooling, &cosmo, &internal_const,
-        cooling.table.element_cooling.H_plus_He_heating, H_plus_He_heat_table,
-        &u_temp, u_ini, dt);
-    construct_1d_table_from_4d(
-        &p, &cooling, &cosmo, &internal_const,
-        cooling.table.element_cooling.H_plus_He_electron_abundance,
-        H_plus_He_electron_abundance_table);
-    construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
-                               cooling.table.element_cooling.temperature,
-                               temp_table);
-    construct_1d_table_from_4d_elements(
-        &p, &cooling, &cosmo, &internal_const,
-        cooling.table.element_cooling.metal_heating, element_cooling_table);
-    construct_1d_table_from_3d(&p, &cooling, &cosmo, &internal_const,
-                               cooling.table.element_cooling.electron_abundance,
-                               element_electron_abundance_table);
-
-    sprintf(output_filename, "%s%d%s", "cooling_integration_output_", j,
-            ".dat");
-    output_file[j] = fopen(output_filename, "w");
-    if (output_file[j] == NULL) {
-      printf("Error opening file!\n");
-      exit(1);
-    }
-    for (int k = 0; k < cooling.N_Temp - 1; k++) {
-      float lambda1 =
-          H_plus_He_heat_table[k] + element_cooling_table[k] *
-                                        H_plus_He_electron_abundance_table[k] /
-                                        element_electron_abundance_table[k];
-      float lambda2 = H_plus_He_heat_table[k + 1] +
-                      element_cooling_table[k + 1] *
-                          H_plus_He_electron_abundance_table[k + 1] /
-                          element_electron_abundance_table[k + 1];
-      // printf("temperature %.5e, internal energy
-      // %.5e\n",pow(10.0,temp_table[k]),
-      // eagle_convert_temp_to_u_1d_table(pow(10.0,temp_table[k]),temp_table,&p,&cooling,&cosmo,&internal_const));
-      float u2 = eagle_convert_temp_to_u_1d_table(pow(10.0, temp_table[k + 1]),
-                                                  temp_table, &p, &cooling,
-                                                  &cosmo, &internal_const);
-      float u1 = eagle_convert_temp_to_u_1d_table(pow(10.0, temp_table[k]),
-                                                  temp_table, &p, &cooling,
-                                                  &cosmo, &internal_const);
-      float delta_u = u2 - u1;
-      // fprintf(output_file[j], "%.5e %.5e %.5e\n", pow(10.0,temp_table[k]), u1
-      // - u_ini - lambda1*ratefact*dt, (lambda2 - lambda1)/delta_u);
-      fprintf(output_file[j], "%.5e %.5e %.5e\n", u1,
-              u1 - u_ini - lambda1 * ratefact * dt,
-              (lambda2 - lambda1) / delta_u);
-    }
-    fclose(output_file[j]);
-
-    LambdaNet = eagle_cooling_rate_1d_table(
-        u_ini, &dLambdaNet_du, H_plus_He_heat_table,
-        H_plus_He_electron_abundance_table, element_cooling_table,
-        element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
-        &internal_const);
-    double u_eq = 5.0e12;
-    double u_temp_guess = u_ini + LambdaNet * ratefact * dt;
-    printf(
-        "u_guess, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du %.5e %.5e %.5e "
-        "%.5e %.5e %.5e %.5e \n",
-        u_temp, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du, ratefact, dt);
-
-    // if ((LambdaNet < 0 && u_temp < u_temp_guess) ||
-    //    (LambdaNet >= 0 && u_temp > u_temp_guess))
-    //  u_temp = u_temp_guess;
-    u_temp = u_temp_guess;
-    if ((LambdaNet < 0 && u_temp < u_eq) || (LambdaNet >= 0 && u_temp > u_eq))
-      u_temp = u_eq;
-
-    x_init = log(u_temp);
-    // int printflag = 1;
-    // x =
-    // newton_iter(x_init,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cosmo,&cooling,&internal_const,dt,&printflag);
-    x = bisection_iter(x_init, u_ini, H_plus_He_heat_table,
-                       H_plus_He_electron_abundance_table,
-                       element_cooling_table, element_electron_abundance_table,
-                       temp_table, &p, &cosmo, &cooling, &internal_const, dt);
-    // printf("testing newton integration, u_ini, u %.5e %.5e, temperature
-    // initial, final %.5e %.5e\n", u_ini, exp(x),
-    // eagle_convert_u_to_temp_1d_table(u_ini,&du,temp_table,&p,&cooling,&cosmo,&internal_const),
-    // eagle_convert_u_to_temp_1d_table(exp(x),&du,temp_table,&p,&cooling,&cosmo,&internal_const));
-
-    u = u_ini;
-    double u_next;
-    int nt = 10000;
-    float dt_sub = dt / nt;
-    for (int t = 0; t < nt; t++) {
-      LambdaNet = eagle_cooling_rate_1d_table(
-          u, &dLambdaNet_du, H_plus_He_heat_table,
-          H_plus_He_electron_abundance_table, element_cooling_table,
-          element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
-          &internal_const);
-      u_next = u + LambdaNet * ratefact * dt_sub;
-      printf(
-          "here u_next u lambda_net ratefact dt_sub, du t %.5e %.5e %.5e %.5e "
-          "%.5e %.5e %d\n",
-          u_next, u, LambdaNet, ratefact, dt_sub, LambdaNet * ratefact * dt_sub,
-          t);
-      u = u_next;
-    }
-    printf(
-        "testing newton integration, u_ini, u, u subcycle %.5e %.5e %.5e, "
-        "temperature initial, final, subcycled %.5e %.5e %.5e\n",
-        u_ini, exp(x), u,
-        eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
-                                         &cosmo, &internal_const),
-        eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
-                                         &cosmo, &internal_const),
-        eagle_convert_u_to_temp_1d_table(u, &du, temp_table, &p, &cooling,
-                                         &cosmo, &internal_const));
-  }
-
-  free(params);
-
-  return 0;
-}