Merge pull request #60 from aymgal/pr-convergence-profile

Add new mass profile: convergence sheet

coolest/api/profiles/mass.py

@@ -5,7 +5,8 @@
 from scipy import special
 
 from coolest.template.classes.profiles.mass import (PEMD as TemplatePEMD,
-                                                    ExternalShear as TemplateExternalShear)
+                                                    ExternalShear as TemplateExternalShear,
+                                                    ConvergenceSheet as TemplateConvergenceSheet)
 from coolest.api.profiles import util
 
 
@@ -151,11 +152,13 @@ def hessian(self, x, y, theta_E=1., gamma=2., phi=0., q=1., center_x=0., center_
 
 
 class ExternalShear(BaseMassProfile):
+    """
+    Coordinates of the origin for the external shear profile are assumed to be (0., 0.).
+    """ 
 
     _template_class = TemplateExternalShear()
 
     def deflection(self, x, y, gamma_ext=0., phi_ext=0.):
-        """coordinates of the origin for the external shear profile assumed to be (0., 0.)""" 
         phi_ext_ = util.eastofnorth2normalradians(phi_ext)
         gamma1 = gamma_ext * np.cos(2.*phi_ext_)
         gamma2 = gamma_ext * np.sin(2.*phi_ext_)
@@ -178,3 +181,35 @@ def hessian(self, x, y, gamma_ext=0., phi_ext=0.):
         H_xy = gamma2
         H_yx = H_xy
         return H_xx, H_xy, H_yx, H_yy
+
+
+class ConvergenceSheet(BaseMassProfile):
+    """
+    Coordinates of the origin for the convergence sheet are assumed to be (0., 0.).
+    """
+
+    _template_class = TemplateConvergenceSheet()
+
+    def potential(self, x, y, kappa_s=0.):
+        x_ = x # no shift
+        y_ = y # no shift
+        r_ = np.hypot(x_, y_)
+        return 0.5 * kappa_s * r_**2
+
+    def deflection(self, x, y, kappa_s=0.):
+        x_ = x # no shift
+        y_ = y # no shift
+        return x_ * kappa_s, y_ * kappa_s
+
+    def convergence(self, x, y, kappa_s=0.):
+        return np.full_like(x, kappa_s)
+
+    def hessian(self, x, y, kappa_s=0.):
+        kappa = np.full_like(x, kappa_s)
+        gamma1 = 0.
+        gamma2 = 0.
+        H_xx = kappa + gamma1
+        H_yy = kappa - gamma1
+        H_xy = gamma2
+        H_yx = H_xy
+        return H_xx, H_xy, H_yx, H_yy

coolest/template/classes/profiles/mass.py

@@ -15,6 +15,7 @@
     'NFW',
     'Chameleon',
     'ExternalShear', 
+    'ConvergenceSheet', 
     'PixelatedRegularGridPotential',
 ]
 SUPPORTED_CHOICES = __all__
@@ -252,6 +253,7 @@ def __init__(self):
 
 class ExternalShear(AnalyticalProfile):
     """External shear defined with a strength and orientation.
+    The 'origin' of the external shear is by convention fixed to coordinates (0, 0).
     
     This profile is described by the following parameters:
 
@@ -271,6 +273,26 @@ def __init__(self):
         super().__init__(parameters)
 
 
+class ConvergenceSheet(AnalyticalProfile):
+    """Convergence 'sheet', infinite and uniform mass density profile.
+    The 'origin' of the convergence sheet is by convention fixed to coordinates (0, 0).
+    
+    This profile is described by the following parameters:
+
+    - 'kappa_s': convergence of the uniform mass density sheet
+    """
+
+    def __init__(self):
+        parameters = {
+            'kappa_s': NonLinearParameter(
+                "Convergence value of the uniform mass density sheet",
+                DefinitionRange(min_value=0., max_value=1e8),
+                latex_str=r"$\kappa_{\rm s}$"
+            ),
+        }
+        super().__init__(parameters)
+
+
 class PixelatedRegularGridPotential(Profile):
     """Lens potential defined on a grid of regular pixels.
 

test/template/json_test.py

@@ -49,9 +49,13 @@ def test_dump_and_read(self):
         # Defines the external shear
         ext_shear = MassField('my lovely external shear', lens_1.redshift,
                               mass_model=MassModel('ExternalShear'))
+
+        # Defines some convergence sheet
+        conv_sheet = MassField('my mass sheet', lens_1.redshift,
+                               mass_model=MassModel('ConvergenceSheet'))
 
         # Put them in a list, which will also create unique IDs for each profile
-        entity_list = LensingEntityList(ext_shear, lens_1, source_1, source_2, source_3)
+        entity_list = LensingEntityList(ext_shear, lens_1, source_1, source_2, source_3, conv_sheet)
 
         # Define the origin of the coordinates system
         origin = CoordinatesOrigin('00h11m20.244s', '-08d45m51.48s')  # <- in degrees (2.83435, )