diff --git a/slsim/Deflectors/DeflectorTypes/deflector_base.py b/slsim/Deflectors/DeflectorTypes/deflector_base.py index 01e6a4314..299badc38 100644 --- a/slsim/Deflectors/DeflectorTypes/deflector_base.py +++ b/slsim/Deflectors/DeflectorTypes/deflector_base.py @@ -89,6 +89,16 @@ def mass_ellipticity(self): ) return e1_mass, e2_mass + @abstractmethod + def mass_model_lenstronomy(self, lens_cosmo): + """Returns lens model instance and parameters in lenstronomy conventions. + + :param lens_cosmo: lens cosmology model + :type lens_cosmo: ~lenstronomy.Cosmo.LensCosmo instance + :return: lens_mass_model_list, kwargs_lens_mass + """ + pass + @abstractmethod def light_model_lenstronomy(self, band=None): """Returns lens model instance and parameters in lenstronomy conventions. diff --git a/slsim/Deflectors/DeflectorTypes/epl_sersic.py b/slsim/Deflectors/DeflectorTypes/epl_sersic.py index d82d2784e..d24e702cf 100644 --- a/slsim/Deflectors/DeflectorTypes/epl_sersic.py +++ b/slsim/Deflectors/DeflectorTypes/epl_sersic.py @@ -37,6 +37,38 @@ def light_ellipticity(self): ) return e1_light, e2_light + def mass_model_lenstronomy(self, lens_cosmo): + """Returns lens model instance and parameters in lenstronomy conventions. + + :param lens_cosmo: lens cosmology model + :type lens_cosmo: ~lenstronomy.Cosmo.LensCosmo instance + :return: lens_mass_model_list, kwargs_lens_mass + """ + if lens_cosmo.z_lens >= lens_cosmo.z_source: + theta_E = 0.0 + else: + theta_E = lens_cosmo.sis_sigma_v2theta_E( + float(self.velocity_dispersion(cosmo=lens_cosmo.background.cosmo)) + ) + gamma = self.halo_properties + e1_mass, e2_mass = self.mass_ellipticity + kwargs_lens_mass = [ + { + "theta_E": theta_E, + "gamma": gamma, + "e1": e1_mass, + "e2": e2_mass, + "center_x": self.deflector_center[0], + "center_y": self.deflector_center[1], + } + ] + if gamma == 2: + lens_mass_model_list = ["SIE"] + kwargs_lens_mass[0].pop("gamma") + else: + lens_mass_model_list = ["EPL"] + return lens_mass_model_list, kwargs_lens_mass + def light_model_lenstronomy(self, band=None): """Returns lens model instance and parameters in lenstronomy conventions. diff --git a/slsim/Deflectors/DeflectorTypes/nfw_cluster.py b/slsim/Deflectors/DeflectorTypes/nfw_cluster.py index c1e3e03c2..09962bf08 100644 --- a/slsim/Deflectors/DeflectorTypes/nfw_cluster.py +++ b/slsim/Deflectors/DeflectorTypes/nfw_cluster.py @@ -48,6 +48,49 @@ def velocity_dispersion(self, cosmo=None): ) return vel_disp + def mass_model_lenstronomy(self, lens_cosmo): + """Returns lens model instance and parameters in lenstronomy conventions. + + :param lens_cosmo: lens cosmology model + :type lens_cosmo: ~lenstronomy.Cosmo.LensCosmo instance + :return: lens_mass_model_list, kwargs_lens_mass + """ + lens_mass_model_list, kwargs_lens_mass = self._halo_mass_model_lenstronomy( + lens_cosmo=lens_cosmo + ) + for subhalo in self._subhalos: + lens_mass_model_list_i, kwargs_lens_mass_i = subhalo.mass_model_lenstronomy( + lens_cosmo=lens_cosmo + ) + lens_mass_model_list += lens_mass_model_list_i + kwargs_lens_mass += kwargs_lens_mass_i + return lens_mass_model_list, kwargs_lens_mass + + def _halo_mass_model_lenstronomy(self, lens_cosmo): + """Returns lens model instance and parameters in lenstronomy conventions for the + main halo. + + :param lens_cosmo: lens cosmology model + :type lens_cosmo: ~lenstronomy.Cosmo.LensCosmo instance + :return: lens_mass_model_list, kwargs_lens_mass + """ + lens_mass_model_list = ["NFW_ELLIPSE_CSE"] + e1_mass, e2_mass = self.mass_ellipticity + center_lens = self.deflector_center + m_halo, c_halo = self.halo_properties + rs_halo, alpha_rs = lens_cosmo.nfw_physical2angle(M=m_halo, c=c_halo) + kwargs_lens_mass = [ + { + "alpha_Rs": alpha_rs, + "Rs": rs_halo, + "e1": e1_mass, + "e2": e2_mass, + "center_x": center_lens[0], + "center_y": center_lens[1], + }, + ] + return lens_mass_model_list, kwargs_lens_mass + def light_model_lenstronomy(self, band=None): """Returns lens model instance and parameters in lenstronomy conventions. diff --git a/slsim/Deflectors/DeflectorTypes/nfw_hernquist.py b/slsim/Deflectors/DeflectorTypes/nfw_hernquist.py index 06be8e714..d79b420e7 100644 --- a/slsim/Deflectors/DeflectorTypes/nfw_hernquist.py +++ b/slsim/Deflectors/DeflectorTypes/nfw_hernquist.py @@ -44,6 +44,43 @@ def velocity_dispersion(self, cosmo=None): ) return vel_disp + def mass_model_lenstronomy(self, lens_cosmo): + """Returns lens model instance and parameters in lenstronomy conventions. + + :param lens_cosmo: lens cosmology model + :type lens_cosmo: ~lenstronomy.Cosmo.LensCosmo instance + :return: lens_mass_model_list, kwargs_lens_mass + """ + lens_mass_model_list = ["NFW_ELLIPSE_CSE", "HERNQUIST_ELLIPSE_CSE"] + e1_light_lens, e2_light_lens = self.light_ellipticity + e1_mass, e2_mass = self.mass_ellipticity + rs_phys = lens_cosmo.dd * self.angular_size_light + sigma0, rs_light_angle = lens_cosmo.hernquist_phys2angular( + mass=self.stellar_mass, rs=rs_phys + ) + # halo mass, concentration, stellar mass + m_halo, c_halo = self.halo_properties + rs_halo, alpha_rs = lens_cosmo.nfw_physical2angle(M=m_halo, c=c_halo) + kwargs_lens_mass = [ + { + "alpha_Rs": alpha_rs, + "Rs": rs_halo, + "e1": e1_mass, + "e2": e2_mass, + "center_x": self.deflector_center[0], + "center_y": self.deflector_center[1], + }, + { + "Rs": rs_light_angle, + "sigma0": sigma0, + "e1": e1_light_lens, + "e2": e2_light_lens, + "center_x": self.deflector_center[0], + "center_y": self.deflector_center[1], + }, + ] + return lens_mass_model_list, kwargs_lens_mass + def light_model_lenstronomy(self, band=None): """Returns lens model instance and parameters in lenstronomy conventions. diff --git a/slsim/Deflectors/deflector.py b/slsim/Deflectors/deflector.py index b44de55e0..88199344e 100644 --- a/slsim/Deflectors/deflector.py +++ b/slsim/Deflectors/deflector.py @@ -99,6 +99,15 @@ def mass_ellipticity(self): """ return self._deflector.mass_ellipticity + def mass_model_lenstronomy(self, lens_cosmo): + """Returns lens model instance and parameters in lenstronomy conventions. + + :param lens_cosmo: lens cosmology model + :type lens_cosmo: ~lenstronomy.Cosmo.LensCosmo instance + :return: lens_mass_model_list, kwargs_lens_mass + """ + return self._deflector.mass_model_lenstronomy(lens_cosmo=lens_cosmo) + def light_model_lenstronomy(self, band=None): """Returns lens model instance and parameters in lenstronomy conventions. diff --git a/slsim/lens.py b/slsim/lens.py index 706a95804..1499231d1 100644 --- a/slsim/lens.py +++ b/slsim/lens.py @@ -567,61 +567,10 @@ def deflector_mass_model_lenstronomy(self): :return: lens_model_list, kwargs_lens """ - if self.deflector.deflector_type in ["EPL"]: - gamma = self.deflector.halo_properties - theta_E = self.einstein_radius_deflector - e1_light_lens, e2_light_lens, e1_mass, e2_mass = ( - self.deflector_ellipticity() + if self.deflector.deflector_type in ["EPL", "NFW_HERNQUIST"]: + lens_mass_model_list, kwargs_lens = self.deflector.mass_model_lenstronomy( + lens_cosmo=self._lens_cosmo ) - center_lens = self.deflector_position - - kwargs_lens = [ - { - "theta_E": theta_E, - "gamma": gamma, - "e1": e1_mass, - "e2": e2_mass, - "center_x": center_lens[0], - "center_y": center_lens[1], - } - ] - if gamma == 2: - lens_mass_model_list = ["SIE"] - kwargs_lens[0].pop("gamma") - else: - lens_mass_model_list = ["EPL"] - - elif self.deflector.deflector_type in ["NFW_HERNQUIST"]: - lens_mass_model_list = ["NFW_ELLIPSE_CSE", "HERNQUIST_ELLIPSE_CSE"] - e1_light_lens, e2_light_lens, e1_mass, e2_mass = ( - self.deflector_ellipticity() - ) - center_lens = self.deflector_position - rs_phys = self._lens_cosmo.dd * self.deflector.angular_size_light - sigma0, rs_light_angle = self._lens_cosmo.hernquist_phys2angular( - mass=self.deflector.stellar_mass, rs=rs_phys - ) - # halo mass, concentration, stellar mass - m_halo, c_halo = self.deflector.halo_properties - rs_halo, alpha_rs = self._lens_cosmo.nfw_physical2angle(M=m_halo, c=c_halo) - kwargs_lens = [ - { - "alpha_Rs": alpha_rs, - "Rs": rs_halo, - "e1": e1_mass, - "e2": e2_mass, - "center_x": center_lens[0], - "center_y": center_lens[1], - }, - { - "Rs": rs_light_angle, - "sigma0": sigma0, - "e1": e1_light_lens, - "e2": e2_light_lens, - "center_x": center_lens[0], - "center_y": center_lens[1], - }, - ] else: raise ValueError( "Deflector model %s not supported for lenstronomy model" diff --git a/tests/test_Deflectors/test_DeflectorTypes/test_epl_sersic.py b/tests/test_Deflectors/test_DeflectorTypes/test_epl_sersic.py index b640cff76..21e7e2983 100644 --- a/tests/test_Deflectors/test_DeflectorTypes/test_epl_sersic.py +++ b/tests/test_Deflectors/test_DeflectorTypes/test_epl_sersic.py @@ -1,4 +1,8 @@ +import pytest + from slsim.Deflectors.DeflectorTypes.epl_sersic import EPLSersic +from astropy.cosmology import FlatLambdaCDM +from lenstronomy.Cosmo.lens_cosmo import LensCosmo class TestEPLSersic(object): @@ -34,6 +38,61 @@ def test_velocity_dispersion(self): vel_disp = self.epl_sersic.velocity_dispersion() assert vel_disp == self.deflector_dict["vel_disp"] + def test_mass_model_lenstronomy_sie(self): + # Should yeld SIE model as gamma = 2 + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + lens_cosmo = LensCosmo( + cosmo=cosmo, z_lens=self.epl_sersic.redshift, z_source=2.0 + ) + lens_mass_model_list, kwargs_lens_mass = self.epl_sersic.mass_model_lenstronomy( + lens_cosmo=lens_cosmo + ) + assert len(lens_mass_model_list) == 1 + assert lens_mass_model_list[0] == "SIE" + + def test_mass_model_no_lensing(self): + # case when z_source < z_lens + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + lens_cosmo = LensCosmo( + cosmo=cosmo, z_lens=self.epl_sersic.redshift, z_source=0.2 + ) + lens_mass_model_list, kwargs_lens_mass = self.epl_sersic.mass_model_lenstronomy( + lens_cosmo=lens_cosmo + ) + assert kwargs_lens_mass[0]["theta_E"] == 0.0 + def test_halo_porperties(self): gamma = self.epl_sersic.halo_properties - assert gamma == 2 + assert gamma == 2.0 + + +@pytest.fixture +def gamma_epl_sersic_instance(): + deflector_dict = { + "vel_disp": 200, + "gamma_pl": 1.9, + "e1_mass": 0.1, + "e2_mass": -0.1, + "angular_size": 0.001, + "n_sersic": 1, + "e1_light": -0.1, + "e2_light": 0.1, + "z": 0.5, + } + return EPLSersic(deflector_dict=deflector_dict) + + +def test_mass_model_lenstronomy_gamma(gamma_epl_sersic_instance): + # case when gamma != 2 + gamma_epl_sersic = gamma_epl_sersic_instance + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + lens_cosmo = LensCosmo(cosmo=cosmo, z_lens=gamma_epl_sersic.redshift, z_source=2.0) + lens_mass_model_list, kwargs_lens_mass = gamma_epl_sersic.mass_model_lenstronomy( + lens_cosmo=lens_cosmo + ) + assert len(lens_mass_model_list) == 1 + assert lens_mass_model_list[0] == "EPL" + + +if __name__ == "__main__": + pytest.main() diff --git a/tests/test_Deflectors/test_DeflectorTypes/test_nfw_cluster.py b/tests/test_Deflectors/test_DeflectorTypes/test_nfw_cluster.py index ff0e23a11..1da9d4e1b 100644 --- a/tests/test_Deflectors/test_DeflectorTypes/test_nfw_cluster.py +++ b/tests/test_Deflectors/test_DeflectorTypes/test_nfw_cluster.py @@ -4,6 +4,7 @@ from astropy.table import Table import os import numpy.testing as npt +from lenstronomy.Cosmo.lens_cosmo import LensCosmo class TestNFWCluster(object): @@ -57,3 +58,11 @@ def test_light_model_lenstronomy(self): # one for each subhalo assert len(lens_light_model_list) == 10 assert len(kwargs_lens_light) == 10 + + def test_mass_model_lenstronomy(self): + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + lens_cosmo = LensCosmo(cosmo=cosmo, z_lens=self.halo_dict["z"], z_source=2.0) + lens_mass_model_list, kwargs_lens_mass = ( + self.nfw_cluster.mass_model_lenstronomy(lens_cosmo=lens_cosmo) + ) + assert len(lens_mass_model_list) == 11 diff --git a/tests/test_Deflectors/test_DeflectorTypes/test_nfw_hernquist.py b/tests/test_Deflectors/test_DeflectorTypes/test_nfw_hernquist.py index 0915c104f..af3ce6b5c 100644 --- a/tests/test_Deflectors/test_DeflectorTypes/test_nfw_hernquist.py +++ b/tests/test_Deflectors/test_DeflectorTypes/test_nfw_hernquist.py @@ -1,6 +1,7 @@ from slsim.Deflectors.DeflectorTypes.nfw_hernquist import NFWHernquist from astropy.cosmology import FlatLambdaCDM import numpy.testing as npt +from lenstronomy.Cosmo.lens_cosmo import LensCosmo class TestNFWHernquist(object): @@ -51,3 +52,13 @@ def test_light_model_lenstronomy(self): self.nfw_hernquist.light_model_lenstronomy(band="g") ) assert len(lens_light_model_list) == 1 + + def test_mass_model_lenstronomy(self): + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + lens_cosmo = LensCosmo( + cosmo=cosmo, z_lens=self.deflector_dict["z"], z_source=2.0 + ) + lens_mass_model_list, kwargs_lens_mass = ( + self.nfw_hernquist.mass_model_lenstronomy(lens_cosmo=lens_cosmo) + ) + assert len(lens_mass_model_list) == 2