Could you tell me how to compute the r_200 and the r_500 of the subhalo? In general, are there formulas to compute those radii which the snapshot do not provide?
These are provided only for (FoF) halos, i.e. "groups", and not for subhalos.
In general, if you want these values for a central subhalo, then you just look up the value of its parent group (SubhaloGrNr).
If you want these values for a satellite subhalo, it is less clear that this makes physical sense, and these values are not in the catalog, but you could compute them yourself. Then I would suggest to make a cumulative enclosed mass radial profile, going outwards until the enclosed density equals the chosen overdensity, i.e. just follow the definition exactly.
When you calculate the r200 of a halo, you consider all the particles within the sphere or only the particles belong to the FoF?
The values in the catalog (such as Group_R_Crit200) are the usual "spherical overdensity" definitions, i.e. they are computed considering all particles/cells in the simulation.
Thank you. I think I can have a good sleep.
I compared the cirtical density and the values calculated from Group_R_Crit200 and Group_M_Crit200. They are not consistent. So, what happen? The following is my code.
import illustris_python as il
import numpy as np
from astropy.cosmology import FlatLambdaCDM
from astropy import units as u
from astropy.constants import k_B, m_p, M_sun, G
with h5py.File(il.snapshot.snapPath(basePath, snapNum), 'r') as f:
header = dict(list(f['Header'].attrs.items()))
h = header['HubbleParam']
Om0 = header['Omega0']
Ob0 = header['OmegaBaryon']
OmL = header['OmegaLambda']
z = header['Redshift']
BoxSize = header['BoxSize']
H = h*100.
H0 = 100.*h*u.km/(u.s*u.Mpc)
cosmo = FlatLambdaCDM(H0=H,Om0=Om0,Ob0=Ob0)
group_fields = il.groupcat.loadSingle(basePath=basePath,snapNum=snapNum,haloID=haloID)
R_Crit200 = group_fields['Group_R_Crit200']
Group_M_Crit200 = group_fields['Group_M_Crit200']
R_200 = R_Crit200/h
M_200 = Group_M_Crit200*1.e10/h
rho_c = (3.*((cosmo.H(z))**2.)/(8.*np.pi*G)).to_value('g/cm3')
Rho = M_200*M_sun.to_value('g')/(4/3*np.pi*((R_200*u.kpc.to('cm'))**3))/200
The output is
Because it seems correct at z=0 but not at high redshift, the issue is usually a missing scalefactor a to account for comoving units.
In this case, Group_R_Crit200 is comoving (as all all lengths).