Merge remote-tracking branch 'origin/epl_theta_E' into epl_theta_E

LSST-strong-lensing · Nov 22, 2024 · c2e6b01 · c2e6b01
2 parents e8f5b99 + 897f3f8
commit c2e6b01
Show file tree

Hide file tree

Showing 5 changed files with 108 additions and 41 deletions.
diff --git a/slsim/Deflectors/DeflectorTypes/epl_sersic.py b/slsim/Deflectors/DeflectorTypes/epl_sersic.py
@@ -16,6 +16,7 @@ class EPLSersic(DeflectorBase):
     - 'e2_light': eccentricity of light
     - 'z': redshift of deflector
     """
+
     def __init__(self, deflector_dict):
         """
 
@@ -62,12 +63,18 @@ def mass_model_lenstronomy(self, lens_cosmo):
             theta_E = 0.0
         else:
             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)
+            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,
+            )
 
         e1_mass, e2_mass = self.mass_ellipticity
         e1_mass_lenstronomy, e2_mass_lenstronomy = ellipticity_slsim_to_lenstronomy(
@@ -130,9 +137,9 @@ def halo_properties(self):
 
         :return: gamma (with =2 is isothermal)
         """
-        #if hasattr(self._deflector_dict, "gamma_pl"):
+        # if hasattr(self._deflector_dict, "gamma_pl"):
         #    return float(self._deflector_dict["gamma_pl"])
-        #else:
+        # else:
         #    # TODO: this can (optionally) be made a function of stellar mass, velocity dispersion etc
         #    return 2
         try:

diff --git a/slsim/Deflectors/velocity_dispersion.py b/slsim/Deflectors/velocity_dispersion.py
@@ -73,9 +73,10 @@ 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
+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
@@ -88,7 +89,14 @@ def vel_disp_power_law(theta_E, gamma, r_half, kwargs_light, light_model_list, l
     # turn physical masses to lenstronomy units
 
     kwargs_mass = [
-        {"theta_E": theta_E, "gamma": gamma, "e1": 0, "e2": 0, "center_x": 0, "center_y": 0},
+        {
+            "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)
@@ -133,27 +141,39 @@ def vel_disp_power_law(theta_E, gamma, r_half, kwargs_light, light_model_list, l
     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
+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 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
+    :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)
+        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))

diff --git a/slsim/lens.py b/slsim/lens.py
@@ -420,11 +420,13 @@ def _einstein_radius(self, source):
             theta_E = 0
             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)
+        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
@@ -434,9 +436,11 @@ def _einstein_radius(self, source):
                 kappa_ext = kappa_ext_convention
             else:
                 # TODO: translate Einstein radius to different source redshift with power-law slope difference
-                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)
+                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
 

diff --git a/tests/test_Deflectors/test_velocity_dispersion.py b/tests/test_Deflectors/test_velocity_dispersion.py
@@ -5,7 +5,7 @@
     vel_disp_nfw_aperture,
     redshifts_from_comoving_density,
     vel_disp_power_law,
-    theta_E_from_vel_disp_epl
+    theta_E_from_vel_disp_epl,
 )
 from lenstronomy.Cosmo.lens_cosmo import LensCosmo
 import numpy as np
@@ -53,14 +53,26 @@ def test_vel_disp_power_law():
     theta_E = 1.5
     gamma = 2
     r_half = 1
-    kwargs_light = [{"amp": 1, "R_sersic": 2, "n_sersic": 1, "e1": 0, "e2": 0, "center_x": 0, "center_y": 0}]
+    kwargs_light = [
+        {
+            "amp": 1,
+            "R_sersic": 2,
+            "n_sersic": 1,
+            "e1": 0,
+            "e2": 0,
+            "center_x": 0,
+            "center_y": 0,
+        }
+    ]
     light_model_list = ["SERSIC_ELLIPSE"]
 
     # kwargs_light = [{"amp": 1, "gamma": 2, "e1": 0, "e2": 0, "center_x": 0, "center_y": 0}]
     # light_model_list = ["POWER_LAW"]
     lens_cosmo = LensCosmo(z_lens=z_lens, z_source=z_source, cosmo=cosmo)
 
-    vel_disp_pl = vel_disp_power_law(theta_E, gamma, r_half, kwargs_light, light_model_list, lens_cosmo)
+    vel_disp_pl = vel_disp_power_law(
+        theta_E, gamma, r_half, kwargs_light, light_model_list, lens_cosmo
+    )
 
     vel_disp_sis = lens_cosmo.sis_theta_E2sigma_v(theta_E)
     npt.assert_almost_equal(vel_disp_pl / vel_disp_sis, 1, decimal=1)
@@ -74,14 +86,34 @@ def test_theta_E_from_vel_disp_epl():
     theta_E = 1.5
     gamma = 2.5
     r_half = 1
-    kwargs_light = [{"amp": 1, "R_sersic": 2, "n_sersic": 1, "e1": 0, "e2": 0, "center_x": 0, "center_y": 0}]
+    kwargs_light = [
+        {
+            "amp": 1,
+            "R_sersic": 2,
+            "n_sersic": 1,
+            "e1": 0,
+            "e2": 0,
+            "center_x": 0,
+            "center_y": 0,
+        }
+    ]
     light_model_list = ["SERSIC_ELLIPSE"]
 
-    vel_disp_pl = vel_disp_power_law(theta_E, gamma, r_half, kwargs_light, light_model_list, lens_cosmo)
+    vel_disp_pl = vel_disp_power_law(
+        theta_E, gamma, r_half, kwargs_light, light_model_list, lens_cosmo
+    )
 
-    theta_E_out = theta_E_from_vel_disp_epl(vel_disp_pl, gamma, r_half, kwargs_light,
-                                            light_model_list, lens_cosmo, kappa_ext=0, sis_convention=False)
-    npt.assert_almost_equal(theta_E_out/theta_E, 1, decimal=3)
+    theta_E_out = theta_E_from_vel_disp_epl(
+        vel_disp_pl,
+        gamma,
+        r_half,
+        kwargs_light,
+        light_model_list,
+        lens_cosmo,
+        kappa_ext=0,
+        sis_convention=False,
+    )
+    npt.assert_almost_equal(theta_E_out / theta_E, 1, decimal=3)
 
 
 def test_schechter_vdf():

diff --git a/tests/test_lens.py b/tests/test_lens.py
@@ -579,14 +579,14 @@ def setup_method(self):
             deflector_class=self.deflector,
             source_class=[self.source1, self.source2],
             cosmo=self.cosmo,
-            lens_equation_solver="lenstronomy_general"
+            lens_equation_solver="lenstronomy_general",
         )
 
         self.lens_class3_analytical = Lens(
             deflector_class=self.deflector,
             source_class=[self.source1, self.source2],
             cosmo=self.cosmo,
-            lens_equation_solver="lenstronomy_analytical"
+            lens_equation_solver="lenstronomy_analytical",
         )
 
     def test_point_source_arrival_time_multi(self):
@@ -637,9 +637,13 @@ def test_image_observer_time_multi(self):
             observation_time
         )
         # Test multisource image observation time
-        npt.assert_almost_equal(image_observation_time1[0],  image_observation_time3[0][0], decimal=5)
+        npt.assert_almost_equal(
+            image_observation_time1[0], image_observation_time3[0][0], decimal=5
+        )
         # assert image_observation_time1[0] == image_observation_time3[0][0]
-        npt.assert_almost_equal(image_observation_time2, image_observation_time3[1], decimal=5)
+        npt.assert_almost_equal(
+            image_observation_time2, image_observation_time3[1], decimal=5
+        )
         # assert np.all(image_observation_time2 == image_observation_time3[1])
         assert len(self.lens_class3.image_observer_times(t_obs=10)) == 2