__author__ = "dgilman", "sibirrer"
import numpy as np
from scipy.integrate import quad
from lenstronomy.LensModel.Profiles.nfw import NFW
from lenstronomy.LensModel.Profiles.base_profile import LensProfileBase
__all__ = ["CNFW"]
[docs]
class CNFW(LensProfileBase):
"""
this class computes the lensing quantities of a cored NFW profile:
rho = rho0 * (r + r_core)^-1 * (r + rs)^-2
alpha_Rs is the normalization equivalent to the deflection angle at rs in the absence of a core
"""
model_name = "CNFW"
_s = 0.001 # numerical limit for minimal radius
param_names = ["Rs", "alpha_Rs", "r_core", "center_x", "center_y"]
lower_limit_default = {
"Rs": 0,
"alpha_Rs": 0,
"r_core": 0,
"center_x": -100,
"center_y": -100,
}
upper_limit_default = {
"Rs": 100,
"alpha_Rs": 10,
"r_core": 100,
"center_x": 100,
"center_y": 100,
}
[docs]
def __init__(self):
""""""
self._nfw = NFW()
super(CNFW, self).__init__()
[docs]
def function(self, x, y, Rs, alpha_Rs, r_core, center_x=0, center_y=0):
"""
:param x: angular position
:param y: angular position
:param Rs: angular turn over point
:param alpha_Rs: deflection at Rs (in the absence of a core
:param r_core: core radius
:param center_x: center of halo
:param center_y: center of halo
:return:
"""
x_ = x - center_x
y_ = y - center_y
r = np.sqrt(x_**2 + y_**2)
r = np.maximum(r, self._s)
rho0 = self._alpha2rho0(alpha_Rs=alpha_Rs, Rs=Rs, r_core=r_core)
if isinstance(r, int) or isinstance(r, float):
return self._num_integral_potential(r, Rs, rho0, r_core)
else:
# TODO: currently the numerical integral is done one by one. More efficient is sorting the radial list and
# then perform one numerical integral reading out to the radial points
f_ = []
for i in range(len(r)):
f_.append(self._num_integral_potential(r[i], Rs, rho0, r_core))
return np.array(f_)
def _num_integral_potential(self, r, Rs, rho0, r_core):
"""
:param r:
:param r_core:
:return:
"""
def _integrand(x):
return self.alpha_r(x, Rs, rho0, r_core)
f_ = quad(_integrand, 0, r)[0]
return f_
[docs]
def derivatives(self, x, y, Rs, alpha_Rs, r_core, center_x=0, center_y=0):
rho0 = self._alpha2rho0(alpha_Rs=alpha_Rs, Rs=Rs, r_core=r_core)
if Rs < 0.0000001:
Rs = 0.0000001
x_ = x - center_x
y_ = y - center_y
R = np.sqrt(x_**2 + y_**2)
R = np.maximum(R, self._s)
f_r = self.alpha_r(R, Rs, rho0, r_core)
f_x = f_r * x_ / R
f_y = f_r * y_ / R
return f_x, f_y
[docs]
def hessian(self, x, y, Rs, alpha_Rs, r_core, center_x=0, center_y=0):
# raise Exception('Hessian for truncated nfw profile not yet implemented.')
"""Returns Hessian matrix of function d^2f/dx^2, d^f/dy^2, d^2/dxdy."""
rho0 = self._alpha2rho0(alpha_Rs=alpha_Rs, Rs=Rs, r_core=r_core)
if Rs < 0.0001:
Rs = 0.0001
x_ = x - center_x
y_ = y - center_y
R = np.sqrt(x_**2 + y_**2)
kappa = self.density_2d(x_, y_, Rs, rho0, r_core)
gamma1, gamma2 = self.cnfw_gamma(R, Rs, rho0, r_core, x_, y_)
f_xx = kappa + gamma1
f_yy = kappa - gamma1
f_xy = gamma2
return f_xx, f_xy, f_xy, f_yy
[docs]
def density(self, R, Rs, rho0, r_core):
"""Three dimensional truncated NFW profile.
:param R: radius of interest
:type R: float/numpy array
:param Rs: scale radius
:type Rs: float
:param rho0: density normalization (central core density)
:type rho0: float
:return: rho(R) density
"""
M0 = 4 * np.pi * rho0 * Rs**3
return (M0 / 4 / np.pi) * ((r_core + R) * (R + Rs) ** 2) ** -1
[docs]
def density_lens(self, R, Rs, alpha_Rs, r_core):
"""Computes the density at 3d radius r given lens model parameterization.
The integral in the LOS projection of this quantity results in the convergence
quantity.
"""
rho0 = self._alpha2rho0(alpha_Rs=alpha_Rs, Rs=Rs, r_core=r_core)
return self.density(R, Rs, rho0, r_core)
[docs]
def density_2d(self, x, y, Rs, rho0, r_core, center_x=0, center_y=0):
"""Projected two dimenstional NFW profile (kappa*Sigma_crit)
:param x: radius of interest
:type x: float/numpy array
:param Rs: scale radius
:type Rs: float
:param rho0: density normalization (characteristic density)
:type rho0: float
:return: Epsilon(R) projected density at radius R
"""
x_ = x - center_x
y_ = y - center_y
R = np.sqrt(x_**2 + y_**2)
b = r_core * Rs**-1
x = R * Rs**-1
Fx = self._F(x, b)
return 2 * rho0 * Rs * Fx
[docs]
def mass_3d(self, R, Rs, rho0, r_core):
"""Mass enclosed a 3d sphere or radius r.
:param R:
:param Rs:
:param rho0:
:param r_core:
:return:
"""
b = r_core * Rs**-1
x = R * Rs**-1
M_0 = 4 * np.pi * Rs**3 * rho0
return M_0 * (
x * (1 + x) ** -1 * (-1 + b) ** -1
+ (-1 + b) ** -2
* ((2 * b - 1) * np.log(1 / (1 + x)) + b**2 * np.log(x / b + 1))
)
[docs]
def mass_3d_lens(self, R, Rs, alpha_Rs, r_core):
"""Mass enclosed a 3d sphere or radius r given a lens parameterization with
angular units.
:return:
"""
rho0 = self._alpha2rho0(alpha_Rs=alpha_Rs, Rs=Rs, r_core=r_core)
return self.mass_3d(R, Rs, rho0, r_core)
[docs]
def alpha_r(self, R, Rs, rho0, r_core):
"""Deflection angel of NFW profile along the radial direction.
:param R: radius of interest
:type R: float/numpy array
:param Rs: scale radius
:type Rs: float
:return: Epsilon(R) projected density at radius R
"""
# R = np.maximum(R, self._s)
x = R / Rs
x = np.maximum(x, self._s)
b = r_core * Rs**-1
# b = max(b, 0.000001)
gx = self._G(x, b)
a = 4 * rho0 * Rs**2 * gx / x
return a
[docs]
def cnfw_gamma(self, R, Rs, rho0, r_core, ax_x, ax_y):
"""Shear gamma of NFW profile (times Sigma_crit) along the projection to
coordinate 'axis'.
:param R: radius of interest
:type R: float/numpy array
:param Rs: scale radius
:type Rs: float
:param rho0: density normalization (characteristic density)
:type rho0: float
:return: Epsilon(R) projected density at radius R
"""
c = 0.000001
if isinstance(R, int) or isinstance(R, float):
R = max(R, c)
else:
R[R <= c] = c
x = R * Rs**-1
b = r_core * Rs**-1
b = max(b, c)
gx = self._G(x, b)
Fx = self._F(x, b)
a = (
2 * rho0 * Rs * (2 * gx / x**2 - Fx)
) # /x #2*rho0*Rs*(2*gx/x**2 - Fx)*axis/x
return a * (ax_y**2 - ax_x**2) / R**2, -a * 2 * (ax_x * ax_y) / R**2
[docs]
def mass_2d(self, R, Rs, rho0, r_core):
"""Analytic solution of the projection integral (convergence)"""
x = R / Rs
b = r_core / Rs
b = max(b, 0.000001)
gx = self._G(x, b)
m_2d = 4 * np.pi * rho0 * Rs * R**2 * gx / x**2
return m_2d
def _alpha2rho0(self, alpha_Rs, Rs, r_core):
b = r_core * Rs**-1
gx = self._G(1.0, b)
rho0 = alpha_Rs * (4 * Rs**2 * gx) ** -1
return rho0
def _rho2alpha(self, rho0, Rs, r_core):
b = r_core * Rs**-1
gx = self._G(1.0, b)
alpha = 4 * Rs**2 * gx * rho0
return alpha
def _nfw_func(self, x):
"""Classic NFW function in terms of arctanh and arctan :param x: r/Rs
:return:"""
# c = 0.000000001
if isinstance(x, np.ndarray):
# x[np.where(x<c)] = c
nfwvals = np.ones_like(x)
inds1 = np.where(x < 1)
inds2 = np.where(x > 1)
nfwvals[inds1] = (1 - x[inds1] ** 2) ** -0.5 * np.arctanh(
(1 - x[inds1] ** 2) ** 0.5
)
nfwvals[inds2] = (x[inds2] ** 2 - 1) ** -0.5 * np.arctan(
(x[inds2] ** 2 - 1) ** 0.5
)
return nfwvals
elif isinstance(x, float) or isinstance(x, int):
# x = max(x, c)
if x == 1:
return 1
if x < 1:
return (1 - x**2) ** -0.5 * np.arctanh((1 - x**2) ** 0.5)
else:
return (x**2 - 1) ** -0.5 * np.arctan((x**2 - 1) ** 0.5)
def _F(self, X, b, c=0.001):
"""Analytic solution of the projection integral.
:param X: a dimensionless quantity, either r/rs or r/rc
:type X: float >0
"""
if b == 1:
b = 1 + c
prefac = (b - 1) ** -2
if isinstance(X, np.ndarray):
X[np.where(X == 1)] = 1 - c
output = np.empty_like(X)
inds1 = np.where(np.absolute(X - b) < c)
output[inds1] = (
prefac * (-2 - b + (1 + b + b**2) * self._nfw_func(b)) * (1 + b) ** -1
)
inds2 = np.where(np.absolute(X - b) >= c)
output[inds2] = prefac * (
(X[inds2] ** 2 - 1) ** -1
* (1 - b - (1 - b * X[inds2] ** 2) * self._nfw_func(X[inds2]))
- self._nfw_func(X[inds2] * b**-1)
)
else:
if X == 1:
X = 1 - c
if np.absolute(X - b) < c:
output = (
prefac
* (-2 - b + (1 + b + b**2) * self._nfw_func(b))
* (1 + b) ** -1
)
else:
output = prefac * (
(X**2 - 1) ** -1 * (1 - b - (1 - b * X**2) * self._nfw_func(X))
- self._nfw_func(X * b**-1)
)
return output
def _G(self, X, b, c=0.000001):
"""Analytic solution of integral for NFW profile to compute deflection angel and
gamma.
:param X: R/Rs
:type X: float >0
"""
if b == 1:
b = 1 + c
b2 = b**2
x2 = X**2
fac = (1 - b) ** 2
prefac = fac**-1
if isinstance(X, np.ndarray):
output = np.ones_like(X)
inds1 = np.where(np.absolute(X - b) <= c)
inds2 = np.where(np.absolute(X - b) > c)
output[inds1] = prefac * (
2 * (1 - 2 * b + b**3) * self._nfw_func(b)
+ fac * (-1.38692 + np.log(b2))
- b2 * np.log(b2)
)
output[inds2] = prefac * (
fac * np.log(0.25 * x2[inds2])
- b2 * np.log(b2)
+ 2 * (b2 - x2[inds2]) * self._nfw_func(X[inds2] * b**-1)
+ 2 * (1 + b * (x2[inds2] - 2)) * self._nfw_func(X[inds2])
)
return 0.5 * output
else:
if np.absolute(X - b) <= c:
output = prefac * (
2 * (1 - 2 * b + b**3) * self._nfw_func(b)
+ fac * (-1.38692 + np.log(b2))
- b2 * np.log(b2)
)
else:
output = prefac * (
fac * np.log(0.25 * x2)
- b2 * np.log(b2)
+ 2 * (b2 - x2) * self._nfw_func(X * b**-1)
+ 2 * (1 + b * (x2 - 2)) * self._nfw_func(X)
)
return 0.5 * output