I am trying to do calculations to find the subhalo velocity dispersion and check it with the SubhaloVelDisp field, however, the values I get do not match the SubhaloVelDisp values. Specifically, how is the SubhaloVelDisp field being calculated?
There are some subtleties of the units, etc, it is computed as (where the num loop is over the member particles):
disp = 0;
mass = 0;
for(i = 0; i < num; i++)
for(j = 0; j < 3; j++)
ddxx = GRAVITY_NEAREST_X(P[p].Pos[j] - SubhaloCM[j]); // periodic distance
dx[j] = current_scalefactor * ddxx;
dv[j] = vel_to_phys * (P[p].Vel[j] - SubhaloCMVel[j]); // units to km/s
dv[j] += H_of_a * dx[j]; // hubble correction
disp += P[p].Mass * dv[j] * dv[j];
mass += P[p].Mass;
veldisp = sqrt(disp / (3 * mass)); /* convert to 1d velocity dispersion */
Thanks for your response. In the code you sent above, do you know how the GRAVITY_NEAREST_X, vel_to_phys, and H_of_a are defined?
GRAVITY_NEAREST_X is the usual periodic distance function, e.g.
if x > BoxHalf:
if x < -BoxHalf:
vel_to_phys is whatever you need to get that quantity into km/s, e.g. likely 1/a.
H_of_a is the usual H(z) = H0 * sqrt(Ωm(1+z)^3 + ΩΛ + Ωk(1+z)^2).
H(z) = H0 * sqrt(Ωm(1+z)^3 + ΩΛ + Ωk(1+z)^2)
Also, for the P[p].Mass above, I figure that the Mass is a list of the masses of each particle in num. I am trying to recreate that whole code in Python, but I am not sure what mass value I would use in replace of P[p].Mass. Do you know what value/field I should use for trying to recreate that code? Thanks.
It should just be the Masses particle field of the stars.