astro_tools.py 文件源码

def nfw(R, norm_constant, Rs, Rcore):
    '''
    2D Navarro-Frenk-White surface radial profile probability density

    See:
    Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
    Bartelmann, M., A&A, 1996, 313, 697
    Rykoff, E.S., et al., ApJ, 746, 178

    @param R: Radius
    @param norm_constant: Normalization constant
    @param Rs: Scale radius
    @param Rcore: Since NFW profile diverges at R=0, the value at the center is held fixed starting at Rcore

    @return probability density of profile at R
    '''

    def compute_nfw(R):
        if R < Rcore:
            R2 = Rcore
        else:
            R2 = R

        if (R2/Rs) < 0.999:
            func = (1 - 2 * math.atanh(math.sqrt((1 - R2/Rs) / (R2/Rs + 1))) / math.sqrt(1 - (R2/Rs)**2)) / ((R2/Rs)**2 - 1)

        # There are some computational issues as R2/Rs -> 1, using taylor expansion of the function
        # around this point
        elif (R2/Rs) < 1.001:
            func = -(20/63)*((R2/Rs)**3 - 1) + (13/35)*((R2/Rs)**2 - 1) - (2/5)*(R2/Rs) + (11/15)

        else:
            func = (1 - 2 * math.atan(math.sqrt((R2/Rs - 1) / (R2/Rs + 1))) / math.sqrt((R2/Rs)**2 - 1)) / ((R2/Rs)**2 - 1)


        return norm_constant * 2 * math.pi * R * func

    if np.isscalar(R):
        return compute_nfw(R)

    else:
        return np.fromiter(map(compute_nfw, R), np.float, count = len(R))