Deflector mass method

slsim/Deflectors/DeflectorTypes/deflector_base.py

@@ -89,6 +89,16 @@
         )
         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.

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.
 

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.
 

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.
 

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.
 

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"

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()

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

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