Merge branch 'main' of github.com:Jammy2211/PyAutoGalaxy

Jammy2211 · Jammy2211 · commit 3ca52aef31fe · 2026-02-13T12:58:35.000Z
diff --git a/autogalaxy/config/priors/mass/dark/cnfw_mcr_scatter.yaml b/autogalaxy/config/priors/mass/dark/cnfw_mcr_scatter.yaml
@@ -0,0 +1,71 @@
+cNFWMCRScatterLudlowSph:
+  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
+  scatter:
+    type: Gaussian
+    mean: 0.0
+    sigma: 3.0
+    width_modifier:
+      type: Absolute
+      value: 1.0
+    limits:
+      lower: -inf
+      upper: inf
diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py
@@ -38,6 +38,7 @@
     NFWVirialMassConcSph,
     cNFWSph,
     cNFWMCRLudlowSph,
+    cNFWMCRScatterLudlowSph,
 )
 from .stellar import (
     Gaussian,
diff --git a/autogalaxy/profiles/mass/dark/__init__.py b/autogalaxy/profiles/mass/dark/__init__.py
@@ -11,3 +11,4 @@
 from .nfw_virial_mass_conc import NFWVirialMassConcSph
 from .cnfw import cNFWSph
 from .cnfw_mcr import cNFWMCRLudlowSph
+from .cnfw_mcr_scatter import cNFWMCRScatterLudlowSph
diff --git a/autogalaxy/profiles/mass/dark/cnfw.py b/autogalaxy/profiles/mass/dark/cnfw.py
@@ -156,18 +156,12 @@ 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."
-        )
+        return xp.zeros(shape=grid.shape[0])
 
     @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."
-        )
+        return xp.zeros(shape=grid.shape[0])
diff --git a/autogalaxy/profiles/mass/dark/cnfw_mcr_scatter.py b/autogalaxy/profiles/mass/dark/cnfw_mcr_scatter.py
@@ -0,0 +1,42 @@
+from typing import Tuple
+
+from autogalaxy.profiles.mass.dark.cnfw import cNFWSph
+
+from autogalaxy.profiles.mass.dark import mcr_util
+
+
+class cNFWMCRScatterLudlowSph(cNFWSph):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        mass_at_200: float = 1e9,
+        scatter_sigma: float = 0.0,
+        f_c=0.01,
+        redshift_object: float = 0.5,
+        redshift_source: float = 1.0,
+    ):
+        self.mass_at_200 = mass_at_200
+        self.scatter_sigma = scatter_sigma
+        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=scatter_sigma,
+            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,
+        )
diff --git a/autogalaxy/profiles/mass/dark/nfw.py b/autogalaxy/profiles/mass/dark/nfw.py
@@ -42,8 +42,63 @@ def __init__(
         )
         super(MassProfileCSE, self).__init__()
 
-    def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike):
-        return self.deflections_2d_via_cse_from(grid=grid)
+    def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        return self.deflections_2d_via_analytic_from(grid=grid, xp=xp, **kwargs)
+
+    @aa.grid_dec.to_vector_yx
+    @aa.grid_dec.transform
+    def deflections_2d_via_analytic_from(
+            self, grid: aa.type.Grid2DLike, xp=np, **kwargs
+    ):
+        """
+        Analytic calculation deflection angles from Heyrovský & Karamazov 2024 via Eq. 30 & 31
+
+        Parameters
+        ----------
+        grid
+            The grid of (y,x) arc-second coordinates the deflection angles are computed on.
+        """
+
+        # Convert e definitions:
+        # from q = (1-e)/(1+e) to q = sqrt(1-e**2)
+
+        e_autolens = xp.sqrt(self.ell_comps[1] ** 2 + self.ell_comps[0] ** 2)
+        e_hk24 = 2 * xp.sqrt(e_autolens) / (1 + e_autolens)
+
+        # Define dimensionless length coords
+
+        x1 = grid.array[:, 1] / self.scale_radius
+        x2 = grid.array[:, 0] / self.scale_radius
+
+        r2 = x1 ** 2 + x2 ** 2
+
+        # Avoid nans
+
+        mask = r2 > 1e-24
+
+        prefactor = xp.where(mask, 4 * self.kappa_s * xp.sqrt(1 - e_hk24 ** 2) / (
+                ((x1 - e_hk24) ** 2 + x2 ** 2) * ((x1 + e_hk24) ** 2 + x2 ** 2)
+        ), 0.0)
+
+        f1 = xp.where(mask, nfw_hk24_util.small_f_1(x1, x2, e_hk24, xp=xp), 0.0)
+        f2 = xp.where(mask, nfw_hk24_util.small_f_2(x1, x2, e_hk24, xp=xp), 0.0)
+        f3 = xp.where(mask, nfw_hk24_util.small_f_3(x1, x2, e_hk24, xp=xp), 0.0)
+
+        deflection_x = (x1 * ((x1**2 - e_hk24**2) * (1 - e_hk24**2) + x2**2 * (1 + e_hk24**2)) * f1
+                        + x1 * (x1**2 + x2**2 - e_hk24**2) * f2
+                        - x2 * (x1**2 + x2**2 + e_hk24**2) * f3)
+        deflection_x *= prefactor
+
+        deflection_y = (x2 * (x1**2 * (1 - 2 * e_hk24**2) + x2**2 + e_hk24**2) * f1
+                        + x2 * (x1**2 + x2**2 + e_hk24**2) * f2
+                        + x1 * (x1**2 + x2**2 - e_hk24**2) * f3)
+        deflection_y *= prefactor
+
+        return self.rotated_grid_from_reference_frame_from(
+            xp.multiply(self.scale_radius, xp.vstack((deflection_y, deflection_x)).T),
+            xp=xp,
+            **kwargs,
+        )
 
     @aa.grid_dec.to_vector_yx
     @aa.grid_dec.transform
@@ -146,7 +201,7 @@ def convergence_func(self, grid_radius: float) -> float:
         return np.real(
             2.0
             * self.kappa_s
-            * np.array(self.coord_func_g(grid_radius=grid_radius, xp=xp))
+            * np.array(self.coord_func_g(grid_radius=grid_radius))
         )
 
     @aa.over_sample
@@ -271,26 +326,26 @@ def shear_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
         # Convert e definitions:
         # from q = (1-e)/(1+e) to q = sqrt(1-e**2)
 
-        e_autolens = np.sqrt(self.ell_comps[1] ** 2 + self.ell_comps[0] ** 2)
-        e_hk24 = 2 * np.sqrt(e_autolens) / np.sqrt(1 + 2 * e_autolens + e_autolens**2)
+        e_autolens = xp.sqrt(self.ell_comps[1] ** 2 + self.ell_comps[0] ** 2)
+        e_hk24 = 2 * xp.sqrt(e_autolens) / (1 + e_autolens)
 
         # Define dimensionless length coords
 
         x1 = grid.array[:, 1] / self.scale_radius
         x2 = grid.array[:, 0] / self.scale_radius
 
         # Avoid nans due to x=0
-        x1 = np.where(np.abs(x1) < 1e-6, 1e-6, x1)
-        x2 = np.where(np.abs(x2) < 1e-6, 1e-6, x2)
+        x1 = xp.where(xp.abs(x1) < 1e-6, 1e-6, x1)
+        x2 = xp.where(xp.abs(x2) < 1e-6, 1e-6, x2)
 
         # Calculate shear from nfw_HK24.py
 
-        g1, g2 = nfw_hk24_util.g1_g2_from(x1=x1, x2=x2, e=e_hk24, k_s=self.kappa_s)
+        g1, g2 = nfw_hk24_util.g1_g2_from(x1=x1, x2=x2, e=e_hk24, k_s=self.kappa_s, xp=xp)
 
         # Rotation for shear
 
         shear_field = self.rotated_grid_from_reference_frame_from(
-            grid=np.vstack((g2, g1)).T, angle=self.angle(xp) * 2
+            grid=xp.vstack((g2, g1)).T, angle=self.angle(xp) * 2, xp=xp
         )
 
         return aa.VectorYX2DIrregular(values=shear_field, grid=grid)
@@ -315,8 +370,8 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
         # Convert e definitions:
         # from q = (1-e)/(1+e) to q = sqrt(1-e**2)
 
-        e_autolens = np.sqrt(self.ell_comps[1] ** 2 + self.ell_comps[0] ** 2)
-        e_hk24 = 2 * np.sqrt(e_autolens) / np.sqrt(1 + 2 * e_autolens + e_autolens**2)
+        e_autolens = xp.sqrt(self.ell_comps[1] ** 2 + self.ell_comps[0] ** 2)
+        e_hk24 = 2 * xp.sqrt(e_autolens) / (1 + e_autolens)
 
         # Define dimensionless length coords
 
@@ -325,13 +380,13 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
 
         # Avoid nans due to x=0
 
-        x1 = np.where(np.abs(x1) < 1e-6, 1e-6, x1)
-        x2 = np.where(np.abs(x2) < 1e-6, 1e-6, x2)
+        x1 = xp.where(xp.abs(x1) < 1e-6, 1e-6, x1)
+        x2 = xp.where(xp.abs(x2) < 1e-6, 1e-6, x2)
 
         # Calculate convergence from nfw_HK24.py
-        a = nfw_hk24_util.semi_major_axis_from(x1, x2, e_hk24)
+        a = nfw_hk24_util.semi_major_axis_from(x1, x2, e_hk24, xp=xp)
 
-        return nfw_hk24_util.kappa_from(k_s=self.kappa_s, a=a)
+        return nfw_hk24_util.kappa_from(k_s=self.kappa_s, a=a, xp=xp)
 
 
 class NFWSph(NFW):
diff --git a/autogalaxy/profiles/mass/dark/nfw_hk24_util.py b/autogalaxy/profiles/mass/dark/nfw_hk24_util.py
diff --git a/test_autogalaxy/profiles/mass/dark/test_nfw.py b/test_autogalaxy/profiles/mass/dark/test_nfw.py
diff --git a/test_autogalaxy/profiles/mass/dark/test_nfw_scatter.py b/test_autogalaxy/profiles/mass/dark/test_nfw_scatter.py

Original file line number	Diff line number	Diff line change
`@@ -38,6 +38,7 @@`
`38`	`38`	`NFWVirialMassConcSph,`
`39`	`39`	`cNFWSph,`
`40`	`40`	`cNFWMCRLudlowSph,`
	`41`	`+ cNFWMCRScatterLudlowSph,`
`41`	`42`	`)`
`42`	`43`	`from .stellar import (`
`43`	`44`	`Gaussian,`