correct temperature calculation

auriga_comparison.py

@@ -93,9 +93,9 @@ for dir in sorted(root_dir.glob("*")):
         if levelmax != 11:
             continue
 
-    input_file = dir / "output_0007.hdf5"
+    input_file = dir / "output_0009.hdf5"
     if mode == Mode.richings:
-        input_file = dir / "output_0002.hdf5"
+        input_file = dir / "output_0004.hdf5"
     if is_by_adrian or is_ramses:
         input_file = dir / "output_0000.hdf5"
         softening_length = None

fix_hdf5_masses.py

@@ -10,8 +10,6 @@ from utils import print_progress
 gamma = 5 / 3
 YHe = 0.245421
 Tcmb0 = 2.7255
-const_proton_mass_cgs = 1.67262192369e-24
-const_boltzmann_k_cgs = 1.380649e-16
 
 
 def calculate_gas_internal_energy(omegab, hubble_param_, zstart_):
@@ -22,7 +20,7 @@ def calculate_gas_internal_energy(omegab, hubble_param_, zstart_):
         npol = 1.0
     unitv = 1e5
     adec = 1.0 / (
-        160.0 * (omegab * hubble_param_ * hubble_param_ / 0.022) ** (2.0 / 5.0)
+            160.0 * (omegab * hubble_param_ * hubble_param_ / 0.022) ** (2.0 / 5.0)
     )
     if astart_ < adec:
         Tini = Tcmb0 / astart_
@@ -84,12 +82,24 @@ hydrogen_mass_function = 1 - const_primordial_He_fraction_cgs
 mu_neutral = 4.0 / (1.0 + 3.0 * hydrogen_mass_function)
 mu_ionised = 4.0 / (8.0 - 5.0 * (1.0 - hydrogen_mass_function))
 T_transition = 1.0e4
+UnitMass_in_cgs = 1.98848e43  # 10^10 M_sun in grams
+UnitLength_in_cgs = 3.08567758e24  # 1 Mpc in centimeters
+UnitVelocity_in_cgs = 1e5  # 1 km/s in centimeters per second
+UnitTime_in_cgs = UnitLength_in_cgs / UnitVelocity_in_cgs
+const_proton_mass_cgs = 1.67262192369e-24
+const_boltzmann_k_cgs = 1.380649e-16
+
+const_proton_mass = const_proton_mass_cgs / UnitMass_in_cgs
+const_boltzmann_k = const_boltzmann_k_cgs / UnitMass_in_cgs / UnitLength_in_cgs ** 2 * (UnitTime_in_cgs ** 2)
+print(const_proton_mass)
+print(const_boltzmann_k)
+print()
 
 
 @njit
 def calculate_T(u):
     T_over_mu = (
-        hydro_gamma_minus_one * u * const_proton_mass_cgs / const_boltzmann_k_cgs
+            hydro_gamma_minus_one * u * const_proton_mass / const_boltzmann_k
     )
     if T_over_mu > (T_transition + 1) / mu_ionised:
         return T_over_mu / mu_ionised
@@ -118,4 +128,9 @@ def add_temperature_column():
 
 if __name__ == "__main__":
     # fix_initial_conditions()
-    add_temperature_column()
+    # add_temperature_column()
+    internal_energies = [6.3726251e+02, 7.7903375e+02, 1.7425287e+04, 6.4113910e+04, 3.8831848e+04,
+                         1.1073163e+03, 7.7394878e+03, 7.5230023e+04, 9.1036992e+04, 2.4060946e+00]
+
+    for u in internal_energies:
+        print(calculate_T(u))

slices.py

@@ -1,3 +1,4 @@
+import random
 from pathlib import Path
 from typing import List
 
@@ -7,25 +8,55 @@ import numpy as np
 from matplotlib.colors import LogNorm
 from scipy.interpolate import griddata
 
+from fix_hdf5_masses import calculate_T
 from utils import create_figure
 
 
 def create_2d_slice(
-    input_file: Path, center: List[float], property: str, axis="Z", thickness=3
+        input_file: Path, center: List[float], property: str, axis="Z", thickness=3, method="nearest"
 ):
     axis_names = ["X", "Y", "Z"]
     cut_axis = axis_names.index(axis)
+    limits = {
+        "X": (46, 52),
+        "Y": (54, 60),
+        "Z": (center[cut_axis] - 10, center[cut_axis] + 10)
+    }
     with h5py.File(input_file) as f:
         pt0 = f["PartType0"]
-        coords = pt0["Coordinates"]
+        coords = pt0["Coordinates"][:]
-        data = pt0[property]
+        energies = pt0["InternalEnergies"][:]
+        entropies = pt0["Entropies"][:]
+
         print((center[cut_axis] - thickness < coords[::, cut_axis]).shape)
         # in_slice = (center[cut_axis] - thickness < coords[::, cut_axis]) & (
         #         coords[::, cut_axis] < center[cut_axis] + thickness)
         # print("got slice")
         # coords_in_slice = coords[in_slice]
         # data_in_slice = data[in_slice]
-        print("stats")
+        filter = (
+                (limits["X"][0] < coords[::, 0]) &
+                (coords[::, 0] < limits["X"][1]) &
+
+                (limits["Y"][0] < coords[::, 1]) &
+                (coords[::, 1] < limits["Y"][1]) &
+
+                (limits["Z"][0] < coords[::, 2]) &
+                (coords[::, 2] < limits["Z"][1])
+        )
+        print("before", coords.shape)
+        energies = energies[filter]
+        entropies = entropies[filter]
+        coords = coords[filter]
+
+        print("after", coords.shape)
+        print("calculating temperatures")
+        print(np.random.choice(energies,10))
+        temperatures = np.array([calculate_T(u) for u in energies])
+        print(temperatures.min(),temperatures.max(),temperatures.mean())
+        print("done")
+        exit()
+
         other_axis = {"X": ("Y", "Z"), "Y": ("X", "Z"), "Z": ("X", "Y")}
         x_axis_label, y_axis_label = other_axis[axis]
         x_axis = axis_names.index(x_axis_label)
@@ -34,7 +65,7 @@ def create_2d_slice(
         yrange = np.linspace(coords[::, y_axis].min(), coords[::, y_axis].max(), 1000)
         gx, gy, gz = np.meshgrid(xrange, yrange, center[cut_axis])
         print("interpolating")
-        grid = griddata(coords, data, (gx, gy, gz), method="linear")[::, ::, 0]
+        grid = griddata(coords, temperatures, (gx, gy, gz), method=method)[::, ::, 0]
         print(grid.shape)
         # stats, x_edge, y_edge, _ = binned_statistic_2d(
         #     coords_in_slice[::, x_axis],
@@ -47,15 +78,18 @@ def create_2d_slice(
         # stats = np.nan_to_num(stats)
         print("plotting")
         img = ax.imshow(
-            grid.T,
+            grid,
             norm=LogNorm(),
             interpolation="nearest",
+            origin="lower",
             extent=[xrange[0], xrange[-1], yrange[0], yrange[-1]],
         )
         ax.set_title(input_file.parent.stem)
         ax.set_xlabel(x_axis_label)
         ax.set_ylabel(y_axis_label)
-        fig.colorbar(img, label=property)
+        ax.set_aspect("equal")
+        fig.colorbar(img, label="Temperatures")
         fig.tight_layout()
+        fig.savefig(Path("~/tmp/slice.png").expanduser(), dpi=300)
         plt.show()
         exit()