Merge branch 'master' of gitlab.cosma.dur.ac.uk:swift/swiftsim

.gitignore

@@ -62,6 +62,7 @@ examples/Cooling/CoolingRates/cooling_rates
 examples/Cooling/CoolingRates/cooling_element_*.dat
 examples/Cooling/CoolingRates/cooling_output.dat
 examples/SubgridTests/StellarEvolution/StellarEvolutionSolution*
+examples/SubgridTests/CosmologicalStellarEvolution/StellarEvolutionSolution*
 
 tests/testActivePair
 tests/testActivePair.sh

examples/SubgridTests/CosmologicalStellarEvolution/README

+Example for testing EAGLE stellar feedback. This consists of a uniform box of gas with a star in the center. The amount of feedback can then be checked by summing over the gas particles in the whole box and comparing to the expected amount of feedback from a single star over a given time period. 
+
+If only mass enrichment is of interest, the box can be run with ICs generated from a smaller glass (eg glassCube_32.hdf5). In this case, however it is necessary to turn off energy feedback (eg. by setting the return value of compute_SNe in src/stars/EAGLE/stars.h to zero) or using a larger glass (glassCube_64.hdf5). 
+
+Use the python script, plot_box_evolution.py to compare total mass evolution of gas particles in the whole box with what is expected based on EAGLE standalone feedback test.
+
+Use plot_paricle_evolution.py to plot the evolution of particles starting in the viscinity of the star at the beginning of the simulation. 

examples/SubgridTests/CosmologicalStellarEvolution/check_continuous_heating.py

+# Script for plotting energy evolution of uniform box of gas with single star in the 
+# centre when running with stochastic energy injection. It also checks that the change
+# in total energy of the gas particles is within a specified tolerance from what is 
+# expected based on the mass of the star particle (Note that this tolerance could be
+# somewhat high because of Poisson noise and the relatively small number of injection
+# events)
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+import numpy as np
+import glob
+
+# Number of snapshots and elements
+newest_snap_name = max(glob.glob('stellar_evolution_*.hdf5'), key=os.path.getctime)
+n_snapshots = int(newest_snap_name.replace('stellar_evolution_','').replace('.hdf5','')) + 1
+n_elements = 9
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left'    : 0.3,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.18,
+'figure.subplot.top'     : 0.92,
+'figure.subplot.wspace'  : 0.21,
+'figure.subplot.hspace'  : 0.19,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+star_initial_mass = sim["/PartType4/Masses"][0]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
+unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
+unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
+unit_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Calculate solar mass in internal units
+const_solar_mass = 1.98848e33 / unit_mass_in_cgs
+
+# Define Gyr
+Gyr_in_cgs = 1e9 * 365 * 24 * 3600.
+
+# Find out how many particles (gas and star) we have
+n_parts = sim["/Header"].attrs["NumPart_Total"][0]
+n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
+
+# Declare arrays for data
+masses = zeros((n_parts,n_snapshots))
+star_masses = zeros((n_sparts,n_snapshots))
+internal_energy = zeros((n_parts,n_snapshots))
+velocity_parts = zeros((n_parts,3,n_snapshots))
+time = zeros(n_snapshots)
+
+# Read fields we are checking from snapshots
+#for i in [0,n_snapshots-1]:
+for i in range(n_snapshots):
+	sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
+	print('reading snapshot '+str(i))
+	masses[:,i] = sim["/PartType0/Masses"]
+	internal_energy[:,i] = sim["/PartType0/InternalEnergies"]
+	velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
+	time[i] = sim["/Header"].attrs["Time"][0]
+
+# Check that the total amount of enrichment is as expected.
+# Define tolerance. Note, relatively high value used due to
+# Poisson noise associated with stochastic energy injection.
+eps = 0.15
+
+# Stochastic heating
+vel2 = zeros((n_parts,n_snapshots))
+vel2[:,:] = velocity_parts[:,0,:]*velocity_parts[:,0,:] + velocity_parts[:,1,:]*velocity_parts[:,1,:] + velocity_parts[:,2,:]*velocity_parts[:,2,:]
+total_kinetic_energy_cgs = np.sum(np.multiply(vel2,masses)*0.5,axis = 0) * unit_energy_in_cgs
+total_energy_cgs = np.sum(np.multiply(internal_energy,masses),axis = 0) * unit_energy_in_cgs
+total_energy_released_cgs = total_energy_cgs[n_snapshots-1] - total_energy_cgs[0] + total_kinetic_energy_cgs[n_snapshots-1] - total_kinetic_energy_cgs[0]
+
+# Calculate energy released
+energy_per_sn = 1.0e51 / unit_energy_in_cgs
+SNIa_efficiency = 2.e-3
+SNIa_timescale_Gyr = 2.0
+expected_energy_released_cgs = np.zeros(n_snapshots)
+for i in range(n_snapshots):
+	age_Gyr = time[i] * unit_time_in_cgs / Gyr_in_cgs
+	total_sn = SNIa_efficiency * (1.0 - np.exp(-age_Gyr/SNIa_timescale_Gyr)) / const_solar_mass
+	expected_energy_released_cgs[i] = total_sn * energy_per_sn * unit_energy_in_cgs
+
+# Did we get it right?
+if abs(total_energy_released_cgs - expected_energy_released_cgs[n_snapshots-1])/expected_energy_released_cgs[n_snapshots-1] < eps:
+	print("total stochastic energy release consistent with expectation. total stochastic energy release "+str(total_energy_released_cgs)+" expected "+ str(expected_energy_released_cgs[n_snapshots-1]) + " initial total internal energy "+ str(total_energy_cgs[0] + total_kinetic_energy_cgs[0]))
+else:
+	print("total stochastic energy release "+str(total_energy_released_cgs)+" expected "+ str(expected_energy_released_cgs[n_snapshots-1]) + " initial total internal energy "+ str(total_energy_cgs[0] + total_kinetic_energy_cgs[0]) + " energy change fraction of total " + str(total_energy_released_cgs/(total_energy_cgs[0]+total_kinetic_energy_cgs[0])))
+
+# Plot the energy evolution
+figure()
+subplot(111)
+plot(time*unit_time_in_cgs/Gyr_in_cgs, total_energy_cgs + total_kinetic_energy_cgs - total_energy_cgs[0] - total_kinetic_energy_cgs[0],color='k', linewidth=0.5, label="SWIFT")
+plot(time*unit_time_in_cgs/Gyr_in_cgs, expected_energy_released_cgs,color = 'r', linewidth=0.5, label="expected")
+xlabel("Time (Gyr)")
+ylabel("Total energy (erg)")
+legend()
+savefig("continuous_energy_evolution.png", dpi=200)
+

examples/SubgridTests/CosmologicalStellarEvolution/check_stellar_evolution.py

+import matplotlib
+
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+import numpy as np
+import glob
+
+# Number of snapshots and elements
+newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
+n_snapshots = (
+    int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
+)
+n_elements = 9
+
+# Plot parameters
+params = {
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 9,
+    "legend.fontsize": 9,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": True,
+    "figure.figsize": (3.15, 3.15),
+    "figure.subplot.left": 0.3,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.18,
+    "figure.subplot.top": 0.92,
+    "figure.subplot.wspace": 0.21,
+    "figure.subplot.hspace": 0.19,
+    "lines.markersize": 6,
+    "lines.linewidth": 2.0,
+    "text.latex.unicode": True,
+}
+
+rcParams.update(params)
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs[
+    "Hydrogen ionization transition temperature"
+][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+star_initial_mass = sim["/PartType4/Masses"][0]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
+unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
+unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
+unit_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Find out how many particles (gas and star) we have
+n_parts = sim["/Header"].attrs["NumPart_Total"][0]
+n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
+
+# Declare arrays for data
+masses = zeros((n_parts, n_snapshots))
+star_masses = zeros((n_sparts, n_snapshots))
+mass_from_AGB = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_AGB = zeros((n_parts, n_snapshots))
+mass_from_SNII = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNII = zeros((n_parts, n_snapshots))
+mass_from_SNIa = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+iron_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+metallicity = zeros((n_parts, n_snapshots))
+abundances = zeros((n_parts, n_elements, n_snapshots))
+internal_energy = zeros((n_parts, n_snapshots))
+coord_parts = zeros((n_parts, 3))
+velocity_parts = zeros((n_parts, 3, n_snapshots))
+coord_sparts = zeros(3)
+time = zeros(n_snapshots)
+
+# Read fields we are checking from snapshots
+for i in range(n_snapshots):
+    sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+    print("reading snapshot " + str(i))
+    abundances[:, :, i] = sim["/PartType0/ElementMassFractions"]
+    metallicity[:, i] = sim["/PartType0/Metallicities"]
+    masses[:, i] = sim["/PartType0/Masses"]
+    star_masses[:, i] = sim["/PartType4/Masses"]
+    mass_from_AGB[:, i] = sim["/PartType0/TotalMassFromAGB"]
+    metal_mass_frac_from_AGB[:, i] = sim["/PartType0/MetalMassFracFromAGB"]
+    mass_from_SNII[:, i] = sim["/PartType0/TotalMassFromSNII"]
+    metal_mass_frac_from_SNII[:, i] = sim["/PartType0/MetalMassFracFromSNII"]
+    mass_from_SNIa[:, i] = sim["/PartType0/TotalMassFromSNIa"]
+    metal_mass_frac_from_SNIa[:, i] = sim["/PartType0/MetalMassFracFromSNIa"]
+    iron_mass_frac_from_SNIa[:, i] = sim["/PartType0/IronMassFracFromSNIa"]
+    internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
+    velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
+    time[i] = sim["/Header"].attrs["Time"][0]
+
+# Define ejecta factor
+ejecta_factor = 1.0e-2
+ejecta_factor_metallicity = 1.0 - 2.0 / n_elements
+ejecta_factor_abundances = 1.0 / n_elements
+ejected_mass = star_initial_mass
+energy_per_SNe = 1.0e51 / unit_energy_in_cgs
+
+# Check that the total amount of enrichment is as expected.
+# Define tolerance
+eps = 0.01
+
+# Total mass
+total_part_mass = np.sum(masses, axis=0)
+if (
+    abs(
+        (total_part_mass[n_snapshots - 1] - total_part_mass[0]) / total_part_mass[0]
+        - ejected_mass / total_part_mass[0]
+    )
+    * total_part_mass[0]
+    / ejected_mass
+    < eps
+):
+    print("total mass released consistent with expectation")
+else:
+    print(
+        "mass increase "
+        + str(total_part_mass[n_snapshots - 1] / total_part_mass[0])
+        + " expected "
+        + str(1.0 + ejected_mass / total_part_mass[0])
+    )
+
+# Check that mass is conserved (i.e. total star mass decreases by same amount as total gas mass increases)
+total_spart_mass = np.sum(star_masses, axis=0)
+if (
+    abs(
+        (total_part_mass[n_snapshots - 1] + total_spart_mass[n_snapshots - 1])
+        / (total_part_mass[0] + total_spart_mass[0])
+        - 1.0
+    )
+    < eps ** 3
+):
+    print("total mass conserved")
+else:
+    print(
+        "initial part, spart mass "
+        + str(total_part_mass[0])
+        + " "
+        + str(total_spart_mass[0])
+        + " final mass "
+        + str(total_part_mass[n_snapshots - 1])
+        + " "
+        + str(total_spart_mass[n_snapshots - 1])
+    )
+
+# Total metal mass from AGB
+total_metal_mass_AGB = np.sum(np.multiply(metal_mass_frac_from_AGB, masses), axis=0)
+expected_metal_mass_AGB = ejecta_factor * ejected_mass
+if (
+    abs(total_metal_mass_AGB[n_snapshots - 1] - expected_metal_mass_AGB)
+    / expected_metal_mass_AGB
+    < eps
+):
+    print("total AGB metal mass released consistent with expectation")
+else:
+    print(
+        "total AGB metal mass "
+        + str(total_metal_mass_AGB[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass_AGB)
+    )
+
+# Total mass from AGB
+total_AGB_mass = np.sum(mass_from_AGB, axis=0)
+expected_AGB_mass = ejecta_factor * ejected_mass
+if abs(total_AGB_mass[n_snapshots - 1] - expected_AGB_mass) / expected_AGB_mass < eps:
+    print("total AGB mass released consistent with expectation")
+else:
+    print(
+        "total AGB mass "
+        + str(total_AGB_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_AGB_mass)
+    )
+
+# Total metal mass from SNII
+total_metal_mass_SNII = np.sum(np.multiply(metal_mass_frac_from_SNII, masses), axis=0)
+expected_metal_mass_SNII = ejecta_factor * ejected_mass
+if (
+    abs(total_metal_mass_SNII[n_snapshots - 1] - expected_metal_mass_SNII)
+    / expected_metal_mass_SNII
+    < eps
+):
+    print("total SNII metal mass released consistent with expectation")
+else:
+    print(
+        "total SNII metal mass "
+        + str(total_metal_mass_SNII[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass_SNII)
+    )
+
+# Total mass from SNII
+total_SNII_mass = np.sum(mass_from_SNII, axis=0)
+expected_SNII_mass = ejecta_factor * ejected_mass
+if (
+    abs(total_SNII_mass[n_snapshots - 1] - expected_SNII_mass) / expected_SNII_mass
+    < eps
+):
+    print("total SNII mass released consistent with expectation")
+else:
+    print(
+        "total SNII mass "
+        + str(total_SNII_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_SNII_mass)
+    )
+
+# Total metal mass from SNIa
+total_metal_mass_SNIa = np.sum(np.multiply(metal_mass_frac_from_SNIa, masses), axis=0)
+expected_metal_mass_SNIa = ejecta_factor * ejected_mass
+if (
+    abs(total_metal_mass_SNIa[n_snapshots - 1] - expected_metal_mass_SNIa)
+    / expected_metal_mass_SNIa
+    < eps
+):
+    print("total SNIa metal mass released consistent with expectation")
+else:
+    print(
+        "total SNIa metal mass "
+        + str(total_metal_mass_SNIa[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass_SNIa)
+    )
+
+# Total iron mass from SNIa
+total_iron_mass_SNIa = np.sum(np.multiply(iron_mass_frac_from_SNIa, masses), axis=0)
+expected_iron_mass_SNIa = ejecta_factor * ejected_mass
+if (
+    abs(total_iron_mass_SNIa[n_snapshots - 1] - expected_iron_mass_SNIa)
+    / expected_iron_mass_SNIa
+    < eps
+):
+    print("total SNIa iron mass released consistent with expectation")
+else:
+    print(
+        "total SNIa iron mass "
+        + str(total_iron_mass_SNIa[n_snapshots - 1])
+        + " expected "
+        + str(expected_iron_mass_SNIa)
+    )
+
+# Total mass from SNIa
+total_SNIa_mass = np.sum(mass_from_SNIa, axis=0)
+expected_SNIa_mass = ejecta_factor * ejected_mass
+if (
+    abs(total_SNIa_mass[n_snapshots - 1] - expected_SNIa_mass) / expected_SNIa_mass
+    < eps
+):
+    print("total SNIa mass released consistent with expectation")
+else:
+    print(
+        "total SNIa mass "
+        + str(total_SNIa_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_SNIa_mass)
+    )
+
+# Total metal mass
+total_metal_mass = np.sum(np.multiply(metallicity, masses), axis=0)
+expected_metal_mass = ejecta_factor_metallicity * ejected_mass
+if (
+    abs(total_metal_mass[n_snapshots - 1] - expected_metal_mass) / expected_metal_mass
+    < eps
+):
+    print("total metal mass released consistent with expectation")
+else:
+    print(
+        "total metal mass "
+        + str(total_metal_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass)
+    )
+
+# Total mass for each element
+expected_element_mass = ejecta_factor_abundances * ejected_mass
+for i in range(n_elements):
+    total_element_mass = np.sum(np.multiply(abundances[:, i, :], masses), axis=0)
+    if (
+        abs(total_element_mass[n_snapshots - 1] - expected_element_mass)
+        / expected_element_mass
+        < eps
+    ):
+        print(
+            "total element mass released consistent with expectation for element "
+            + str(i)
+        )
+    else:
+        print(
+            "total element mass "
+            + str(total_element_mass[n_snapshots - 1])
+            + " expected "
+            + str(expected_element_mass)
+            + " for element "
+            + str(i)
+        )
+

examples/SubgridTests/CosmologicalStellarEvolution/check_stochastic_heating.py

+# Script for plotting energy evolution of uniform box of gas with single star in the 
+# centre when running with stochastic energy injection. It also checks that the change
+# in total energy of the gas particles is within a specified tolerance from what is 
+# expected based on the mass of the star particle (Note that this tolerance could be
+# somewhat high because of Poisson noise and the relatively small number of injection
+# events)
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+import numpy as np
+import glob
+
+# Number of snapshots and elements
+newest_snap_name = max(glob.glob('stellar_evolution_*.hdf5'), key=os.path.getctime)
+n_snapshots = int(newest_snap_name.replace('stellar_evolution_','').replace('.hdf5','')) + 1
+n_elements = 9
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left'    : 0.3,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.18,
+'figure.subplot.top'     : 0.92,
+'figure.subplot.wspace'  : 0.21,
+'figure.subplot.hspace'  : 0.19,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+star_initial_mass = sim["/PartType4/Masses"][0]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
+unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
+unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
+unit_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Calculate solar mass in internal units
+const_solar_mass = 1.98848e33 / unit_mass_in_cgs
+
+# Define Gyr
+Gyr_in_cgs = 1e9 * 365 * 24 * 3600.
+
+# Find out how many particles (gas and star) we have
+n_parts = sim["/Header"].attrs["NumPart_Total"][0]
+n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
+
+# Declare arrays for data
+masses = zeros((n_parts,n_snapshots))
+star_masses = zeros((n_sparts,n_snapshots))
+internal_energy = zeros((n_parts,n_snapshots))
+velocity_parts = zeros((n_parts,3,n_snapshots))
+time = zeros(n_snapshots)
+
+# Read fields we are checking from snapshots