Miscentering integrand

[eduardojsbarroso]

In the theory modeling/miscentering.py, there function:

def integrand_surface_density_nfw(theta, r_proj, r_mis, r_s):
    r"""Computes integrand for surface mass density with the NFW profile.

    Parameters
    ----------
    theta : float
        Angle of polar coordinates of the miscentering direction.
    r_proj : float
        Projected radial position from the cluster center in :math:`M\!pc`.
    r_mis : float
        Projected miscenter distance in :math:`M\!pc`.
    r_s : float
        Scale radius

    Returns
    -------
    float
        2D projected density in units of :math:`M_\odot\ Mpc^{-2}`.
    """
    r_norm = np.sqrt(r_proj**2.0 + r_mis**2.0 - 2.0 * r_proj * r_mis * np.cos(theta)) / r_s

    r2m1 = r_norm**2.0 - 1.0
    if r_norm < 1:
        sqrt_r2m1 = np.sqrt(-r2m1)
        res = np.arcsinh(sqrt_r2m1 / r_norm) / (-r2m1) ** (3.0 / 2.0) + 1.0 / r2m1
    elif r_norm > 1:
        sqrt_r2m1 = np.sqrt(r2m1)
        res = -np.arcsin(sqrt_r2m1 / r_norm) / (r2m1) ** (3.0 / 2.0) + 1.0 / r2m1
    else:
        res = 1.0 / 3.0
    return res

Breaks when theta=0 and r_mis = r_proj. When integrating this function sometimes this error breaks everything (np.arcsinh(sqrt_r2m1 / r_norm) goes to infinity).