import h5py
import numpy
import astropy.cosmology

h = 0.7
omega_m = 0.3
omega_l = 0.7

c = astropy.cosmology.FlatLambdaCDM(100 * h, omega_l)
rho0 = c.critical_density0.to('M_sun/kpc**3').value * omega_m / h**2
with h5py.File('fdm_profiles_p8640_b10000_m7e-23.hdf5', 'r') as f:
	p = f['rho_ms_kpc3'][...]
	p[numpy.isnan(p)] = 0
	r = f['r_bins_kpc'][...]

for halo in [0,1,36,62,75]:
	spherical = numpy.cumsum(p[halo] * 4/3 * numpy.pi * (r[1:]**3 - r[:-1]**3)) / (
		4/3 * numpy.pi * r[1:]**3
	)
	overdense_mask = spherical > 200 * rho0
	overdense_bins, = numpy.nonzero(overdense_mask)
	contiguous = numpy.arange(len(overdense_bins))
	overdense_bins = overdense_bins[overdense_bins == contiguous]
	i = overdense_bins[-1]
	mass = spherical[i] * 4/3 * numpy.pi * r[i + 1]**3
	print(f'{halo}: {mass / h:e}')