diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
index 199e6f25f2297d3f989145067af37c1fde4468b3..c886c58d032b662724df4da81f567705033c3f15 100644
--- a/examples/CoolingBox/analytical_test.py
+++ b/examples/CoolingBox/analytical_test.py
@@ -57,8 +57,8 @@ T_init = 1.0e7
 
 # Read the initial state of the gas
 f = h5.File(snap_filename,'r')
-rho = np.mean(f["/PartType0/Density"])
-pressure = np.mean(f["/PartType0/Pressure"])
+#rho = np.mean(f["/PartType0/Density"])
+#pressure = np.mean(f["/PartType0/Pressure"])
 
 # Read the units parameters from the snapshot
 units = f["InternalCodeUnits"]
@@ -79,14 +79,14 @@ u_snapshots_cgs = zeros(nsnap)
 u_part_snapshots_cgs = zeros((nsnap,npart))
 t_snapshots_cgs = zeros(nsnap)
 scale_factor = zeros(nsnap)
+rho = zeros(nsnap)
 for i in range(nsnap):
     snap = h5.File("coolingBox_%0.4d.hdf5"%i,'r')
     u_part_snapshots_cgs[i][:] = snap["/PartType0/InternalEnergy"][:]*(unit_length**2)/(unit_time**2)
     t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
     scale_factor[i] = snap["/Header"].attrs["Scale-factor"]
+    rho[i] = np.mean(snap["/PartType0/Density"])*unit_mass/(unit_length**3)/(scale_factor[i]**3)
 
-print(rho, rho*unit_mass/(unit_length**3)/(scale_factor[0]**3))
-rho = rho*unit_mass/(unit_length**3)/(scale_factor[0]**3)
 
 # calculate analytic solution. Since cooling rate is constant,
 # temperature decrease in linear in time.
@@ -101,8 +101,7 @@ for line in file_in:
 
 tfinal = t_snapshots_cgs[nsnap-1]
 tfirst = t_snapshots_cgs[0]
-#nt = nsnap*10
-nt = 128
+nt = nsnap
 dt = (tfinal-tfirst)/nt
 
 t_sol = np.zeros(int(nt))
@@ -116,11 +115,12 @@ for j in range(int(nt-1)):
 	t_sol[j+1] = t_sol[j] + dt
 	Lambda_net = interpol_lambda(internal_energy,cooling_rate,u_sol[j])
 	if int(sys.argv[4]) == 1:
-		nH = 0.702*rho/m_p
+		nH = 0.702*rho[j]/m_p
 		ratefact = nH*0.702/m_p
 	else:
-		nH = 0.752*rho/m_p
+		nH = 0.752*rho[j]/m_p
 		ratefact = nH*0.752/m_p
+	print(rho[j],nH,ratefact,u,Lambda_net)
 	u_next = u - Lambda_net*ratefact*dt
 	print(u_next, u, ratefact, Lambda_net, dt, nH)
 	u_sol[j+1] = u_next
diff --git a/examples/CoolingBox/coolingBox_non-cosmo.yml b/examples/CoolingBox/coolingBox_non-cosmo.yml
index d8627e8bd2072bfcb20e7abd5aedda7b336566cb..2c10909022f0291bf60b4424a53ed27efa47d98a 100644
--- a/examples/CoolingBox/coolingBox_non-cosmo.yml
+++ b/examples/CoolingBox/coolingBox_non-cosmo.yml
@@ -31,6 +31,7 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./coolingBox.hdf5     # The file to read
+  periodic:   1
 
 # Dimensionless pre-factor for the time-step condition
 LambdaCooling:
diff --git a/examples/CoolingBox/coolingBox_primordial.yml b/examples/CoolingBox/coolingBox_primordial.yml
index 6b4875bdae5e8ff0bcb8a27889eae4129f3f6245..e352dbdb9f1469dfd939c2ff5fb1431f68b9cc5f 100644
--- a/examples/CoolingBox/coolingBox_primordial.yml
+++ b/examples/CoolingBox/coolingBox_primordial.yml
@@ -20,7 +20,7 @@ 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-6  # The maximal time-step size of the simulation (in internal units).
+  dt_max:     1e-9  # The maximal time-step size of the simulation (in internal units).
   
 Scheduler:
   max_top_level_cells:    15
@@ -55,6 +55,7 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./coolingBox.hdf5     # The file to read
+  periodic:   1
 
 # Dimensionless pre-factor for the time-step condition
 LambdaCooling:
diff --git a/examples/CoolingBox/test.sh b/examples/CoolingBox/test.sh
index ff635f83ed546dc052632f72fa9345770120426e..8f746d22f194fc40502c24d693e50a75e3dd1ee1 100755
--- a/examples/CoolingBox/test.sh
+++ b/examples/CoolingBox/test.sh
@@ -22,15 +22,20 @@ for z in 0.1; do
         
       # set starting, ending redshift, how often to write to file
       a_begin=$(python -c "print 1.0/(1.0+$z)")
-      a_end=$(python -c "print min(1.0, $a_begin*1.001)")
+      a_end=$(python -c "print min(1.0, $a_begin*1.0001)")
       delta_a=$(python -c "print 1.0 + ($a_end/$a_begin - 1.0)/50.")
   
       # change hydrogen number density
-      nh=$(python -c "print 10.0**$nh_exp/(1.0+$z)**3")
+      if [ $solar == 0 ]; 
+      then
+        rho=$(python -c "print 10.0**$nh_exp/0.7*1.6726e-24")
+      else
+        rho=$(python -c "print 10.0**$nh_exp/0.752*1.6726e-24")
+      fi
       for pressure_index in 8; do
         pressure=$(python -c "print 6.68e-14")
         
-	python makeIC.py # $nh $pressure
+	python makeIC.py $rho $pressure
 	rm coolingBox_*hdf5
   
         # run cooling box
diff --git a/examples/CoolingBox/testCooling.c b/examples/CoolingBox/testCooling.c
index 148c36cee0fb185c9a770245b62021974b85f8e7..2a5969e6464befa3681427175e94f356f8575dfd 100644
--- a/examples/CoolingBox/testCooling.c
+++ b/examples/CoolingBox/testCooling.c
@@ -30,11 +30,14 @@
  * hydrogen number density and internal energy specified.
  *
  * @param p Particle data structure
+ * @param xp extra particle structure
+ * @param us unit system struct
  * @param cooling Cooling function data structure
  * @param cosmo Cosmology data structure
  * @param phys_const Physical constants data structure
- * @param nh Hydrogen number density (cgs units)
+ * @param nh_cgs Hydrogen number density (cgs units)
  * @param u Internal energy (cgs units)
+ * @param ti_current integertime to set cosmo quantities
  */
 void set_quantities(struct part *restrict p,
 		    struct xpart *restrict xp,
diff --git a/examples/CoolingRates/cooling_rates_plot.py b/examples/CoolingRates/cooling_rates_plot.py
index 1cd372f05e3d1fb6d7d2db92db84ec76deca0176..4ae1f4a19ddea7fd106f3b2fd540081a663b092b 100644
--- a/examples/CoolingRates/cooling_rates_plot.py
+++ b/examples/CoolingRates/cooling_rates_plot.py
@@ -5,125 +5,47 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-k_in_J_K = 1.38064852e-23
-mH_in_kg = 1.6737236e-27
-erg_to_J = 1.0e-7
-gas_gamma=5.0/3.0
-
-# Primoridal ean molecular weight as a function of temperature
-def mu(T):
-    H_frac=0.752
-    T_trans=10000.0
-    if T > T_trans:
-        return 4. / (8. - 5. * (1. - H_frac))
-    else:
-        return 4. / (1. + 3. * H_frac)
-
-# Temperature of some primoridal gas with a given internal energy
-def T(u):
-    H_frac=0.752
-    T_trans=10000.0
-    T_over_mu = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
-    ret = np.ones(np.size(u)) * T_trans
-
-    # Enough energy to be ionized?
-    mask_ionized = (T_over_mu > (T_trans+1) / mu(T_trans+1))
-    if np.sum(mask_ionized)  > 0:
-        ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans*10)
-
-    # Neutral gas?
-    mask_neutral = (T_over_mu < (T_trans-1) / mu((T_trans-1)))
-    if np.sum(mask_neutral)  > 0:
-        ret[mask_neutral] = T_over_mu[mask_neutral] * mu(0)
-
-    return ret
-
 # Number of metals tracked by EAGLE cooling
 elements = 11
 
 # Declare arrays of internal energy and cooling rate
-u = [[] for i in range(elements+1)]
+u = []
 cooling_rate = [[] for i in range(elements+1)]
 Temperature = [[] for i in range(elements+1)]
 
 # Read in total cooling rate
-#file_in = open('cooling_output.dat', 'r')
-#for line in file_in:
-#	data = line.split()
-#	u.append(float(data[0]))
-#	cooling_rate[0].append(-float(data[1]))
-#
-#file_in.close()
+file_in = open('cooling_output.dat', 'r')
+for line in file_in:
+	data = line.split()
+	u.append(float(data[0]))
+	cooling_rate[0].append(-float(data[1]))
+file_in.close()
 
-n_iz = 10
-for iz in range(n_iz):
-	file_in = open('cooling_output_'+str(iz)+'.dat','r')
-	z = float(file_in.readline())
+# Read in contributions to cooling rates from each of the elements
+for elem in range(elements):
+	file_in = open('cooling_element_'+str(elem)+'.dat','r')
 	for line in file_in:
-		data = line.split() 
-		u[iz+1].append(float(data[0]))
-		cooling_rate[iz+1].append(-float(data[1]))
+        	data = line.split()
+        	cooling_rate[elem+1].append(-float(data[0]))
 	file_in.close()
-	a = 1.0/(1.0 + z)
-	#u[iz+1] = np.asarray(u[iz+1]) / a**(3 * (gas_gamma - 1.))
-	#Temperature[iz+1] = T(u[iz+1] * erg_to_J)
 
 # Plot
 ax = plt.subplot(111)
-#p0, = plt.loglog(u, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1,  = plt.loglog(u[1], cooling_rate[1], linewidth = 0.5, color = 'k',             label = 'z = 10')
-p2,  = plt.loglog(u[2], cooling_rate[2], linewidth = 0.5, color = 'b',             label = 'z = 9')
-p3,  = plt.loglog(u[3], cooling_rate[3], linewidth = 0.5, color = 'g',             label = 'z = 8')
-p4,  = plt.loglog(u[4], cooling_rate[4], linewidth = 0.5, color = 'r',             label = 'z = 7')
-p5,  = plt.loglog(u[5], cooling_rate[5], linewidth = 0.5, color = 'c',             label = 'z = 6')
-p6,  = plt.loglog(u[6], cooling_rate[6], linewidth = 0.5, color = 'm',             label = 'z = 5')
-p7,  = plt.loglog(u[7], cooling_rate[7], linewidth = 0.5, color = 'y',             label = 'z = 4')
-p8,  = plt.loglog(u[8], cooling_rate[8], linewidth = 0.5, color = 'lightgray',     label = 'z = 3')
-p9,  = plt.loglog(u[9], cooling_rate[9], linewidth = 0.5, color = 'olive',         label = 'z = 2')
-p10, = plt.loglog(u[10], cooling_rate[10], linewidth = 0.5, color = 'saddlebrown', label = 'z = 1')
-p11, = plt.loglog(u[1],  -np.asarray(cooling_rate[1]),  linewidth = 0.5, linestyle = '--', color = 'k')
-p12, = plt.loglog(u[2],  -np.asarray(cooling_rate[2]),  linewidth = 0.5, linestyle = '--', color = 'b')
-p13, = plt.loglog(u[3],  -np.asarray(cooling_rate[3]),  linewidth = 0.5, linestyle = '--', color = 'g')
-p14, = plt.loglog(u[4],  -np.asarray(cooling_rate[4]),  linewidth = 0.5, linestyle = '--', color = 'r')
-p15, = plt.loglog(u[5],  -np.asarray(cooling_rate[5]),  linewidth = 0.5, linestyle = '--', color = 'c')
-p16, = plt.loglog(u[6],  -np.asarray(cooling_rate[6]),  linewidth = 0.5, linestyle = '--', color = 'm')
-p17, = plt.loglog(u[7],  -np.asarray(cooling_rate[7]),  linewidth = 0.5, linestyle = '--', color = 'y')
-p18, = plt.loglog(u[8],  -np.asarray(cooling_rate[8]),  linewidth = 0.5, linestyle = '--', color = 'lightgray')
-p19, = plt.loglog(u[9],  -np.asarray(cooling_rate[9]),  linewidth = 0.5, linestyle = '--', color = 'olive')
-p20, = plt.loglog(u[10], -np.asarray(cooling_rate[10]), linewidth = 0.5, linestyle = '--', color = 'saddlebrown')
+p0, = plt.loglog(u, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
+p1, = plt.loglog(u, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
+p2, = plt.loglog(u, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
+p3, = plt.loglog(u, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
+p4, = plt.loglog(u, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
+p5, = plt.loglog(u, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
+p6, = plt.loglog(u, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
+p7, = plt.loglog(u, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
+p8, = plt.loglog(u, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
+p9, = plt.loglog(u, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
+p10, = plt.loglog(u, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
 ax.set_position([0.15,0.15,0.75,0.75])
-#plt.xlim([1e10,1e17])
-#plt.ylim([1e-24,1e-21])
-plt.xlabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
+plt.xlim([1e3,1e8])
+plt.ylim([1e-24,1e-21])
+plt.xlabel("Temperature ${\\rm{[K]}}$", fontsize = 14)
 plt.ylabel("${\Lambda/n_h^2 }$ ${\\rm{[erg \cdot cm^3 \cdot s^{-1}]}}$", fontsize = 14)
-plt.legend(handles = [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
-plt.show()
-
-# Read in contributions to cooling rates from each of the elements
-#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()
-
-## Plot
-#ax = plt.subplot(111)
-#p0, = plt.loglog(u, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-#p1, = plt.loglog(u, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
-#p2, = plt.loglog(u, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
-#p3, = plt.loglog(u, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
-#p4, = plt.loglog(u, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
-#p5, = plt.loglog(u, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
-#p6, = plt.loglog(u, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
-#p7, = plt.loglog(u, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
-#p8, = plt.loglog(u, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
-#p9, = plt.loglog(u, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
-#p10, = plt.loglog(u, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
-#ax.set_position([0.15,0.15,0.75,0.75])
-#plt.xlim([1e12,1e17])
-#plt.ylim([1e-24,1e-21])
-#plt.xlabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
-#plt.ylabel("${\Lambda/n_h^2 }$ ${\\rm{[erg \cdot cm^3 \cdot s^{-1}]}}$", fontsize = 14)
-##plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
-#plt.show()
+plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
+plt.savefig("cooling_rates",dpi=200)
diff --git a/examples/CoolingRates/testCooling.c b/examples/CoolingRates/testCooling.c
index c7ea961f953169723cc078b4cba0a102eae59dbd..87b07b1520bb2eb683820c9bee86f00554e9a188 100644
--- a/examples/CoolingRates/testCooling.c
+++ b/examples/CoolingRates/testCooling.c
@@ -153,13 +153,11 @@ int main(int argc, char **argv) {
 
   // Calculate contributions from metals to cooling rate
   // open file
-  const char *output_filename = "cooling_output.dat";
-  FILE *output_file = fopen(output_filename, "w");
+  FILE *output_file = fopen("cooling_output.dat", "w");
   if (output_file == NULL) {
     printf("Error opening file!\n");
     exit(1);
   }
-  fprintf(output_file, "%.5e\n", cosmo.z);
 
   // set hydrogen number density
   if (log_10_nh == 100) {
@@ -186,16 +184,12 @@ int main(int argc, char **argv) {
     u = hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
         cooling.internal_energy_scale;
     float cooling_du_dt;
-    double dLambdaNet_du;
 
     // calculate cooling rates
-    cooling_du_dt = eagle_cooling_rate(
-        log(u), &dLambdaNet_du, n_h_i, d_n_h, He_i, d_He, &p, &cooling, &cosmo, &internal_const,
-        abundance_ratio);
     temperature = eagle_convert_u_to_temp(log10(u), &du, n_h_i, He_i, d_n_h, d_He, &cooling, &cosmo);
-    //cooling_du_dt = eagle_print_metal_cooling_rate(
-    //    n_h_i, d_n_h, He_i, d_He, &p, &xp, &cooling, &cosmo, &internal_const,
-    //    abundance_ratio);
+    cooling_du_dt = eagle_print_metal_cooling_rate(
+        n_h_i, d_n_h, He_i, d_He, &p, &xp, &cooling, &cosmo, &internal_const,
+        abundance_ratio);
     fprintf(output_file, "%.5e %.5e\n", exp(M_LN10*temperature), cooling_du_dt);
   }
   fclose(output_file);
diff --git a/examples/CoolingRates/testCooling.yml b/examples/CoolingRates/testCooling.yml
index 912cc05387f991b307c369dcf1233651a2858769..63f3e381d73f706a2dcd97cf7882f3136eb2e35a 100644
--- a/examples/CoolingRates/testCooling.yml
+++ b/examples/CoolingRates/testCooling.yml
@@ -58,18 +58,19 @@ 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
-  CalciumOverSilicon: 0.0941736
-  SulphurOverSilicon: 0.6054160
+  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
+
 
 EagleCooling:
   filename:                /cosma5/data/Eagle/BG_Tables/CoolingTables/