Merge pull request #266 from Jammy2211/feature/cored_NFW

Feature/cored nfw
the cored spherical NFW deflection angles. Best used with the profile cNFWMCRLudlowSph

Author: NiekWielders

autogalaxy/config/priors/mass/dark/cnfw.yaml

@@ -0,0 +1,51 @@
+cNFWSph:
+  centre_0:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  centre_1:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  kappa_s:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 1.0
+    width_modifier:
+      type: Relative
+      value: 0.2
+    limits:
+      lower: 0.0
+      upper: inf
+  scale_radius:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 30.0
+    width_modifier:
+      type: Relative
+      value: 0.2
+    limits:
+      lower: 0.0
+      upper: inf
+  core_radius:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 15.0
+    width_modifier:
+      type: Relative
+      value: 0.2
+    limits:
+      lower: 0.0
+      upper: inf

autogalaxy/config/priors/mass/dark/cnfw_mcr.yaml

@@ -0,0 +1,61 @@
+cNFWMCRLudlowSph:
+  centre_0:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  centre_1:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  mass_at_200:
+    type: LogUniform
+    lower_limit: 100000000.0
+    upper_limit: 1000000000000000.0
+    width_modifier:
+      type: Relative
+      value: 0.5
+    limits:
+      lower: 0.0
+      upper: inf
+  f_c:
+    type: Uniform
+    lower_limit: 0.0001
+    upper_limit: 0.5
+    width_modifier:
+      type: Relative
+      value: 0.2
+    limits:
+      lower: 0.0001
+      upper: inf
+  redshift_object:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 1.0
+    width_modifier:
+      type: Relative
+      value: 0.5
+    limits:
+      lower: 0.0
+      upper: inf
+  redshift_source:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 1.0
+    width_modifier:
+      type: Relative
+      value: 0.5
+    limits:
+      lower: 0.0
+      upper: inf

autogalaxy/profiles/mass/__init__.py

@@ -36,6 +36,8 @@
     NFWMCRLudlow,
     gNFWMCRLudlow,
     NFWVirialMassConcSph,
+    cNFWSph,
+    cNFWMCRLudlowSph,
 )
 from .stellar import (
     Gaussian,

autogalaxy/profiles/mass/dark/__init__.py

@@ -9,3 +9,5 @@
 from .nfw_truncated_mcr import NFWTruncatedMCRLudlowSph, NFWTruncatedMCRDuffySph
 from .nfw_truncated_mcr_scatter import NFWTruncatedMCRScatterLudlowSph
 from .nfw_virial_mass_conc import NFWVirialMassConcSph
+from .cnfw import cNFWSph
+from .cnfw_mcr import cNFWMCRLudlowSph

autogalaxy/profiles/mass/dark/cnfw.py

@@ -0,0 +1,175 @@
+import numpy as np
+
+from typing import Tuple
+
+from autogalaxy.profiles.mass.dark.abstract import AbstractgNFW
+
+import autoarray as aa
+
+class cNFWSph(AbstractgNFW):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        kappa_s: float = 0.05,
+        scale_radius: float = 1.0,
+        core_radius: float = 0.01,
+    ):
+        """
+                Represents a spherical cored NFW density distribution
+
+                Parameters
+                ----------
+                centre
+                    The (y,x) arc-second coordinates of the profile centre.
+                kappa_s
+                    The overall normalization of the dark matter halo \|
+                    (kappa_s = (rho_0 * D_d * scale_radius)/lensing_critical_density)
+                scale_radius
+                    The cored NFW scale radius `theta_s`, as an angle on the sky in arcseconds.
+                core_radius
+                    The cored NFW core radius `theta_c`, as an angle on the sky in arcseconds.
+                """
+
+        super().__init__(
+            centre=centre,
+            ell_comps=(0.0, 0.0))
+
+        self.kappa_s = kappa_s
+        self.scale_radius = scale_radius
+        self.core_radius = core_radius
+
+
+    @aa.grid_dec.to_vector_yx
+    @aa.grid_dec.transform
+    def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        """
+        Calculate the deflection angles on a grid of (y,x) arc-second coordinates.
+
+        The input grid of (y,x) coordinates are transformed to a coordinate system centred on the profile centre with
+        and rotated based on the position angle defined from its `ell_comps` (this is described fully below).
+
+        The numerical backend can be selected via the ``xp`` argument, allowing this
+        method to be used with both NumPy and JAX (e.g. inside ``jax.jit``-compiled
+        code). This is described fully later in this example.
+
+        Parameters
+        ----------
+        grid
+            The grid of (y,x) arc-second coordinates the deflection angles are computed on.
+        xp
+            The numerical backend to use, either `numpy` or `jax.numpy`.
+        """
+        theta = self.radial_grid_from(grid=grid, xp=xp, **kwargs).array
+        theta = xp.maximum(theta, 1e-8)
+
+        factor = (
+            4.0
+            * self.kappa_s
+            * self.scale_radius**2
+        )
+
+        deflection_r = (
+            factor
+            * (self.F_func(theta, self.scale_radius, xp=xp) - self.F_func(theta, self.core_radius, xp=xp)
+               - (self.scale_radius - self.core_radius) * self.dev_F_func(theta, self.scale_radius, xp=xp)
+            )
+            / (theta * (self.scale_radius - self.core_radius)**2)
+        )
+
+
+        return self._cartesian_grid_via_radial_from(
+            grid=grid,
+            radius=deflection_r,
+            xp=xp,
+            **kwargs,
+        )
+
+    def F_func(self, theta, radius, xp=np):
+
+        F = theta * 0.0
+
+        # theta < radius
+        mask1 = (theta > 0) & (theta < radius)
+
+        # theta > radius
+        mask2 = theta > radius
+
+        F = xp.where(
+            mask1,
+            (
+                    radius / 2 * xp.log(2 * radius / theta)
+                    - xp.sqrt(radius ** 2 - theta ** 2)
+                    * xp.arctanh(xp.sqrt((radius - theta) / (radius + theta)))
+            ),
+            F,
+        )
+
+        F = xp.where(
+            mask2,
+            (
+                    radius / 2 * xp.log(2 * radius / theta)
+                    + xp.sqrt(theta ** 2 - radius ** 2)
+                    * xp.arctan(xp.sqrt((theta - radius) / (theta + radius)))
+            ),
+            F,
+        )
+
+        return 2 * radius * F
+
+    def dev_F_func(self, theta, radius, xp=np):
+
+        dev_F = theta * 0.0
+
+        mask1 = (theta > 0) & (theta < radius)
+        mask2 = theta == radius
+        mask3 = theta > radius
+
+        dev_F = xp.where(
+            mask1,
+            (
+                    radius * xp.log(2 * radius / theta)
+                    - (2 * radius ** 2 - theta ** 2) / xp.sqrt(radius ** 2 - theta ** 2)
+                    * xp.arctanh(xp.sqrt((radius - theta) / (radius + theta)))
+            ),
+            dev_F,
+        )
+
+        dev_F = xp.where(
+            mask2,
+            radius * (xp.log(2) - 1 / 2),
+            dev_F,
+        )
+
+        dev_F = xp.where(
+            mask3,
+            (
+                    radius * xp.log(2 * radius / theta)
+                    + (theta ** 2 - 2 * radius ** 2) / xp.sqrt(theta ** 2 - radius ** 2)
+                    * xp.arctan(xp.sqrt((theta - radius) / (theta + radius)))
+            ),
+            dev_F,
+        )
+
+        return 2 * dev_F
+
+    @aa.grid_dec.to_array
+    def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        """
+        Convergence (dimensionless surface mass density) for the cored NFW profile.
+        This is not yet implemented for `cNFWSph`.
+        """
+        raise NotImplementedError(
+            "convergence_2d_from is not implemented for cNFWSph; a physical cNFW "
+            "convergence expression must be added before use."
+        )
+
+    @aa.grid_dec.to_array
+    def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        """
+        Lensing potential for the cored NFW profile.
+        This is not yet implemented for `cNFWSph`.
+        """
+        raise NotImplementedError(
+            "potential_2d_from is not implemented for cNFWSph; a physical cNFW "
+            "potential expression must be added before use."
+        )

autogalaxy/profiles/mass/dark/cnfw_mcr.py

@@ -0,0 +1,34 @@
+from typing import Tuple
+
+from autogalaxy.profiles.mass.dark.cnfw import cNFWSph
+
+from autogalaxy.profiles.mass.dark import mcr_util
+
+class cNFWMCRLudlowSph(cNFWSph):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        mass_at_200: float = 1e9,
+        f_c = 0.01,
+        redshift_object: float = 0.5,
+        redshift_source: float = 1.0,
+    ):
+        self.mass_at_200 = mass_at_200
+        self.f_c = f_c
+        self.redshift_object = redshift_object
+        self.redshift_source = redshift_source
+
+        (
+            kappa_s,
+            scale_radius,
+            core_radius,
+            radius_at_200,
+        ) = mcr_util.kappa_s_scale_radius_and_core_radius_for_ludlow(
+            mass_at_200=mass_at_200,
+            scatter_sigma=0.0,
+            f_c=f_c,
+            redshift_object=redshift_object,
+            redshift_source=redshift_source,
+        )
+
+        super().__init__(centre=centre, kappa_s=kappa_s, scale_radius=scale_radius, core_radius=core_radius)