import yt

if __name__ == "__main__":
    ds = yt.load("/disk12/legacy/GVD_C700_l1600n2048_SLEGAC/dm_gadget/covering_grids/snapshot_101_covering_grid.h5")

    c = ds.arr([0.70095276, 0.84325974, 0.81088409], "unitary")
    w = ds.quan(83, "Mpccm/h")
    box = ds.box(c-w/2, c+w/2)

    V = box.quantities.total_quantity(("index", "cell_volume"))
    M = box.quantities.total_quantity(("grid", "nbody_mass"))
    rho_ave = ds.cosmology.critical_density(0) * (1 + ds.current_redshift)**3 * ds.omega_matter

    delta = (M/V) / rho_ave - 1
    print (f"{delta=}")