import h5py
import numpy as np

# constants and load data
h      = 0.7
rho_c  =  2.77536627e11 * 0.7**2;
Omegam = 0.3
rho_h  = Omegam * rho_c * 200 # at z=3
f      = h5py.File('fdm_profiles_p8640_b10000_m7e-23.hdf5', 'r')
allrho    = f['rho_ms_kpc3']
r      = f['r_bins_kpc'][:]

# change units
allrho   = allrho[:].T * 1e9  * h**2 # in unit M_solar/Mpc**3
rmid     = (r[1:] + r[:-1])/2 / h  / 1e3 # in unit Mpc
dr       = (r[1:] - r[:-1]) / h / 1e3 # in unit Mpc

#filter
select_simon_halo = 0
rho_halo = allrho[:,select_simon_halo]

nanmask  = np.isfinite(rho_halo)
rho_halo = rho_halo[nanmask]
rmid     = rmid[nanmask]
dr       = dr[nanmask]

count = f['count'][select_simon_halo][nanmask]
cell_len_Mpc = 10 / h / 8640
bin_masses = rho_halo * count * cell_len_Mpc**3
sphere_masses = np.cumsum(bin_masses)

mask     = np.where(sphere_masses / (4/3 * np.pi * rmid**3) <= rho_h)[0][0]
rho_halo = rho_halo[0:mask]
rmid     = rmid[0:mask]
dr       = dr[0:mask]

#size     = mask.shape[0]
size     = mask
rho_halo = rho_halo.reshape(size,1)
rmid     = rmid.reshape(size,1)
dr       = dr.reshape(size,1)

Mh        = sphere_masses[mask]
print ('Mh = %e'%(Mh))