Cartesian outputs

Note: The IonEnergy and IonFlux fields are not densities. They are also affected by the reduced speed of light \(\tilde{c} = 0.2\,c\). For example, after converting to cgs units the photo-ionization integral would be \(\Gamma_{\rm HI} = \tilde{c} \sigma_{\rm HI} {\tt IonEnergy} / {\tt MeanPhotonEnergy} / {\tt VoxelVolume}\).

Note: The IonLuminosity, LyaLuminosity, and StarFormationRate fields are not densities.

import numpy as np
import h5py

def read_render(snap=300, ren_dir='Thesan-1/output'):
    # Read header info
    with h5py.File(f'{ren_dir}/cartesian_{snap:03d}/cartesian_{snap:03d}.000.hdf5', 'r') as f:
        g = f['Header'] # Access the HDF5 Header group
        NumFiles = g.attrs['NumFiles'] # Number of subfiles per output
        NumPixels = g.attrs['NumPixels'] # Number of cells in each Cartesian dimension
        NumPixels3 = NumPixels**3 # Total number of Cartesian cells
        BoxSize = g.attrs['BoxSize'] # Spatial extent of the periodic box [ckpc/h]
        HubbleParam = g.attrs['HubbleParam'] # Hubble constant [100 km/s/Mpc]
        Time = g.attrs['Time'] # Scale factor a=1/(1+z) of the current output
        Redshift = g.attrs['Redshift'] # Redshift of the current output
        Omega0 = g.attrs['Omega0'] # Cosmological matter density [rho_crit_0]
        OmegaBaryon = g.attrs['OmegaBaryon'] # Cosmological baryonic matter density [rho_crit_0]
        OmegaLambda = g.attrs['OmegaLambda'] # Cosmological constant density [rho_crit_0]
        UnitLength_in_cm = g.attrs['UnitLength_in_cm'] # Code length [cm]
        UnitMass_in_g = g.attrs['UnitMass_in_g'] # Code mass [g]
        UnitVelocity_in_cm_per_s = g.attrs['UnitVelocity_in_cm_per_s'] # Code velocity [cm/s]

    # Define various convenience factors
    a,z,h = Time, Redshift, HubbleParam
    length_to_cgs = a * UnitLength_in_cm / h
    volume_to_cgs = length_to_cgs**3
    mass_to_cgs = UnitMass_in_g / h
    density_to_cgs = mass_to_cgs / volume_to_cgs
    velocity_to_cgs = np.sqrt(a) * UnitVelocity_in_cm_per_s
    cell_volume = (length_to_cgs * BoxSize / float(NumPixels))**3
    print(f'snap = {snap}  |  BoxSize = {BoxSize} ckpc/h  |  NumPixels = {NumPixels}')

    # Read full data
    offset = 0
    Density = np.zeros(NumPixels3, dtype=np.float32)
    for sub in range(NumFiles):
        filename = f'{ren_dir}/cartesian_{snap:03d}/cartesian_{snap:03d}.{sub:03d}.hdf5'
        with h5py.File(filename, 'r') as f:
            n_buf = len(f['Density'])
            next_offset = offset + n_buf
            Density[offset:next_offset] = f['Density'][:]
            offset = next_offset # Keep track of file offset
    Density = Density.reshape([NumPixels,NumPixels,NumPixels]) # 3D density grid