from ...backend_obj import backend
from ...constants import arcsec_to_rad, G_over_c2, rad_to_arcsec
from ...utils import translate_rotate
[docs]
def concentration_tnfw(mass, Rs, critical_density, d_l, DELTA=200.0):
"""
Compute the concentration parameter "c" for a TNFW profile.
Parameters
----------
mass: ArrayLike
The mass of the lens.
*Unit: Msun*
Rs: ArrayLike
The scale radius of the TNFW lens.
*Unit: arcsec*
critical_density: ArrayLike
The critical density of the universe.
*Unit: Msun / Mpc^3*
d_l: ArrayLike
The angular diameter distance to the lens.
*Unit: Mpc*
DELTA: float
The overdensity parameter.
*Unit: unitless*
Returns
-------
ArrayLike
The concentration parameter "c" for a TNFW profile.
*Unit: unitless*
"""
r_delta = (3 * mass / (4 * backend.pi * DELTA * critical_density)) ** (1 / 3) # fmt: skip
return r_delta / (Rs * d_l * arcsec_to_rad) # fmt: skip
def _F_tnfw(x):
"""
Compute the function F(x) for a TNFW profile.
Helper method from Baltz et al. 2009 equation A.5
"""
x_gt1 = backend.clamp(x, min=1 + 1e-3)
x_lt1 = backend.clamp(x, max=1 - 1e-3)
return backend.where(
x < 1 - 1e-2,
backend.arccosh(1 / x_lt1) / backend.sqrt(1.0 - x_lt1**2),
backend.where(
x > 1 + 1e-2,
backend.arccos(1 / x_gt1) / backend.sqrt(x_gt1**2 - 1.0),
1 - 2 * (x - 1) / 3 + 7 * (x - 1) ** 2 / 15, # where: x == 1
),
)
def _L_tnfw(x, tau):
"""
Helper method from Baltz et al. 2009 equation A.6
"""
return backend.log(x / (tau + backend.sqrt(tau**2 + x**2))) # fmt: skip
[docs]
def scale_density_tnfw(c, critical_density, DELTA=200.0):
return DELTA / 3 * critical_density * c**3 / (backend.log(1 + c) - c / (1 + c)) # fmt: skip
[docs]
def M0_totmass_tnfw(mass, tau):
"""
Rearranged from Baltz et al. 2009 equation A.4
"""
return mass * (tau**2 + 1) ** 2 / (tau**2 * ((tau**2 - 1) * backend.log(tau) + backend.pi * tau - (tau**2 + 1))) # fmt: skip
[docs]
def M0_scalemass_tnfw(Rs, c, critical_density, d_l, DELTA=200.0):
"""What M0 would be for an NFW
Parameters
----------
Rs: ArrayLike
The scale radius of the TNFW lens.
*Unit: arcsec*
c: ArrayLike
The concentration parameter of an NFW lens with the same parameters.
*Unit: unitless*
critical_density: ArrayLike
The critical density of the universe.
*Unit: Msun / Mpc^3*
d_l: ArrayLike
The angular diameter distance to the lens.
*Unit: Mpc*
DELTA: float
The overdensity parameter.
*Unit: unitless*
"""
return 4 * backend.pi * (Rs * d_l * arcsec_to_rad) ** 3 * scale_density_tnfw(c, critical_density, DELTA) # fmt: skip
[docs]
def mass_enclosed_2d_tnfw(
r,
Rs,
tau,
M0,
):
"""
Total projected mass (Msun) within a radius r (arcsec). Given in Baltz et
al. 2009 equation A.11
Parameters
----------
r: ArrayLike
Radius at which to compute the enclosed mass.
*Unit: arcsec*
mass: ArrayLike
Mass of the lens.
*Unit: Msun*
Rs: ArrayLike
Scale radius of the TNFW lens.
*Unit: arcsec*
tau: ArrayLike
Truncation scale. Ratio of truncation radius to scale radius.
*Unit: unitless*
Returns
-------
ArrayLike
Integrated mass projected in infinite cylinder within radius r.
*Unit: Msun*
"""
g = r / Rs
t2 = tau**2
F = _F_tnfw(g)
L = _L_tnfw(g, tau)
a1 = t2 / (t2 + 1) ** 2
a2 = (t2 + 1 + 2 * (g**2 - 1)) * F
a3 = tau * backend.pi
a4 = (t2 - 1) * backend.log(tau)
a5 = backend.sqrt(t2 + g**2) * (-backend.pi + (t2 - 1) * L / tau) # fmt: skip
return M0 * a1 * (a2 + a3 + a4 + a5)
[docs]
def physical_deflection_angle_tnfw(
x0,
y0,
Rs,
tau,
x,
y,
M0,
d_l,
s=0.0,
):
"""
Compute the physical deflection angle for a TNFW profile. Converted from
Baltz et al. 2009 equation A.18
Parameters
----------
x0: ArrayLike
The x-coordinate of the lens center.
*Unit: arcsec*
y0: ArrayLike
The y-coordinate of the lens center.
*Unit: arcsec*
Rs: ArrayLike
The scale radius of the TNFW lens.
*Unit: arcsec*
tau: ArrayLike
The truncation scale. Ratio of truncation radius to scale radius.
*Unit: unitless*
x: ArrayLike
The x-coordinate in the lens plane.
*Unit: arcsec*
y: ArrayLike
The y-coordinate in the lens plane.
*Unit: arcsec*
M0: ArrayLike
The mass normalization constant. See `M0_totmass_tnfw` and
`M0_scalemass_tnfw`.
*Unit: Msun*
d_l: ArrayLike
The angular diameter distance to the lens.
*Unit: Mpc*
s: float
Softening parameter to prevent numerical instabilities.
*Unit: arcsec*
"""
x, y = translate_rotate(x, y, x0, y0)
r = backend.sqrt(x**2 + y**2) + s
theta = backend.arctan2(y, x)
# The below actually equally comes from eq 2.13 in Meneghetti notes
dr = mass_enclosed_2d_tnfw(r, Rs, tau, M0) / (
r * d_l * arcsec_to_rad
) # note dpsi(u)/du = 2x*dpsi(x)/dx when u = x^2 # fmt: skip
S = 4 * G_over_c2 * rad_to_arcsec
return S * dr * backend.cos(theta), S * dr * backend.sin(theta)
[docs]
def convergence_tnfw(
x0,
y0,
Rs,
tau,
x,
y,
critical_density,
M0,
d_l,
s=0.0,
):
"""
Compute the dimensionless convergence for the TNFW. See Baltz et al. 2009
equation A.8
Parameters
----------
x0: ArrayLike
The x-coordinate of the lens center.
*Unit: arcsec*
y0: ArrayLike
The y-coordinate of the lens center.
*Unit: arcsec*
Rs: ArrayLike
The scale radius of the TNFW lens.
*Unit: arcsec*
tau: ArrayLike
The truncation scale. Ratio of truncation radius to scale radius.
*Unit: unitless*
x: ArrayLike
The x-coordinate in the lens plane.
*Unit: arcsec*
y: ArrayLike
The y-coordinate in the lens plane.
*Unit: arcsec*
M0: ArrayLike
The mass normalization constant. See `M0_totmass_tnfw` and
`M0_scalemass_tnfw`.
*Unit: Msun*
d_l: ArrayLike
The angular diameter distance to the lens.
*Unit: Mpc*
s: float
Softening parameter to prevent numerical instabilities.
*Unit: arcsec*
"""
x, y = translate_rotate(x, y, x0, y0)
r = backend.sqrt(x**2 + y**2) + s
g = r / Rs
F = _F_tnfw(g)
L = _L_tnfw(g, tau)
S = M0 / (2 * backend.pi * (Rs * d_l * arcsec_to_rad) ** 2) # fmt: skip
t2 = tau**2
a1 = t2 / (t2 + 1) ** 2
a2 = backend.where(g == 1, (t2 + 1) / 3.0, (t2 + 1) * (1 - F) / (g**2 - 1)) # fmt: skip
a3 = 2 * F
a4 = -backend.pi / backend.sqrt(t2 + g**2)
a5 = (t2 - 1) * L / (tau * backend.sqrt(t2 + g**2))
return a1 * (a2 + a3 + a4 + a5) * S / critical_density # fmt: skip
def _P_tnfw(x):
"""
Compute the function F(x) for a TNFW profile.
Helper method from Baltz et al. 2009 equation A.5
"""
x_gt1 = backend.clamp(x, min=1 + 1e-4)
x_lt1 = backend.clamp(x, max=1 - 1e-4)
return backend.where(
x < 1 - 1e-3,
-backend.arccosh(1 / x_lt1) ** 2,
backend.where(
x > 1 + 1e-3,
backend.arccos(1 / x_gt1) ** 2,
(2 * (x - 1) - 5 * (x - 1) ** 2 / 3), # where: x == 1 # fixme
),
)
[docs]
def potential_tnfw(
x0,
y0,
Rs,
tau,
x,
y,
M0,
d_l,
d_s,
d_ls,
s=0.0,
):
"""
Compute the lensing potential for a TNFW profile. See Baltz et al. 2009
equation A.14
Parameters
----------
x0: ArrayLike
The x-coordinate of the lens center.
*Unit: arcsec*
y0: ArrayLike
The y-coordinate of the lens center.
*Unit: arcsec*
Rs: ArrayLike
The scale radius of the TNFW lens.
*Unit: arcsec*
tau: ArrayLike
The truncation scale. Ratio of truncation radius to scale radius.
*Unit: unitless*
x: ArrayLike
The x-coordinate in the lens plane.
*Unit: arcsec*
y: ArrayLike
The y-coordinate in the lens plane.
*Unit: arcsec*
M0: ArrayLike
The mass normalization constant. See `M0_totmass_tnfw` and
`M0_scalemass_tnfw`.
*Unit: Msun*
d_l: ArrayLike
The angular diameter distance to the lens.
*Unit: Mpc*
d_s: ArrayLike
The angular diameter distance to the source.
*Unit: Mpc*
d_ls: ArrayLike
The angular diameter distance between the lens and the source.
*Unit: Mpc*
s: float
Softening parameter to prevent numerical instabilities.
*Unit: arcsec*
"""
x, y = translate_rotate(x, y, x0, y0)
r = backend.sqrt(x**2 + y**2) + s
g = r / Rs
t2 = tau**2
u = g**2
F = _F_tnfw(g)
L = _L_tnfw(g, tau)
# fmt: off
S = 2 * M0 * G_over_c2 * (d_ls / d_s) / (d_l * arcsec_to_rad**2)
a1 = 1 / (t2 + 1) ** 2
a2 = 2 * backend.pi * t2 * (tau - backend.sqrt(t2 + u) + tau * backend.log(tau + backend.sqrt(t2 + u)))
a3 = 2 * (t2 - 1) * tau * backend.sqrt(t2 + u) * L
a4 = t2 * (t2 - 1) * L**2
a5 = 4 * t2 * (u - 1) * F
a6 = t2 * (t2 - 1) * _P_tnfw(g)
a7 = t2 * ((t2 - 1) * backend.log(tau) - t2 - 1) * backend.log(u)
a8 = t2 * ((t2 - 1) * backend.log(tau) * backend.log(4 * tau) + 2 * backend.log(tau / 2) - 2 * tau * (tau - backend.pi) * backend.log(2 * tau))
return S * a1 * (a2 + a3 + a4 + a5 + a6 + a7 - a8)
