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

Question

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.pyplot as plt
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[i]

    # 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()