def METHOD_NAME(r, sigma0, r_core):
    """Rho(r) =  2/pi * Sigma_crit R_c**3 * (R_c**2 + r**2)**(-2)
    :param r: radius (angular scale)
    :param sigma0: convergence in the core
    :param r_core: core radius
    :return: density at radius r
    """
    return 2 / np.pi * sigma0 * r_core**3 * (r_core**2 + r**2) ** (-2)