Question on NFWPhaesSpace.

[kgb0255]

Hi, I've been using halotools to sample from NFW distribution using NFWPhaseSpace in halotools.empirical_models. Following this link,

from halotools.empirical_models import NFWPhaseSpace

nfw = NFWPhaseSpace()
data = nfw.mc_generate_nfw_phase_space_points(Ngals=1, mass=1e13, conc=10)
data

gives the following (reformatted):

x: -0.03284552276264489
y: -0.01479777085244665
z: 0.025154780583077484
vx: -67.91742242067731
vy: 86.03371314561191	
vz: -109.0031809751228	
radial_position: 0.04393819948445447
radial_velocity: -84.20070041741873			

By definition radial_velocity must be equal to (x*vx + y*vy + z*vz)/radial_position, which is not what the output gives (the spatial unit Mpc/h cancels). Instead,

(x*vx + y*vy + z*vz)/radial_position = -40.608742502827226

I looked at the source code and found that nfw.mc_generate_nfw_phase_space_points uses _relative_positions_and_velocities to compute relative positions and velocities.

This one sometimes gives an wrong relative velocity. For example, if x1=1, x2=2, v1=1, v2=0, the relative position of particle 1 with respect to particle 2 should be -1, and the relative velocity should be 1. Instead,
_relative_positions_and_velocities(1,2,v1=1,v2=0)
gives
(array([-1.]), array([-1.]))

Should _relative_positions_and_velocities be modified to something like

def _relative_positions_and_velocities(x1, x2, period=None, **kwargs):
    s = _sign_pbc(x1, x2, period=period, equality_fill_val=1.0)
    absd = np.abs(x1 - x2)
    if period is None:
        xrel = s * absd
    else:
        xrel = s * np.where(absd > period / 2.0, period - absd, absd)

    try:
        v1 = kwargs["v1"]
        v2 = kwargs["v2"]
        s_v = _sign_pbc(v1, v2, period=period, equality_fill_val=1.0) #added
        return xrel, s_v * np.abs(v1 - v2) #added
    except KeyError:
        return xrel

[aphearin]

Thanks for catching this @kgb0255. Your example does seem to indicate a problem with this function - I'll take a look ASAP and try out your proposed fix when I do.

[aphearin]

Hi @kgb0255 my apologies for letting this issue get forgotten. But I'm now turning attention to a new release of halotools, and so I've been spending some time this morning sorting through the sign convention issue you raised. I think this may just be due to a misunderstanding of the definition of the sign convention adopted by halotools when defining relative velocities.

Consider the example you gave:

For example, if x1=1, x2=2, v1=1, v2=0, the relative position of particle 1 with respect to particle 2 should be -1, and the relative velocity should be 1

In this example, the particles are moving towards each other, and so their relative velocity should be negative, as returned by halotools, not positive. That convention is explained in this line of the docstring. Does this clarify the issue?