feat(mass): add ExternalPotential mass profile (Powell 2022 Eq 4)

New `ag.mp.ExternalPotential` line-of-sight profile that generalises
`ExternalShear` by adding the τ (linear surface-mass-density gradient,
spin-1) and δ (spin-3 generalised-shear) terms from Powell et al. 2022
Eq 4. Six free parameters plus a free centre (unlike ExternalShear's
fixed (0,0), since τ and δ have radial deflections). γ-only reduces
exactly to ExternalShear.

Includes a `from_magnitudes_and_angles` classmethod and per-term
magnitude/angle accessors. Convergence κ = τ₁·x + τ₂·y is non-zero for
the τ term (the prototype's κ=0 placeholder is wrong — derived from
Laplacian of ψ). 14 new unit tests; full mass-profiles suite (405
tests) still green.

Closes #419

Author: Jammy2211

autogalaxy/profiles/mass/__init__.py

@@ -60,4 +60,4 @@
     Chameleon,
     ChameleonSph,
 )
-from .sheets import ExternalShear, MassSheet
+from .sheets import ExternalPotential, ExternalShear, MassSheet

autogalaxy/profiles/mass/sheets/__init__.py

@@ -1,2 +1,3 @@
 from .mass_sheet import MassSheet
 from .external_shear import ExternalShear
+from .external_potential import ExternalPotential

autogalaxy/profiles/mass/sheets/external_potential.py

@@ -0,0 +1,224 @@
+from typing import Tuple
+
+import numpy as np
+
+import autoarray as aa
+
+from autogalaxy.profiles.mass.abstract.abstract import MassProfile
+
+
+class ExternalPotential(MassProfile):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        gamma_1: float = 0.0,
+        gamma_2: float = 0.0,
+        tau_1: float = 0.0,
+        tau_2: float = 0.0,
+        delta_1: float = 0.0,
+        delta_2: float = 0.0,
+    ):
+        r"""
+        A line-of-sight external potential that generalises the constant external shear used in
+        strong-lens models by adding the two next-order terms from Powell et al. 2022 (Eq. 4):
+
+        .. math::
+
+            \psi(\mathbf{r}) =
+                \tfrac{1}{2}\, r^2 \big(\gamma_1 \cos 2\theta + \gamma_2 \sin 2\theta\big)
+                + \tfrac{1}{4}\, r^3 \big(\tau_1 \cos\theta + \tau_2 \sin\theta\big)
+                + \tfrac{1}{6}\, r^3 \big(\delta_1 \cos 3\theta + \delta_2 \sin 3\theta\big)
+
+        where :math:`(r, \theta)` are polar coordinates centred on the profile's ``centre``.
+
+        Term-by-term:
+
+        - :math:`\gamma_1, \gamma_2` — the constant external shear contribution. With
+          :math:`\tau_i = \delta_i = 0` and ``centre = (0, 0)`` this reduces exactly to
+          :class:`ExternalShear`.
+        - :math:`\tau_1, \tau_2` — a linear gradient in the surface mass density (spin-1), giving a
+          non-zero convergence :math:`\kappa(x, y) = \tau_1 x + \tau_2 y`.
+        - :math:`\delta_1, \delta_2` — a higher-order spin-3 generalised-shear term.
+
+        Unlike :class:`ExternalShear`, where the deflection field is a constant in the lens plane
+        and the source position is degenerate with ``centre``, the :math:`\tau` and :math:`\delta`
+        deflections have explicit radial dependence — so ``centre`` is a free parameter (typically
+        tied to the primary lens centre when modelling).
+
+        Parameters
+        ----------
+        centre
+            The (y, x) arc-second coordinates of the profile centre.
+        gamma_1, gamma_2
+            The two components of the constant external shear (spin-2).
+        tau_1, tau_2
+            The two components of the linear surface-mass-density gradient (spin-1).
+        delta_1, delta_2
+            The two components of the higher-order spin-3 generalised shear.
+        """
+
+        super().__init__(centre=centre, ell_comps=(0.0, 0.0))
+        self.gamma_1 = gamma_1
+        self.gamma_2 = gamma_2
+        self.tau_1 = tau_1
+        self.tau_2 = tau_2
+        self.delta_1 = delta_1
+        self.delta_2 = delta_2
+
+    @staticmethod
+    def _magnitude_from(c1, c2, xp=np):
+        return xp.sqrt(c1 * c1 + c2 * c2)
+
+    @staticmethod
+    def _angle_from(c1, c2, harmonic: int, xp=np):
+        r"""
+        Return the principal angle in degrees for a harmonic of order ``harmonic``.
+
+        - ``harmonic = 1`` -> angle in [0, 360)
+        - ``harmonic = 2`` -> angle in [0, 180) (shear convention)
+        - ``harmonic = 3`` -> angle in [0, 120)
+        """
+        angle = xp.rad2deg(xp.arctan2(c2, c1)) / harmonic
+        period = 360.0 / harmonic
+        return angle % period
+
+    @classmethod
+    def from_magnitudes_and_angles(
+        cls,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        gamma: float = 0.0,
+        theta_gamma: float = 0.0,
+        tau: float = 0.0,
+        theta_tau: float = 0.0,
+        delta: float = 0.0,
+        theta_delta: float = 0.0,
+    ):
+        r"""
+        Build an :class:`ExternalPotential` from per-term magnitudes and position angles, matching
+        the paper-style parameterisation. Angles are in degrees, anticlockwise from the +x axis.
+        """
+        tg = np.deg2rad(theta_gamma)
+        tt = np.deg2rad(theta_tau)
+        td = np.deg2rad(theta_delta)
+
+        gamma_1 = gamma * np.cos(2.0 * tg)
+        gamma_2 = gamma * np.sin(2.0 * tg)
+        tau_1 = tau * np.cos(tt)
+        tau_2 = tau * np.sin(tt)
+        delta_1 = delta * np.cos(3.0 * td)
+        delta_2 = delta * np.sin(3.0 * td)
+
+        return cls(
+            centre=centre,
+            gamma_1=gamma_1,
+            gamma_2=gamma_2,
+            tau_1=tau_1,
+            tau_2=tau_2,
+            delta_1=delta_1,
+            delta_2=delta_2,
+        )
+
+    def gamma_magnitude(self, xp=np):
+        return self._magnitude_from(self.gamma_1, self.gamma_2, xp=xp)
+
+    def gamma_angle(self, xp=np):
+        return self._angle_from(self.gamma_1, self.gamma_2, harmonic=2, xp=xp)
+
+    def tau_magnitude(self, xp=np):
+        return self._magnitude_from(self.tau_1, self.tau_2, xp=xp)
+
+    def tau_angle(self, xp=np):
+        return self._angle_from(self.tau_1, self.tau_2, harmonic=1, xp=xp)
+
+    def delta_magnitude(self, xp=np):
+        return self._magnitude_from(self.delta_1, self.delta_2, xp=xp)
+
+    def delta_angle(self, xp=np):
+        return self._angle_from(self.delta_1, self.delta_2, harmonic=3, xp=xp)
+
+    def convergence_func(self, grid_radius: float) -> float:
+        return 0.0
+
+    def average_convergence_of_1_radius(self):
+        return 0.0
+
+    @aa.decorators.to_array
+    @aa.decorators.transform
+    def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        r"""
+        Returns the convergence :math:`\kappa = \tfrac{1}{2}\nabla^2 \psi` at each grid point.
+
+        Only the :math:`\tau` term contributes; the :math:`\gamma` (spin-2 shear) and
+        :math:`\delta` (spin-3) terms are harmonic and yield zero convergence:
+
+        .. math::
+
+            \kappa(x, y) = \tau_1 \, x + \tau_2 \, y
+
+        where :math:`(x, y)` are coordinates relative to ``centre``.
+        """
+        y = grid.array[:, 0]
+        x = grid.array[:, 1]
+        return self.tau_1 * x + self.tau_2 * y
+
+    @aa.decorators.to_array
+    @aa.decorators.transform
+    def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        r"""
+        Returns the lensing potential of the external potential at each grid point, following
+        Powell et al. 2022 Eq. 4.
+        """
+        y = grid.array[:, 0]
+        x = grid.array[:, 1]
+        r2 = x * x + y * y
+        r = xp.sqrt(r2)
+        theta = xp.arctan2(y, x)
+
+        gamma_term = 0.5 * r2 * (
+            self.gamma_1 * xp.cos(2.0 * theta) + self.gamma_2 * xp.sin(2.0 * theta)
+        )
+        tau_term = 0.25 * r2 * r * (
+            self.tau_1 * xp.cos(theta) + self.tau_2 * xp.sin(theta)
+        )
+        delta_term = (1.0 / 6.0) * r2 * r * (
+            self.delta_1 * xp.cos(3.0 * theta) + self.delta_2 * xp.sin(3.0 * theta)
+        )
+
+        return gamma_term + tau_term + delta_term
+
+    @aa.decorators.to_vector_yx
+    @aa.decorators.transform
+    def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        r"""
+        Returns the deflection vector :math:`\boldsymbol{\alpha} = \nabla \psi` at each grid point.
+
+        Computed in polar form (`alpha_r = d\psi/dr`, `alpha_theta = (1/r)\,d\psi/d\theta`) and
+        then projected back to ``(y, x)`` Cartesian components. ``centre`` enters through the
+        ``@transform`` decorator, which shifts the grid before the function body runs.
+        """
+        y = grid.array[:, 0]
+        x = grid.array[:, 1]
+        r = xp.sqrt(x * x + y * y)
+        theta = xp.arctan2(y, x)
+        cos_t = xp.cos(theta)
+        sin_t = xp.sin(theta)
+        cos_2t = xp.cos(2.0 * theta)
+        sin_2t = xp.sin(2.0 * theta)
+        cos_3t = xp.cos(3.0 * theta)
+        sin_3t = xp.sin(3.0 * theta)
+
+        alpha_r = (
+            r * (self.gamma_1 * cos_2t + self.gamma_2 * sin_2t)
+            + 0.75 * r * r * (self.tau_1 * cos_t + self.tau_2 * sin_t)
+            + 0.5 * r * r * (self.delta_1 * cos_3t + self.delta_2 * sin_3t)
+        )
+        alpha_theta = (
+            r * (-self.gamma_1 * sin_2t + self.gamma_2 * cos_2t)
+            + 0.25 * r * r * (-self.tau_1 * sin_t + self.tau_2 * cos_t)
+            + 0.5 * r * r * (-self.delta_1 * sin_3t + self.delta_2 * cos_3t)
+        )
+
+        alpha_y = sin_t * alpha_r + cos_t * alpha_theta
+        alpha_x = cos_t * alpha_r - sin_t * alpha_theta
+
+        return xp.vstack((alpha_y, alpha_x)).T

test_autogalaxy/profiles/mass/sheets/test_external_potential.py

@@ -0,0 +1,126 @@
+import autogalaxy as ag
+import numpy as np
+import pytest
+
+
+grid_diag = ag.Grid2DIrregular([[0.1, 0.1]])
+grid_unit_x = ag.Grid2DIrregular([[0.0, 1.0]])
+grid_multi = ag.Grid2DIrregular([[0.1, 0.1], [0.2, 0.2], [0.3, 0.3]])
+
+
+def test__convergence_2d_from__gamma_only_is_zero():
+    mp = ag.mp.ExternalPotential(gamma_1=0.1, gamma_2=-0.05)
+    convergence = mp.convergence_2d_from(grid=grid_multi)
+    assert convergence == pytest.approx(np.zeros(3), 1.0e-8)
+
+
+def test__convergence_2d_from__delta_only_is_zero():
+    mp = ag.mp.ExternalPotential(delta_1=0.07, delta_2=0.04)
+    convergence = mp.convergence_2d_from(grid=grid_multi)
+    assert convergence == pytest.approx(np.zeros(3), 1.0e-8)
+
+
+def test__convergence_2d_from__tau_only_is_linear_in_xy():
+    mp = ag.mp.ExternalPotential(tau_1=0.1, tau_2=0.2)
+    grid = ag.Grid2DIrregular([[0.0, 1.0], [0.5, 1.0], [1.0, 0.0]])
+    convergence = mp.convergence_2d_from(grid=grid)
+    # kappa(x, y) = tau_1 * x + tau_2 * y
+    expected = np.array([0.1 * 1.0 + 0.2 * 0.0, 0.1 * 1.0 + 0.2 * 0.5, 0.1 * 0.0 + 0.2 * 1.0])
+    assert convergence == pytest.approx(expected, 1.0e-6)
+
+
+def test__convergence_2d_from__tau_with_nonzero_centre_shifts_origin():
+    mp = ag.mp.ExternalPotential(centre=(1.0, 1.0), tau_1=0.1, tau_2=0.0)
+    grid = ag.Grid2DIrregular([[1.0, 1.0], [1.0, 2.0]])
+    convergence = mp.convergence_2d_from(grid=grid)
+    assert convergence == pytest.approx(np.array([0.0, 0.1]), 1.0e-6)
+
+
+def test__potential_2d_from__gamma_only_matches_external_shear():
+    shear = ag.mp.ExternalShear(gamma_1=0.1, gamma_2=-0.05)
+    pot = ag.mp.ExternalPotential(gamma_1=0.1, gamma_2=-0.05)
+    expected = shear.potential_2d_from(grid=grid_multi)
+    actual = pot.potential_2d_from(grid=grid_multi)
+    assert actual == pytest.approx(np.asarray(expected), 1.0e-6)
+
+
+def test__potential_2d_from__tau_only_unit_x_axis():
+    mp = ag.mp.ExternalPotential(tau_1=0.05, tau_2=0.0)
+    # at (y=0, x=1): r=1, theta=0 -> psi = 0.25 * 1^3 * (0.05 * 1 + 0) = 0.0125
+    potential = mp.potential_2d_from(grid=grid_unit_x)
+    assert potential == pytest.approx(np.array([0.0125]), 1.0e-6)
+
+
+def test__potential_2d_from__delta_only_unit_x_axis():
+    mp = ag.mp.ExternalPotential(delta_1=0.1, delta_2=0.0)
+    # at (y=0, x=1): r=1, theta=0 -> psi = (1/6) * 1^3 * (0.1 * 1 + 0) = 0.1/6
+    potential = mp.potential_2d_from(grid=grid_unit_x)
+    assert potential == pytest.approx(np.array([0.1 / 6.0]), 1.0e-6)
+
+
+def test__deflections_yx_2d_from__gamma_only_matches_external_shear():
+    shear = ag.mp.ExternalShear(gamma_1=-0.17320, gamma_2=0.1)
+    pot = ag.mp.ExternalPotential(gamma_1=-0.17320, gamma_2=0.1)
+    grid = ag.Grid2DIrregular([[0.1625, 0.1625]])
+    expected = np.asarray(shear.deflections_yx_2d_from(grid=grid))
+    actual = np.asarray(pot.deflections_yx_2d_from(grid=grid))
+    assert actual == pytest.approx(expected, 1.0e-5)
+
+
+def test__deflections_yx_2d_from__tau_only_radial_unit_x_axis():
+    mp = ag.mp.ExternalPotential(tau_1=0.05, tau_2=0.0)
+    # at (y=0, x=1): r=1, theta=0
+    #   alpha_r = 0.75 * 1 * 0.05 = 0.0375 ; alpha_theta = 0
+    #   alpha_y = 0, alpha_x = 0.0375
+    deflections = mp.deflections_yx_2d_from(grid=grid_unit_x)
+    assert deflections[0, 0] == pytest.approx(0.0, 1.0e-6)
+    assert deflections[0, 1] == pytest.approx(0.0375, 1.0e-6)
+
+
+def test__deflections_yx_2d_from__delta_only_unit_x_axis():
+    mp = ag.mp.ExternalPotential(delta_1=0.1, delta_2=0.0)
+    # at (y=0, x=1): alpha_r = 0.5 * 1 * 0.1 = 0.05 ; alpha_theta = 0
+    deflections = mp.deflections_yx_2d_from(grid=grid_unit_x)
+    assert deflections[0, 0] == pytest.approx(0.0, 1.0e-6)
+    assert deflections[0, 1] == pytest.approx(0.05, 1.0e-6)
+
+
+def test__deflections_yx_2d_from__nonzero_centre_shifts_origin():
+    mp = ag.mp.ExternalPotential(centre=(1.0, 1.0), gamma_1=0.1, gamma_2=0.0)
+    # at (y=1, x=2) post-shift becomes (y=0, x=1): r=1, theta=0
+    #   alpha_r = 1 * (0.1 * cos(0) + 0) = 0.1 ; alpha_theta = 0
+    #   alpha_y = 0, alpha_x = 0.1
+    grid = ag.Grid2DIrregular([[1.0, 2.0]])
+    deflections = mp.deflections_yx_2d_from(grid=grid)
+    assert deflections[0, 0] == pytest.approx(0.0, 1.0e-6)
+    assert deflections[0, 1] == pytest.approx(0.1, 1.0e-6)
+
+
+def test__from_magnitudes_and_angles__roundtrip_gamma():
+    gamma = 0.1
+    theta_gamma = 30.0  # degrees
+    mp = ag.mp.ExternalPotential.from_magnitudes_and_angles(
+        gamma=gamma, theta_gamma=theta_gamma
+    )
+    assert mp.gamma_magnitude() == pytest.approx(gamma, 1.0e-6)
+    assert mp.gamma_angle() == pytest.approx(theta_gamma, 1.0e-6)
+
+
+def test__from_magnitudes_and_angles__roundtrip_tau():
+    tau = 0.05
+    theta_tau = 200.0  # degrees, spin-1 so [0, 360)
+    mp = ag.mp.ExternalPotential.from_magnitudes_and_angles(
+        tau=tau, theta_tau=theta_tau
+    )
+    assert mp.tau_magnitude() == pytest.approx(tau, 1.0e-6)
+    assert mp.tau_angle() == pytest.approx(theta_tau, 1.0e-6)
+
+
+def test__from_magnitudes_and_angles__roundtrip_delta():
+    delta = 0.02
+    theta_delta = 40.0  # degrees, spin-3 so [0, 120)
+    mp = ag.mp.ExternalPotential.from_magnitudes_and_angles(
+        delta=delta, theta_delta=theta_delta
+    )
+    assert mp.delta_magnitude() == pytest.approx(delta, 1.0e-6)
+    assert mp.delta_angle() == pytest.approx(theta_delta, 1.0e-6)