[Python] Im not able to find the correct conversion factor for my code

I have a code for finding out the lymann alpha surface brightness maps. I have used the surface brightness value. However, the slice data value which should have the entries as 3.94x10^62 has the value 5.02x10^62 instead. Please let me know if there is something wrong with either the formula or the code.

import matplotlib.colors as colors
from astropy.cosmology import Planck15 as cosmo
from astropy import units as u


# Load data
with h5py.File(file_path, 'r') as f:
lum_dens_data = f['maps/TNG50_total'][:] # Luminosity density in erg/s/pMpc³
width_arcsec = f.attrs['width_arcsec'] # Angular width in arcseconds
depth_angstrom = f.attrs['depth_angstrom'] # Observed depth in Ångströms

# Constants
z = 2.0 # Redshift
lambda_rest = 1215.67 * u.Angstrom # Rest-frame Lyα wavelength
c = 3e10 * u.cm / u.s # Speed of light in cm/s
H_z = cosmo.H(z).to(u.s**-1).value # Hubble parameter at redshift z in s⁻¹
d_lum = cosmo.luminosity_distance(z).to(u.cm).value # Luminosity distance in cm
d_A = cosmo.angular_diameter_distance(z).to(u.cm).value # Angular diameter distance in cm

# Corrected y(z) factor
y_factor = (c / H_z) / lambda_rest.to(u.cm).value # Dimensionless factor

# Conversion factors
str_to_arcsec2 = (1 * u.sr).to(u.arcsec**2).value # Steradian to arcseconds squared
conversion_factor = (u.Mpc.to(u.cm)) ** -3 # Convert pMpc³ to cm³

# Convert luminosity density from erg/s/pMpc³ to erg/s/cm³
lum_dens_data_cm3 = lum_dens_data * conversion_factor

# Apply scaling factor (1e6)
lum_dens_data_cm3 *= 1e-6 # Scale luminosity density

# Surface brightness calculation using corrected Eq. (10)
surface_brightness = (
lum_dens_data_cm3 / (4 * np.pi * d_lum**2) * d_A**2 * y_factor
) # SB in erg/s/cm²/str


# Convert steradians to arcseconds squared
surface_brightness /= str_to_arcsec2 # SB in erg/s/cm²/arcsec²

# Plotting with fixed range for each slice
for i in range(surface_brightness.shape[0]):
slice_data = surface_brightness

# Fixed color scale limits
vmin = 1e-24 # Minimum surface brightness
vmax = 1e-18 # Maximum surface brightness

plt.figure(figsize=(10, 10))
plt.imshow(slice_data, norm=colors.LogNorm(vmin=vmin, vmax=vmax),
cmap='RdBu', extent=(0, width_arcsec, 0, width_arcsec)) # Use a softer colormap like 'plasma'
plt.colorbar(label=r'Surface Brightness $\gamma$ (erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$)')
plt.title(f"Lyman-alpha Surface Brightness Map (z = 2, Slice {i+1})")
plt.xlabel("Arcsec")
plt.ylabel("Arcsec")
plt.savefig(f'Lyman_alpha_surface_brightness_slice_{i+1}.png', dpi=300)
plt.show()

Continue reading...