Merge pull request #287 from sibirrer/epl_theta_E

Fixed multi-source plane conventions and solver and added EPL velocity dispersions

requirements.txt

@@ -1,5 +1,5 @@
 skypy
-lenstronomy>=1.12.1
+lenstronomy>=1.12.2
 astropy>=5.2
 numpy
 scipy<=1.13.1

slsim/Deflectors/DeflectorTypes/epl_sersic.py

@@ -1,5 +1,6 @@
 from slsim.Deflectors.DeflectorTypes.deflector_base import DeflectorBase
 from slsim.Util.param_util import ellipticity_slsim_to_lenstronomy
+from slsim.Deflectors.velocity_dispersion import theta_E_from_vel_disp_epl
 
 
 class EPLSersic(DeflectorBase):
@@ -16,6 +17,18 @@ class EPLSersic(DeflectorBase):
     - 'z': redshift of deflector
     """
 
+    def __init__(self, deflector_dict):
+        """
+
+        :param deflector_dict: dictionary of deflector quantities
+        """
+        super().__init__(deflector_dict=deflector_dict)
+        try:
+            sis_convention = deflector_dict["sis_convention"]
+        except KeyError:
+            sis_convention = True
+        self._sis_convention = sis_convention
+
     def velocity_dispersion(self, cosmo=None):
         """Velocity dispersion of deflector.
 
@@ -45,13 +58,24 @@ def mass_model_lenstronomy(self, lens_cosmo):
         :type lens_cosmo: ~lenstronomy.Cosmo.LensCosmo instance
         :return: lens_mass_model_list, kwargs_lens_mass
         """
+        gamma = self.halo_properties
         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))
+            lens_light_model_list, kwargs_lens_light = self.light_model_lenstronomy()
+            theta_E = theta_E_from_vel_disp_epl(
+                vel_disp=float(
+                    self.velocity_dispersion(cosmo=lens_cosmo.background.cosmo)
+                ),
+                gamma=gamma,
+                r_half=self.angular_size_light,
+                kwargs_light=kwargs_lens_light,
+                light_model_list=lens_light_model_list,
+                lens_cosmo=lens_cosmo,
+                kappa_ext=0,
+                sis_convention=self._sis_convention,
             )
-        gamma = self.halo_properties
+
         e1_mass, e2_mass = self.mass_ellipticity
         e1_mass_lenstronomy, e2_mass_lenstronomy = ellipticity_slsim_to_lenstronomy(
             e1_slsim=e1_mass, e2_slsim=e2_mass
@@ -113,8 +137,8 @@ def halo_properties(self):
 
         :return: gamma (with =2 is isothermal)
         """
-
         try:
             return float(self._deflector_dict["gamma_pl"])
         except KeyError:
+            # TODO: this can (optionally) be made a function of stellar mass, velocity dispersion etc
             return 2

slsim/Deflectors/deflector.py

@@ -36,7 +36,7 @@ def redshift(self):
 
         :return: redshift
         """
-        return self._deflector.redshift
+        return float(self._deflector.redshift)
 
     def velocity_dispersion(self, cosmo=None):
         """Velocity dispersion of deflector.

slsim/Deflectors/velocity_dispersion.py

@@ -3,6 +3,10 @@
 from scipy import interpolate
 import copy
 
+from lenstronomy.Cosmo.lens_cosmo import LensCosmo
+from lenstronomy.Analysis.light_profile import LightProfileAnalysis
+from lenstronomy.LightModel.light_model import LightModel
+
 """
 This module provides functions to compute velocity dispersion using schechter function.
 """
@@ -35,8 +39,6 @@ def vel_disp_composite_model(r, m_star, rs_star, m_halo, c_halo, cosmo, z_lens):
     }
 
     # turn physical masses to lenstronomy units
-    from lenstronomy.Cosmo.lens_cosmo import LensCosmo
-
     lens_cosmo = LensCosmo(z_lens=z_lens, z_source=10, cosmo=cosmo)
     # Hernquist profile
     sigma0, rs_angle_hernquist = lens_cosmo.hernquist_phys2angular(
@@ -71,6 +73,113 @@ def vel_disp_composite_model(r, m_star, rs_star, m_halo, c_halo, cosmo, z_lens):
     return vel_disp
 
 
+def vel_disp_power_law(
+    theta_E, gamma, r_half, kwargs_light, light_model_list, lens_cosmo
+):
+    """Velocity dispersion for a power-law mass density profile.
+
+    :param theta_E: Einstein radius [arc seconds]
+    :param gamma: power-law slope of deflector
+    :param r_half: half light radius of deflector
+    :param kwargs_light: list of dict for light model parameters
+    :param light_model_list: list of light models
+    :param lens_cosmo: ~LensCosmo instance
+    :return: half light radius averaged velocity dispersion
+    """
+    # turn physical masses to lenstronomy units
+
+    kwargs_mass = [
+        {
+            "theta_E": theta_E,
+            "gamma": gamma,
+            "e1": 0,
+            "e2": 0,
+            "center_x": 0,
+            "center_y": 0,
+        },
+    ]
+    kwargs_anisotropy = {"beta": 0}
+    light_model = LightModel(light_model_list=light_model_list)
+    lensLightProfile = LightProfileAnalysis(light_model=light_model)
+
+    (
+        amps,
+        sigmas,
+        center_x,
+        center_y,
+    ) = lensLightProfile.multi_gaussian_decomposition(
+        kwargs_light,
+        r_h=r_half,
+        n_comp=20,
+    )
+    light_profile_list = ["MULTI_GAUSSIAN"]
+    kwargs_model = {
+        "mass_profile_list": ["EPL"],
+        "light_profile_list": light_profile_list,
+        "anisotropy_model": "const",
+    }
+
+    kwargs_light_mge = [{"amp": amps, "sigma": sigmas}]
+
+    from lenstronomy.GalKin.numeric_kinematics import NumericKinematics
+
+    kwargs_numerics = {
+        "interpol_grid_num": 1000,
+        "log_integration": True,
+        "max_integrate": 1000,
+        "min_integrate": 0.0001,
+        "max_light_draw": None,
+        "lum_weight_int_method": True,
+    }
+
+    kwargs_cosmo = {"d_d": lens_cosmo.dd, "d_s": lens_cosmo.ds, "d_ds": lens_cosmo.dds}
+
+    num_kin = NumericKinematics(kwargs_model, kwargs_cosmo, **kwargs_numerics)
+    vel_disp = num_kin.lum_weighted_vel_disp(
+        r_half, kwargs_mass, kwargs_light_mge, kwargs_anisotropy
+    )
+    return vel_disp
+
+
+def theta_E_from_vel_disp_epl(
+    vel_disp,
+    gamma,
+    r_half,
+    kwargs_light,
+    light_model_list,
+    lens_cosmo,
+    kappa_ext=0,
+    sis_convention=True,
+):
+    """Calculates Einstein radius given measured aperture averaged velocity
+    dispersion and given power-law slope.
+
+    :param vel_disp: velocity dispersion measured within an aperture
+        radius [km/s]
+    :param gamma: power-law slope
+    :param r_half: half light radius (aperture radius) [arc seconds]
+    :param kwargs_light: list of dict for light model parameters
+    :param light_model_list: list of light models
+    :param lens_cosmo: ~LensCosmo instance
+    :param kappa_ext: external convergence
+    :param sis_convention: it True, uses velocity dispersion not as
+        measured one but as the SIS equivalent velocity dispersion
+    :return: Einstein radius matching the velocity dispersion and
+        external convergence
+    """
+    if gamma == 2 or sis_convention is True:
+        theta_E = lens_cosmo.sis_sigma_v2theta_E(vel_disp)
+    else:
+        theta_E_0 = 1
+        vel_disp_0 = vel_disp_power_law(
+            theta_E_0, gamma, r_half, kwargs_light, light_model_list, lens_cosmo
+        )
+        # transform theta_E from vel_disp_0 prediction
+        theta_E = theta_E_0 * (vel_disp / vel_disp_0) ** (2 / (gamma - 1))
+    theta_E /= (1 - kappa_ext) ** (1.0 / (gamma - 1))
+    return theta_E
+
+
 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 =
@@ -190,8 +299,6 @@ def vel_disp_nfw(m_halo, c_halo, cosmo, z_lens):
     """Computes vel_disp_nfw_aperture using the characteristic radius rs of the
     NFW as aperture (which is independent of the source redshift).
 
-    :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
@@ -403,7 +510,7 @@ def schechter_vel_disp_redshift(
         bounds.
     sky_area : `~astropy.units.Quantity`
         Sky area over which galaxies are sampled. Must be in units of solid angle.
-    cosmology : `~astropy.cosmology`
+    cosmo : `~astropy.cosmology`
         `astropy.cosmology` object to calculate comoving densities.
     noise : bool, optional
         Poisson-sample the number of galaxies. Default is `True`.

slsim/Sources/source.py

@@ -343,19 +343,19 @@ def kwargs_variability_extracted(self):
     def redshift(self):
         """Returns source redshift."""
 
-        return self.source_dict["z"]
+        return float(self.source_dict["z"])
 
     @property
     def n_sersic(self):
         """Returns sersic index of the source."""
 
-        return self.source_dict["n_sersic"]
+        return float(self.source_dict["n_sersic"])
 
     @property
     def angular_size(self):
         """Returns angular size of the source."""
 
-        return self.source_dict["angular_size"]
+        return float(self.source_dict["angular_size"])
 
     @property
     def ellipticity(self):

slsim/lens.py

@@ -99,8 +99,8 @@ def __init__(
         self._source_type = self.max_redshift_source_class.source_type
         # we conventionally use highest source redshift in the lens cosmo.
         self._lens_cosmo = LensCosmo(
-            z_lens=float(self.deflector.redshift),
-            z_source=float(self.max_redshift_source_class.redshift),
+            z_lens=self.deflector.redshift,
+            z_source=self.max_redshift_source_class.redshift,
             cosmo=self.cosmo,
         )
 
@@ -210,6 +210,7 @@ def _point_source_image_positions(self, source):
         point_source_pos_x, point_source_pos_y = source.point_source_position(
             center_lens=self.deflector_position, draw_area=self.test_area
         )
+
         # uses analytical lens equation solver in case it is supported by lenstronomy for speed-up
         if (
             self._lens_equation_solver == "lenstronomy_analytical"
@@ -417,27 +418,34 @@ def _einstein_radius(self, source):
         """
         if self.deflector.redshift >= source.redshift:
             theta_E = 0
-        elif self.deflector.deflector_type in ["EPL"]:
-            _lens_cosmo = LensCosmo(
-                z_lens=float(self.deflector.redshift),
-                z_source=float(source.redshift),
-                cosmo=self.cosmo,
-            )
-            kappa_ext = self.los_class.convergence
+            return theta_E
+        lens_model_list, kwargs_lens = self.deflector_mass_model_lenstronomy()
+        lens_model = LensModel(
+            lens_model_list=lens_model_list,
+            z_lens=self.deflector_redshift,
+            z_source_convention=self.max_redshift_source_class.redshift,
+            multi_plane=False,
+            z_source=source.redshift,
+        )
+        if self.deflector.deflector_type in ["EPL"]:
+            kappa_ext_convention = self.los_class.convergence
             gamma_pl = self.deflector.halo_properties
-            theta_E = _lens_cosmo.sis_sigma_v2theta_E(
-                float(self.deflector.velocity_dispersion(cosmo=self.cosmo))
-            ) / (1 - kappa_ext) ** (1.0 / (gamma_pl - 1))
+            theta_E_convention = kwargs_lens[0]["theta_E"]
+            if source.redshift == self.max_redshift_source_class.redshift:
+                theta_E = theta_E_convention
+                kappa_ext = kappa_ext_convention
+            else:
+                beta = self._lens_cosmo.beta_double_source_plane(
+                    z_lens=self.deflector_redshift,
+                    z_source_2=self.max_redshift_source_class.redshift,
+                    z_source_1=source.redshift,
+                )
+                theta_E = theta_E_convention * beta ** (1.0 / (gamma_pl - 1))
+                kappa_ext = kappa_ext_convention * beta
+
+            theta_E /= (1 - kappa_ext) ** (1.0 / (gamma_pl - 1))
         else:
             # numerical solution for the Einstein radius
-            lens_model_list, kwargs_lens = self.deflector_mass_model_lenstronomy()
-            lens_model = LensModel(
-                lens_model_list=lens_model_list,
-                z_lens=self.deflector_redshift,
-                z_source_convention=self.max_redshift_source_class.redshift,
-                multi_plane=False,
-                z_source=source.redshift,
-            )
             lens_analysis = LensProfileAnalysis(lens_model=lens_model)
             theta_E = lens_analysis.effective_einstein_radius(
                 kwargs_lens, r_min=1e-3, r_max=5e1, num_points=50
@@ -817,6 +825,8 @@ def deflector_mass_model_lenstronomy(self):
 
         :return: lens_model_list, kwargs_lens
         """
+        if hasattr(self, "_lens_mass_model_list") and hasattr(self, "_kwargs_lens"):
+            return self._lens_mass_model_list, self._kwargs_lens
         if self.deflector.deflector_type in ["EPL", "NFW_HERNQUIST", "NFW_CLUSTER"]:
             lens_mass_model_list, kwargs_lens = self.deflector.mass_model_lenstronomy(
                 lens_cosmo=self._lens_cosmo
@@ -827,7 +837,7 @@ def deflector_mass_model_lenstronomy(self):
                 % self.deflector.deflector_type
             )
         # adding line-of-sight structure
-        gamma1, gamma2, kappa_ext = self.los_linear_distortions
+        kappa_ext, gamma1, gamma2 = self.los_linear_distortions
         gamma1_lenstronomy, gamma2_lenstronomy = ellipticity_slsim_to_lenstronomy(
             e1_slsim=gamma1, e2_slsim=gamma2
         )
@@ -842,6 +852,8 @@ def deflector_mass_model_lenstronomy(self):
         kwargs_lens.append({"kappa": kappa_ext, "ra_0": 0, "dec_0": 0})
         lens_mass_model_list.append("SHEAR")
         lens_mass_model_list.append("CONVERGENCE")
+        self._kwargs_lens = kwargs_lens
+        self._lens_mass_model_list = lens_mass_model_list
 
         return lens_mass_model_list, kwargs_lens