Моя рассчитанная доля плотности HI из субобъема Галформа на несколько порядков слишком мала. Что мне не хватает? ⇐ Python

[Anonymous]

Я работаю с данными о галактиках из модели Galform, запущенной в рамках симуляции P-Millennium. Файл, который я использую, содержит данные, составляющие 1/1000 от полного поля моделирования, и я пытаюсь вычислить поле плотности HI и среднее значение \Omega_HI. Вот фрагмент моего кода:

Код: Выделить всё

#Open the file
with h5.File("file.hdf5", "r") as gals:
group = gals["group"]
x = group["xgal"][:]  #galaxy positions in Mpc/h
y = group["ygal"][:]
z = group["zgal"][:]
mHI = group["mcold_atom"][:] # Msun/h

#Shift positions so min=0
x -= np.min(x)
y -= np.min(y)
z -= np.min(z)
galaxy_pos = np.column_stack((x, y, z)) #in Mpc/h

#Define grid parameters
box_dim = np.array([542.16, 542.16, 542.16])  # Mpc/h
N_grids = np.array([150, 150, 150]) #number of pixels
V_cell = np.prod(box_dim/N_grids)

#Project HI onto grid
mHI_per_cell = project_particle_to_regular_grid(
particle_pos=galaxy_pos,
box_len=box_dim,
box_ndim=N_grids,
particle_mass=mHI,
average=False
)[0] #in Msun/h

rhoHI_cell = mHI_per_cell/V_cell #HI mass density per cell in h^2 Msun/Mpc^3
rho_crit = 2.775e11  #crtical density of the Universe today in h^2 Msun/Mpc^3
OmegaHI_cell = rhoHI_cell/rho_crit

print("average OmegaHI =", OmegaHI_cell.mean())

Проблема в том, что мой средний показатель \Omega_HI составляет около 3e-7, что на несколько порядков меньше ожидаемого (~3-5e-4 для z = 0).
Сначала я подозревал, что единицы измерения или координатная сетка, но после проверки:

Код: Выделить всё

print("HI mass in catalogue =", np.sum(mHI))
print("HI mass on grid =", np.sum(mHI_per_cell))
print("Ratio =", np.sum(mHI_per_cell)/np.sum(mHI))

масса идеально сохраняется (отношение ~ 1), поэтому при построении сетки ничего не теряется.
Что мне не хватает? Пожалуйста, помогите.

Подробнее здесь: https://stackoverflow.com/questions/798 ... of-magnitu