fix ramses a and cache centering results

auriga_comparison.py

@@ -114,8 +114,10 @@ for dir in sorted(root_dir.glob("*")):
         center = np.array([60.7, 29, 64]) / h
         softening_length = None
     elif "ramses" in dir.name:
+        h = 0.6777
         hr_coordinates, particles_meta, center = load_ramses_data(dir / "output_00007")
         df = pd.DataFrame(hr_coordinates, columns=["X", "Y", "Z"])
+        center = center
         softening_length = None
     else:
         df, particles_meta = read_file(input_file)

find_center.py

@@ -1,12 +1,26 @@
+import hashlib
+import json
+
 import numpy as np
 import pandas as pd
 
 from utils import print_progress
 
+cache_file = "center_cache.json"
+
+try:
+    with open(cache_file, "r") as f:
+        center_cache = json.load(f)
+except FileNotFoundError:
+    center_cache = {}
+
 
 def find_center(df: pd.DataFrame, center: np.ndarray, initial_radius=1):
     # plt.figure()
     all_particles = df[["X", "Y", "Z"]].to_numpy()
+    hash = hashlib.sha256(np.ascontiguousarray(all_particles).data).hexdigest()
+    if hash in center_cache:
+        return np.array(center_cache[hash])
     radius = initial_radius
     center_history = []
     i = 0
@@ -30,4 +44,7 @@ def find_center(df: pd.DataFrame, center: np.ndarray, initial_radius=1):
     # plt.colorbar(label="step")
     # plt.show()
     print()
+    center_cache[hash] = center.tolist()
+    with open(cache_file, "w") as f:
+        json.dump(center_cache, f)
     return center

ramses.py

@@ -12,14 +12,15 @@ def load_ramses_data(ramses_dir: Path):
     s: RamsesSnap = pynbody.load(str(ramses_dir))
     mass_array: SimArray = s.dm["mass"]
     coord_array: SimArray = s.dm["pos"]
-    print("RAMSES a", s.properties["a"])
+    a = s.properties["a"]
+    print("RAMSES a", a)
 
     masses = np.asarray(mass_array.in_units("1e10 Msol"))
     high_res_mass = np.amin(np.unique(masses))  # get lowest mass of particles
     is_high_res_particle = masses == high_res_mass
 
     coordinates = np.asarray(coord_array.in_units("Mpc"))
-    hr_coordinates = coordinates[is_high_res_particle]
+    hr_coordinates = coordinates[is_high_res_particle] / a
 
     particles_meta = ParticlesMeta(particle_mass=high_res_mass)
     center = np.median(hr_coordinates, axis=0)