__author__ = "sibirrer"
# this file contains a class to convert lensing and physical units
import numpy as np
import lenstronomy.Util.constants as const
from lenstronomy.Cosmo.background import Background
from lenstronomy.Cosmo.nfw_param import NFWParam
from lenstronomy.Cosmo.gnfw_param import GNFWParam
from lenstronomy.LensModel.Profiles.gnfw import GNFW
__all__ = ["LensCosmo"]
[docs]
class LensCosmo(object):
"""Class to manage the physical units and distances present in a single plane lens
with fixed input cosmology."""
[docs]
def __init__(self, z_lens, z_source, cosmo=None):
"""
:param z_lens: redshift of lens
:param z_source: redshift of source
:param cosmo: ~astropy.cosmology instance
"""
self.z_lens = z_lens
self.z_source = z_source
self.background = Background(cosmo=cosmo)
self.nfw_param = NFWParam(cosmo=cosmo)
self.gnfw_param = GNFWParam(cosmo=cosmo)
@property
def h(self):
return self.background.cosmo.H(0).value / 100.0
@property
def dd(self):
"""
:return: angular diameter distance to the deflector [Mpc]
"""
return self.background.d_xy(0, self.z_lens)
@property
def ds(self):
"""
:return: angular diameter distance to the source [Mpc]
"""
return self.background.d_xy(0, self.z_source)
@property
def dds(self):
"""
:return: angular diameter distance from deflector to source [Mpc]
"""
return self.background.d_xy(self.z_lens, self.z_source)
@property
def ddt(self):
"""
:return: time delay distance [Mpc]
"""
return (1 + self.z_lens) * self.dd * self.ds / self.dds
@property
def sigma_crit(self):
"""Returns the critical projected lensing mass density in units of M_sun/Mpc^2.
:return: critical projected lensing mass density
"""
if not hasattr(self, "_sigma_crit_mpc"):
const_SI = const.c**2 / (
4 * np.pi * const.G
) # c^2/(4*pi*G) in units of [kg/m]
conversion = const.Mpc / const.M_sun # converts [kg/m] to [M_sun/Mpc]
factor = const_SI * conversion # c^2/(4*pi*G) in units of [M_sun/Mpc]
self._sigma_crit_mpc = (
self.ds / (self.dd * self.dds) * factor
) # [M_sun/Mpc^2]
return self._sigma_crit_mpc
@property
def sigma_crit_angle(self):
"""Returns the critical surface density in units of M_sun/arcsec^2 (in physical
solar mass units) when provided a physical mass per physical Mpc^2.
:return: critical projected mass density
"""
if not hasattr(self, "_sigma_crit_arcsec"):
const_SI = const.c**2 / (
4 * np.pi * const.G
) # c^2/(4*pi*G) in units of [kg/m]
conversion = const.Mpc / const.M_sun # converts [kg/m] to [M_sun/Mpc]
factor = const_SI * conversion # c^2/(4*pi*G) in units of [M_sun/Mpc]
self._sigma_crit_arcsec = (
self.ds / (self.dd * self.dds) * factor * (self.dd * const.arcsec) ** 2
) # [M_sun/arcsec^2]
return self._sigma_crit_arcsec
[docs]
def phys2arcsec_lens(self, phys):
"""Convert physical Mpc into arc seconds.
:param phys: physical distance [Mpc]
:return: angular diameter [arcsec]
"""
return phys / self.dd / const.arcsec
[docs]
def arcsec2phys_lens(self, arcsec):
"""Convert angular to physical quantities for lens plane.
:param arcsec: angular size at lens plane [arcsec]
:return: physical size at lens plane [Mpc]
"""
return arcsec * const.arcsec * self.dd
[docs]
def arcsec2phys_source(self, arcsec):
"""Convert angular to physical quantities for source plane.
:param arcsec: angular size at source plane [arcsec]
:return: physical size at source plane [Mpc]
"""
return arcsec * const.arcsec * self.ds
[docs]
def kappa2proj_mass(self, kappa):
"""Convert convergence to projected mass M_sun/Mpc^2.
:param kappa: lensing convergence
:return: projected mass [M_sun/Mpc^2]
"""
return kappa * self.sigma_crit
[docs]
def mass_in_theta_E(self, theta_E):
"""Mass within Einstein radius (area * epsilon crit) [M_sun]
:param theta_E: Einstein radius [arcsec]
:return: mass within Einstein radius [M_sun]
"""
mass = self.arcsec2phys_lens(theta_E) ** 2 * np.pi * self.sigma_crit
return mass
[docs]
def mass_in_coin(self, theta_E):
"""
:param theta_E: Einstein radius [arcsec]
:return: mass in coin calculated in mean density of the universe
"""
chi_L = self.background.T_xy(0, self.z_lens)
chi_S = self.background.T_xy(0, self.z_source)
return (
1.0
/ 3
* np.pi
* (chi_L * theta_E * const.arcsec) ** 2
* chi_S
* self.background.rho_crit
) # [M_sun/Mpc**3]
[docs]
def time_delay_units(self, fermat_pot, kappa_ext=0):
"""
:param fermat_pot: in units of arcsec^2 (e.g. Fermat potential)
:param kappa_ext: unit-less external shear not accounted for in the Fermat potential
:return: time delay in days
"""
D_dt = self.ddt * (1.0 - kappa_ext) * const.Mpc # eqn 7 in Suyu et al.
return (
D_dt / const.c * fermat_pot / const.day_s * const.arcsec**2
) # * self.arcsec2phys_lens(1.)**2
[docs]
def time_delay2fermat_pot(self, dt):
"""
:param dt: time delay in units of days
:return: Fermat potential in units arcsec**2 for a given cosmology
"""
D_dt = self.ddt * const.Mpc
return dt * const.c * const.day_s / D_dt / const.arcsec**2
[docs]
def nfw_angle2physical(self, Rs_angle, alpha_Rs):
"""Converts the angular parameters into the physical ones for an NFW profile.
:param alpha_Rs: observed bending angle at the scale radius in units of arcsec
:param Rs_angle: scale radius in units of arcsec
:return: rho0 [Msun/Mpc^3], Rs [Mpc], c, r200 [Mpc], M200 [Msun]
"""
Rs = Rs_angle * const.arcsec * self.dd
theta_scaled = alpha_Rs * self.sigma_crit * self.dd * const.arcsec
rho0 = theta_scaled / (4 * Rs**2 * (1 + np.log(1.0 / 2.0)))
rho0_com = rho0 / self.h**2
c = self.nfw_param.c_rho0(rho0_com, self.z_lens)
r200 = c * Rs
M200 = self.nfw_param.M_r200(r200 * self.h, self.z_lens) / self.h
return rho0, Rs, c, r200, M200
[docs]
def gnfw_angle2physical(self, Rs_angle, alpha_Rs, gamma_in):
"""Converts the angular parameters into the physical ones for a gNFW profile.
:param alpha_Rs: observed bending angle at the scale radius in units of arcsec
:param Rs_angle: scale radius in units of arcsec
:param gamma_in: inner slope of the gNFW profile
:return: rho0 [Msun/Mpc^3], Rs [Mpc], c, r200 [Mpc], M200 [Msun]
"""
if not hasattr(self, "self._gnfw"):
self._gnfw = GNFW()
Rs = Rs_angle * const.arcsec * self.dd
theta_scaled = alpha_Rs * self.sigma_crit * self.dd * const.arcsec
factor = self._gnfw.get_alpha_Rs_for_kappa_s_1(1, gamma_in) / 4.0
rho0 = theta_scaled / (4 * Rs**2 * factor)
rho0_com = rho0 / self.h**2
c = self.gnfw_param.c_rho0(rho0_com, self.z_lens, gamma_in)
r200 = c * Rs
M200 = self.gnfw_param.M_r200(r200 * self.h, self.z_lens) / self.h
return rho0, Rs, c, r200, M200
[docs]
def nfw_physical2angle(self, M, c):
"""Converts the physical mass and concentration parameter of an NFW profile into
the lensing quantities.
:param M: mass enclosed 200 rho_crit in units of M_sun (physical units, meaning
no little h)
:param c: NFW concentration parameter (r200/r_s)
:return: Rs_angle (angle at scale radius) (in units of arcsec), alpha_Rs
(observed bending angle at the scale radius)
"""
rho0, Rs, r200 = self.nfwParam_physical(M, c)
Rs_angle = Rs / self.dd / const.arcsec # Rs in arcsec
alpha_Rs = rho0 * (4 * Rs**2 * (1 + np.log(1.0 / 2.0)))
return Rs_angle, alpha_Rs / self.sigma_crit / self.dd / const.arcsec
[docs]
def gnfw_physical2angle(self, M, c, gamma_in):
"""Converts the physical mass and concentration parameter of a gNFW profile into
the lensing quantities.
:param M: mass enclosed 200 rho_crit in units of M_sun (physical units, meaning
no little h)
:param c: NFW concentration parameter (r200/r_s)
:return: Rs_angle (angle at scale radius) (in units of arcsec), alpha_Rs
(observed bending angle at the scale radius
"""
if not hasattr(self, "self._gnfw"):
self._gnfw = GNFW()
rho0, Rs, r200 = self.gnfwParam_physical(M, c, gamma_in)
Rs_angle = Rs / self.dd / const.arcsec # Rs in arcsec
factor = self._gnfw.get_alpha_Rs_for_kappa_s_1(1, gamma_in) / 4.0
alpha_Rs = rho0 * (4 * Rs**2 * factor)
return Rs_angle, alpha_Rs / self.sigma_crit / self.dd / const.arcsec
[docs]
def nfwParam_physical(self, M, c):
"""Returns the NFW parameters in physical units.
:param M: physical mass in M_sun in definition m200
:param c: concentration
:return: rho0 [Msun/Mpc^3], Rs [Mpc], r200 [Mpc]
"""
r200 = (
self.nfw_param.r200_M(M * self.h, self.z_lens) / self.h
) # physical radius r200
rho0 = (
self.nfw_param.rho0_c(c, self.z_lens) * self.h**2
) # physical density in M_sun/Mpc**3
Rs = r200 / c
return rho0, Rs, r200
[docs]
def gnfwParam_physical(self, M, c, gamma_in):
"""Returns the gNFW parameters in physical units.
:param M: physical mass in M_sun in definition m200
:param c: concentration
:return: rho0 [Msun/Mpc^3], Rs [Mpc], r200 [Mpc]
"""
r200 = (
self.gnfw_param.r200_M(M * self.h, self.z_lens) / self.h
) # physical radius r200
rho0 = (
self.gnfw_param.rho0_c(c, self.z_lens, gamma_in) * self.h**2
) # physical density in M_sun/Mpc**3
Rs = r200 / c
return rho0, Rs, r200
[docs]
def nfw_M_theta_r200(self, M):
"""Returns r200 radius in angular units of arc seconds on the sky.
:param M: physical mass in M_sun
:return: angle (in arc seconds) of the r200 radius
"""
r200 = (
self.nfw_param.r200_M(M * self.h, self.z_lens) / self.h
) # physical radius r200
theta_r200 = r200 / self.dd / const.arcsec
return theta_r200
[docs]
def gnfw_M_theta_r200(self, M):
"""Returns r200 radius in angular units of arc seconds on the sky.
:param M: physical mass in M_sun
:return: angle (in arc seconds) of the r200 radius
"""
r200 = (
self.gnfw_param.r200_M(M * self.h, self.z_lens) / self.h
) # physical radius r200
theta_r200 = r200 / self.dd / const.arcsec
return theta_r200
[docs]
def sis_theta_E2sigma_v(self, theta_E):
"""Converts the lensing Einstein radius into a physical velocity dispersion.
:param theta_E: Einstein radius (in arcsec)
:return: velocity dispersion in units (km/s)
"""
v_sigma_c2 = theta_E * const.arcsec / (4 * np.pi) * self.ds / self.dds
return np.sqrt(v_sigma_c2) * const.c / 1000
[docs]
def sis_sigma_v2theta_E(self, v_sigma):
"""Converts the velocity dispersion into an Einstein radius for a SIS profile.
:param v_sigma: velocity dispersion (km/s)
:return: theta_E (arcsec)
"""
theta_E = (
4
* np.pi
* (v_sigma * 1000.0 / const.c) ** 2
* self.dds
/ self.ds
/ const.arcsec
)
return theta_E
[docs]
def hernquist_phys2angular(self, mass, rs):
"""Translates physical mass definitions of the Hernquist profile to the angular
units used in the Hernquist lens profile of lenstronomy.
'sigma0' is defined such that the deflection at projected RS leads to alpha =
2./3 * Rs * sigma0
:param mass: A spherical overdensity mass in M_sun corresponding to the mass
definition mdef at redshift z
:param rs: rs in units of physical Mpc
:return: sigma0, Rs_angle
"""
rs_angle = rs / self.dd / const.arcsec # Rs in arcsec
rhos = mass / (2 * np.pi) / rs**3 # units of M_sun / Mpc^3
sigma0 = rhos * rs # units of M_sun / Mpc^2
sigma0 /= self.sigma_crit
return sigma0, rs_angle
[docs]
def hernquist_angular2phys(self, sigma0, rs_angle):
"""'sigma0' is defined such that the deflection at projected RS leads to alpha =
2./3 * Rs * sigma0.
:param sigma0: convergence normalization
:param rs_angle: rs in angular units [arcseconds]
:return: mass [M_sun], rs [Mpc]
"""
rs = rs_angle * self.dd * const.arcsec # units of Mpc
rhos = sigma0 / rs * self.sigma_crit
m_tot = 2 * np.pi * rhos * rs**3
return m_tot, rs
[docs]
def uldm_angular2phys(self, kappa_0, theta_c):
"""Converts the anguar parameters entering the LensModel Uldm() (Ultra Light
Dark Matter) class in physical masses, i.e. the total soliton mass and the mass
of the particle.
:param kappa_0: central convergence of profile
:param theta_c: core radius (in arcseconds)
:return: m_eV_log10, M_sol_log10, the log10 of the masses, m in eV and M in
M_sun
"""
D_Lens = self.dd * 10**6 # in parsec
Sigma_c = self.sigma_crit * 10 ** (-12) # in M_sun / parsec^2
r_c = theta_c * const.arcsec * D_Lens
rho0 = 2048 * np.sqrt(0.091) * kappa_0 * Sigma_c / (429 * np.pi * r_c)
m_log10 = -22 + 0.5 * np.log10(190 / rho0 * (r_c / 100) ** (-4))
M_log10 = 9 + np.log10(160 * 1.4 / r_c) - 2 * (m_log10 + 22)
return m_log10, M_log10
[docs]
def uldm_mphys2angular(self, m_log10, M_log10):
"""Converts physical ULDM mass in the ones, in angular units, that enter the
LensModel Uldm() class.
:param m_log10: exponent of ULDM mass in eV
:param M_log10: exponent of soliton mass in M_sun
:return: kappa_0, theta_c, the central convergence and core radius (in
arcseconds)
"""
D_Lens = self.dd * 10**6 # in parsec
Sigma_c = self.sigma_crit * 10 ** (-12) # in M_sun/parsec^2
m22 = 10 ** (m_log10 + 22)
M9 = 10 ** (M_log10 - 9)
r_c = 160 * 1.4 * m22 ** (-2) * M9 ** (-1) # core radius in parsec
rho0 = (
190 * m22 ** (-2) * (r_c / 100) ** (-4)
) # central density in M_sun/parsec^3
kappa_0 = 429 * np.pi * rho0 * r_c / (2048 * np.sqrt(0.091) * Sigma_c)
theta_c = r_c / D_Lens / const.arcsec
return kappa_0, theta_c
[docs]
def sersic_m_star2k_eff(self, m_star, R_sersic, n_sersic):
"""Translates a total stellar mass into 'k_eff', the convergence at 'R_sersic'
(effective radius or half-light radius) for a Sersic profile.
:param m_star: total stellar mass in physical Msun
:param R_sersic: half-light radius in arc seconds
:param n_sersic: Sersic index
:return: k_eff
"""
# compute mass integral
from lenstronomy.LensModel.Profiles.sersic_utils import SersicUtil
sersic_util = SersicUtil()
norm_integral = sersic_util.total_flux(
amp=1, R_sersic=R_sersic, n_sersic=n_sersic
)
# compute total kappa normalization and re
k_eff = m_star / self.sigma_crit_angle
# renormalize
k_eff /= norm_integral
return k_eff
[docs]
def vel_disp_dPIED_sigma0(self, vel_disp, Ra, Rs):
"""Sigma0 value in the convention of the lenstronomy pseudo_jaffe lens model.
lenstronomy conventions:
.. math::
\\Sigma(R) = \\Sigma_0 \\frac{Ra Rs}{Rs-Ra}\\left(\\frac{1}{\\sqrt{Ra^2+R^2}} - \\frac{1}{\\sqrt{Rs^2+R^2}} \\right)
with
.. math::
\\Sigma_0 = \\pi \\rho_0 \\frac{Ra Rs}{Rs + Ra}
In the lensing parameterization,
.. math::
\\sigma_0 = \\frac{\\Sigma_0}{\\Sigma_{\\rm crit}}
SIS profile:
.. math::
\\theta_{\\rm E, SIS} = 4\\pi \\frac{D_{\\rm LS}}{D_{\\rm S}} \\left(\\frac{\\sigma_v}{c} \\right)^2
.. math::
\\kappa_{\\rm SIS} = \\frac{1}{2} \\frac{\\theta_E}{R}
.. math::
\\kappa_{\\rm dPIED} \\approx \\sigma_0 \\frac{R_s R_a}{R_s - R_a} \\left(\\sim \\frac{1}{R} \\right))
relation then between velocity dispersion and sigma0:
.. math::
\\sigma_v = c \\sqrt{\\sigma_0 \\frac{R_s R_a}{R_s - R_a} \\frac{D_{\\rm S}}{D_{\\rm LS}} \\frac{1}{2\\pi} }
:param vel_disp: SIS equivalent velocity dispersion (km/s)
:return: sigma0 value in the convention of the lenstronomy pseudo_jaffe lens model
"""
sigma_0 = (
(vel_disp * 1000 / const.c) ** 2
* 2
* np.pi
* self.dds
/ self.ds
* ((Rs - Ra) / (Rs * Ra))
/ const.arcsec
)
return sigma_0
[docs]
def sersic_k_eff2m_star(self, k_eff, R_sersic, n_sersic):
"""Translates convergence at half-light radius to total integrated physical
stellar mass for a Sersic profile.
:param k_eff: lensing convergence at half-light radius
:param R_sersic: half-light radius in arc seconds
:param n_sersic: Sersic index
:return: stellar mass in physical Msun
"""
from lenstronomy.LensModel.Profiles.sersic_utils import SersicUtil
sersic_util = SersicUtil()
norm_integral = sersic_util.total_flux(
amp=1, R_sersic=R_sersic, n_sersic=n_sersic
)
m_star = k_eff * self.sigma_crit_angle * norm_integral
return m_star
[docs]
def beta_double_source_plane(self, z_lens, z_source_1, z_source_2):
"""Model prediction of ratio of scaled deflection angles.
:param z_lens: lens redshift
:param z_source_1: source_1 redshift
:param z_source_2: source_2 redshift
:param cosmo: ~astropy.cosmology instance
:return: beta
"""
return self.background.beta_double_source_plane(
z_lens=z_lens, z_source_1=z_source_1, z_source_2=z_source_2
)
[docs]
def theta_E_power_law_scaling(
self,
theta_E_convention,
kappa_ext_convention,
gamma_pl,
z_lens,
z_source_convention,
z_source,
):
"""Maps Einstein radius of a power-law profile with external convergence to
different source redshifts.
:param theta_E_convention: Einstein radius for the lens when a source is at
z_source_conventions coming from the main deflector (excluding external
convergence)
:param kappa_ext_convention: external convergence for z_source_convention
:param gamma_pl: power-law slope of the deflector
:param z_lens: lens redshift
:param z_source_convention: source redshift for lens model conventions
:param z_source: source redshift
:return: Einstein radius for a source at redshift z_source
"""
if z_source == z_source_convention:
theta_E = theta_E_convention
kappa_ext = kappa_ext_convention
else:
beta = self.beta_double_source_plane(
z_lens=z_lens,
z_source_2=z_source_convention,
z_source_1=z_source,
)
kappa_ext = kappa_ext_convention * beta
theta_E = theta_E_convention * beta ** (1.0 / (gamma_pl - 1.0))
# theta_E = theta_E_convention * (beta + kappa_ext_convention * (1 - beta)) ** (1. / (gamma_pl - 1.))
# β − (1 − λ)(1 − β) ** (1.0 / (gamma_pl - 1))
theta_E /= (1 - kappa_ext) ** (1.0 / (gamma_pl - 1))
return theta_E