Merge branch 'main' of github.com:Findus23/halo_comparison

comparison_plot.py

@@ -181,7 +181,7 @@ def plot_comparison_hist(ax: Axes, file: Path, property: str, m_min=None, m_max=
     if property == "match":
         labels = {
             (-inf, 30): "$M<30$",
-            (None, None): "$M$",
+            (None, None): "all $M$ $[10^{10}$ $\mathrm{M}_\odot]$",
             (30, 100): "$30<M<100$",
             (100, inf): "$100<M$",
         }

readfiles.py

@@ -73,3 +73,24 @@ def read_halo_file(file: Path) -> DataFrame:
 
 def read_fof_file(path: Path):
     file = path / ""
+
+def read_g4_file(file: Path, zoom_type: str) -> Tuple[np.ndarray, np.ndarray]:
+    reference_file = h5py.File(file)
+
+    hubble_param = reference_file["Parameters"].attrs["HubbleParam"]
+    if zoom_type == 'pbh':
+        highres_parttype = 'PartType0'
+        lowres_parttype = 'PartType1'
+    elif zoom_type == 'cdm':
+        highres_parttype = 'PartType1'
+        lowres_parttype = 'PartType2'
+        highres_mass = reference_file['Header'].attrs['MassTable'][1] / hubble_param
+    else:
+        raise ValueError('Please select pbh or cdm as zoom_type!')
+
+    highres_coordinates = reference_file[highres_parttype]["Coordinates"][:] #all coordinates in Mpc/h without adaption!
+    lowres_coordinates = reference_file[lowres_parttype]["Coordinates"][:]
+
+    reference_file.close()
+
+    return highres_coordinates, lowres_coordinates

threed_pbh_zoom.py

@@ -0,0 +1,53 @@
+import sys
+from pathlib import Path
+import numpy as np
+
+import pyvista
+import vtk
+from pandas import DataFrame
+from pyvista import Plotter
+
+from paths import cdm_zoom_dir, pbh_zoom_dir
+from readfiles import read_g4_file
+
+pyvista.set_plot_theme("dark")  # "paraview"
+
+def plotnp3d(pl: Plotter, data: np.ndarray, color="white"):
+    pdata = pyvista.PointSet(data[::, 0:3])
+    # pdata.plot(point_size=1, scalars=data[::, -1], render_points_as_spheres=False, parallel_projection=True,
+    #            anti_aliasing=True, opacity=0.2)
+    # pl.add_points(a) would be equivalent to pl.add_mesh(style="points")
+    pl.add_mesh(
+        pdata,
+        point_size=1,
+        style="points",
+        opacity=0.1,
+        color=color,
+        # scalars=data[::, -1],
+    )
+    # pl.camera_position
+    # pl.enable_stereo_render()
+    # pl.ren_win.SetStereoTypeToSplitViewportHorizontal()
+    # pl.show_grid()
+    pl.enable_terrain_style()
+    pl.enable_parallel_projection()
+    pl.camera.zoom(2)
+    renderer = pl.renderer
+    basic_passes = vtk.vtkRenderStepsPass()
+    blur_pass = vtk.vtkGaussianBlurPass()
+    blur_pass.SetDelegatePass(basic_passes)
+
+    glrenderer = vtk.vtkOpenGLRenderer.SafeDownCast(renderer)
+    glrenderer.SetPass(blur_pass)
+
+if __name__ == "__main__":
+    input_file = cdm_zoom_dir / f"snapshot_{int(sys.argv[1]):03d}.hdf5"
+    highres_coords, lowres_coords = read_g4_file(input_file, 'cdm')
+
+    random_highres_selection = highres_coords[np.random.choice(highres_coords.shape[0], 100000, replace=False)]
+    random_lowres_selection = lowres_coords[np.random.choice(lowres_coords.shape[0], 100000, replace=False)]
+
+    pl = Plotter()
+    plotnp3d(pl, random_lowres_selection)
+    plotnp3d(pl, random_highres_selection)
+    pl.show()