Skip to content

Commit

Permalink
Merge pull request #148 from furcelay/cluster_lenses
Browse files Browse the repository at this point in the history 
Skeleton for cluster lens simulations
  • Loading branch information
@nkhadka21
nkhadka21 authored May 15, 2024
2 parents 4d7ca84 + 537690b commit 9fb1f8f
Show file tree
Hide file tree
Showing 6 changed files with 268 additions and 0 deletions.
73 changes: 73 additions & 0 deletions slsim/Deflectors/DeflectorTypes/nfw_cluster.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
from slsim.Deflectors.DeflectorTypes.deflector_base import DeflectorBase
from slsim.Deflectors.velocity_dispersion import vel_disp_nfw_aperture
from lenstronomy.Cosmo.lens_cosmo import LensCosmo


class NFWCluster(DeflectorBase):
"""Class of a NFW halo lens model with subhalos. Each subhalo is a Deflector
instance with its own mass and light.
required quantities in dictionary:
- 'halo_mass': halo mass in physical M_sol
- 'concentration': halo concentration
- 'e1_mass': eccentricity of NFW profile
- 'e2_mass': eccentricity of NFW profile
- 'z': redshift of deflector
"""

def __init__(self, deflector_dict, subhalos_list):
"""
:param deflector_dict: parameters of the cluster halo
:type deflector_dict: dict
:param subhalos_list: list with Deflector instances as cluster subhalos
:type subhalos_list: list[Deflector]
"""
super(NFWCluster, self).__init__(deflector_dict)
self._subhalos = subhalos_list

def velocity_dispersion(self, cosmo=None):
"""Velocity dispersion of deflector. Simplified assumptions on anisotropy and
averaged over the characteristic radius.
:param cosmo: cosmology
:type cosmo: ~astropy.cosmology class
:return: velocity dispersion [km/s]
"""
# convert radian to arc seconds
lens_cosmo = LensCosmo(z_lens=self.redshift, z_source=10, cosmo=cosmo)

m_halo, c_halo = self.halo_properties
rs_arcsec, _ = lens_cosmo.nfw_physical2angle(m_halo, c_halo)
vel_disp = vel_disp_nfw_aperture(
r=rs_arcsec,
m_halo=m_halo,
c_halo=c_halo,
cosmo=cosmo,
z_lens=self.redshift,
)
return vel_disp

def light_model_lenstronomy(self, band=None):
"""Returns lens model instance and parameters in lenstronomy conventions.
:param band: imaging band
:type band: str
:return: lens_light_model_list, kwargs_lens_light
"""
lens_light_model_list, kwargs_lens_light = [], []
for subhalo in self._subhalos:
lens_light_model_list_i, kwargs_lens_light_i = (
subhalo.light_model_lenstronomy(band=band)
)
lens_light_model_list += lens_light_model_list_i
kwargs_lens_light += kwargs_lens_light_i
return lens_light_model_list, kwargs_lens_light

@property
def halo_properties(self):
"""Properties of the NFW halo.
:return: halo mass M200 [physical M_sol], concentration r200/rs
"""
return self._deflector_dict["halo_mass"], self._deflector_dict["concentration"]
112 changes: 112 additions & 0 deletions slsim/Deflectors/velocity_dispersion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,118 @@ def vel_disp_composite_model(r, m_star, rs_star, m_halo, c_halo, cosmo, z_lens):
return vel_disp


def vel_disp_nfw_3d(r, m_halo, c_halo, cosmo, z_lens):
"""Computes the unweighted velocity dispersion at 3D radius r for a deflector with a
NFW halo profile, assuming isotropic anisotropy (beta = 0).
Based on equation (14) of Lokas and Mamon 2001 (
https://arxiv.org/abs/astro-ph/0002395)
:param r: radius of the unweighted velocity dispersion [arcsec]
:param m_halo: Halo mass [physical M_sun]
:param c_halo: halo concentration
:param cosmo: cosmology
:type cosmo: ~astropy.cosmology class
:param z_lens: redshift of the deflector
:return: velocity dispersion [km/s]
"""

from scipy.special import spence
from lenstronomy.Cosmo.lens_cosmo import LensCosmo
from astropy.constants import G
from lenstronomy.Util.constants import arcsec

lens_cosmo = LensCosmo(z_lens=z_lens, z_source=10, cosmo=cosmo)
r_halo = lens_cosmo.nfw_M_theta_r200(m_halo)
s = r / r_halo
vel2 = (
G.to("km2 Mpc / M_sun s2").value * m_halo / (r_halo * arcsec * lens_cosmo.dd)
) # km^2 / s^2
cs = c_halo * s
g_c = 1 / (np.log(1 + c_halo) - c_halo / (1 + c_halo))
vel_disp2 = vel2 * (
1
/ 2
* c_halo**2
* g_c
* s
* (1 + cs) ** 2
* (
np.pi**2
- np.log(cs)
- 1 / cs
- 1 / (1 + cs) ** 2
- 6 / (1 + cs)
+ (1 + 1 / cs**2 - 4 / cs - 2 / (1 + cs)) * np.log(1 + cs)
+ 3 * np.log(1 + cs) ** 2
+ 6 * spence(1 + cs)
)
)
return np.sqrt(vel_disp2)


def vel_disp_nfw_aperture(r, m_halo, c_halo, cosmo, z_lens):
"""Computes the average line-of-sight velocity dispersion in an aperture r for a
deflector with a NFW halo profile, assuming isotropic anisotropy (beta = 0).
Based on equation (48) of Lokas & Mamon 2001 (
https://arxiv.org/abs/astro-ph/0002395)
:param r: radius of the aperture for the velocity dispersion [arcsec]
:param m_halo: Halo mass [physical M_sun]
:param c_halo: halo concentration
:param cosmo: cosmology
:type cosmo: ~astropy.cosmology class
:param z_lens: redshift of the deflector
:return: velocity dispersion [km/s]
"""
from lenstronomy.Cosmo.lens_cosmo import LensCosmo
from lenstronomy.LensModel.Profiles.nfw import NFW

def _log_integrate(func, xmin, xmax, n_grid=200):
min_log = np.log(xmin)
max_log = np.log(xmax)
dlogx = (max_log - min_log) / (n_grid - 1)
x = np.logspace(
min_log + dlogx / 2.0,
max_log + dlogx / 2.0,
n_grid,
base=np.e,
)
dlog_x = np.log(x[2]) - np.log(x[1])
y = func(x)
return np.sum(y * dlog_x * x)

lens_cosmo = LensCosmo(z_lens=z_lens, z_source=10, cosmo=cosmo)

g_c = 1 / (np.log(1 + c_halo) - c_halo / (1 + c_halo))
rs, alpha_rs = lens_cosmo.nfw_physical2angle(m_halo, c_halo)
r_halo = rs * c_halo
rmin = 1e-3 * rs
rmax = 1e3 * rs
nfw = NFW()
m_2d_r = nfw.mass_2d_lens(r, rs, alpha_rs) * lens_cosmo.sigma_crit_angle
int1 = _log_integrate(
lambda r_: vel_disp_nfw_3d(r_, m_halo, c_halo, cosmo, z_lens) ** 2
* r_
/ r_halo
/ (1 + c_halo * r_ / r_halo) ** 2,
rmin,
rmax,
)
int2 = _log_integrate(
lambda r_: (
vel_disp_nfw_3d(r_, m_halo, c_halo, cosmo, z_lens) ** 2
/ (1 + c_halo * r_ / r_halo) ** 2
* np.sqrt((r_ / r_halo) ** 2 - (r / r_halo) ** 2)
),
r,
rmax,
)
vel_disp2 = c_halo**2 * g_c * m_halo / m_2d_r * (int1 - int2) / r_halo
return np.sqrt(vel_disp2)


def vel_disp_sdss(sky_area, redshift, vd_min, vd_max, cosmology, noise=True):
"""Velocity dispersion function in a cone matched by SDSS measurements.
Expand Down
Binary file added tests/TestData/halo_NFW.fits
View file
Binary file not shown.
7 changes: 7 additions & 0 deletions tests/TestData/subhalos_table.fits
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 208 / length of dimension 1 NAXIS2 = 10 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 22 / number of table fields TTYPE1 = 'z ' TFORM1 = 'D ' TTYPE2 = 'M ' TFORM2 = 'D ' TTYPE3 = 'coeff ' TFORM3 = '5D ' TDIM3 = '(5) ' TTYPE4 = 'ellipticity' TFORM4 = 'D ' TTYPE5 = 'physical_size' TFORM5 = 'D ' TUNIT5 = 'kpc ' TTYPE6 = 'stellar_mass' TFORM6 = 'D ' TTYPE7 = 'angular_size' TFORM7 = 'D ' TUNIT7 = 'rad ' TTYPE8 = 'mag_g ' TFORM8 = 'D ' TTYPE9 = 'mag_r ' TFORM9 = 'D ' TTYPE10 = 'mag_i ' TFORM10 = 'D ' TTYPE11 = 'mag_z ' TFORM11 = 'D ' TTYPE12 = 'mag_Y ' TFORM12 = 'D ' TTYPE13 = 'vel_disp' TFORM13 = 'D ' TTYPE14 = 'e1_light' TFORM14 = 'D ' TTYPE15 = 'e2_light' TFORM15 = 'D ' TTYPE16 = 'e1_mass ' TFORM16 = 'D ' TTYPE17 = 'e2_mass ' TFORM17 = 'D ' TTYPE18 = 'n_sersic' TFORM18 = 'D ' TTYPE19 = 'center_x_mass' TFORM19 = 'D ' TTYPE20 = 'center_y_mass' TFORM20 = 'D ' TTYPE21 = 'center_x_light' TFORM21 = 'D ' TTYPE22 = 'center_y_light' TFORM22 = 'D ' END ?�7�P��1� ����?�����>��D6O�?�$L�ֲ?�ZE>�?\‹n6s?Ҡ�=�?�͏>4wA�ڃG�Ó>�H�Q)%=@8��獺@7���j�@7>�@6ąo9�@6��V:p@m���I��/����?�"��𰀿`K�b���?�y4��^a@@����[@[*o3GR@����[@[*o3GR?�Xl9p!$�5U���a�?��<��#>�@gd,&=?�9���?�F�Y��?F�~5�?�ln%��@�$8���BL �r��_>ˉ��rk@5�G��5V@4���@3\��zsN@2�Ks u�@2���]�@q09�\���,�!K�?�l�p�6K���S���?�4��
@?��t4e+?�~RÅ�=?��t4e+?�~RÅ�=?��?�㽢�1�=�!�x?�Ë�?l�J@�?h�B�V�?�dh���?PS����?��x,��C?��r�t^�A�rEn'�>�>�*K�@8����@7,�)&l?@6�,��@6E��L @6ЬA:�@m��"�)?��}�K�c?�x�~/6}?����b�?Σ&@~�@����:�����k�"����:�����k�"?۹Ǐ˽��5��Ԉ�?��fk��#?!G�:.:0?u(y�$�}?���~�B?M�'^�A&?���֕}�@"�B�BNE���{�>�C�3�@@5$��P@3�B�ݩ�@3T���H@3�Y���@2�L�w�H@qN�A7k�?�����?����T�?�O�{Ί?La���=�@?��4J��-@"�F*�.�?��4J��-@"�F*�.�?�x!N �4R�7s�?���=�>�CB�Z�?�v�hK��?�-n�T�?E���s^?ԣ�O>�r?�z�\��B3.��`>���N@6xc儗X@4���%��@4e1Po
@4�P%�@3�`)j@p�{�.N������k�?�����p翓��;d?��T?߱@@zbn��b@�@=@zbn��b@�@=?��g��c�2�)i�_?�/B�5�*>�Ѳ�r��?��I��;�?���HU?[=3�7�?��3�t??�9�a���B�i�3 >�klG�2@87�< ��@6̒���@6%wYw�2@5�@��|@5��p��@n{b9�E��NN��?���&1֍�����Z��?�S��\@��u8U �ts�0���u8U �ts�0�?�=a�*��0��ԭ�?���m��>�/1��?i���7h�?��s�Z?,��}ڣ?�˪Q�.?���=5A�Ύ�{>{�& �q@:e?
�c@8ԎoVE@8.�URi-@7�+5��@7��z �@m���[�v���!5$S?mN 0��?�ѕүj6����� �e@@��９@"���0��@��９@"���0��?�j��,�/qTEv?�Z{�S�~?I��E?r���a�x?�l����*?1�G����?�e<��?��bZ�dA�Lg�mV>��SG߶@;
��)�<@9�oY��p@9��[@8�U����@8{)9�T@m~H6�뿣�vߕ{?��B�B�����
���?�a��9 �@@0�}F��b�;��� @0�}F��b�;��� ?���/�T��3�-/��?� (����? ��3�?�?��@-�?����4t?P�Ex�?��+�,N?�q�V�G�B%��=�>�n��3�@6�m���z@5�=�՞�@5g�@4�!DBֲ@4��;�2j@o6���?��Ks��?�k����?��� u:?ď��-/`@����F��4��@O�����F��4��@O�?ܛ��9���5B�~)�?�L^s� I>���͛0t?QB=}�^?��9*��f?@�i8�?��]�
�@>T�yPBI�<�"��>��U�vod@5�!x{ʻ@4l�3��Q@3�.��x�@3ZL�Z�@3&�H�^W@q86e)?�ȏ�i?��H-��%?��I�ی�?����j=�@��y���Y@ g`���V��y���Y@ g`���V
Expand Down
Loading

0 comments on commit 9fb1f8f

Please sign in to comment.