External potential

[Sketos]

Overview

This profile is described in Powell 2022 and function form of the potential is given in Eq 4. The first term is equivalent to the external shear that we use in autolens. The other terms correspond to a gradient in the surface mass density and gradient of the external shear.

Plan

I think the best approach for the code provided is to commit the first term of the ExternalPotential as this is exactly the ExternalShear in autolens. The deflection angles for the two other terms have a radial dependance so the centre should also be a free parameter (which would often be set to external_potential.centre = mass.centre), unlike the ExternalShear where deflection angles are constant and therefore the position of the source is degenerate with them (centre can be fixed to 0,0).

Detailed implementation plan

Add a file external_potential.py

Affected Repositories

profiles/mass/sheets

Implementation Steps

Key Files

Example Code

class ExternalPotential(MassProfile):
    def __init__(
        self,
        centre=(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,
    ):
        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 radial_and_angle_grid_from(
        grid: aa.type.Grid2DLike, centre: Tuple[float, float] = (0.0, 0.0), xp=np
    ) -> Tuple[np.ndarray, np.ndarray]:
        y, x = grid.array.T

        x_shifted = xp.subtract(x, centre[1])
        y_shifted = xp.subtract(y, centre[0])

        radial_grid = xp.sqrt(x_shifted**2 + y_shifted**2)

        angle_grid = xp.arctan2(y_shifted, x_shifted)

        return radial_grid, angle_grid

    @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):
        """
        Return the principal angle in degrees for a harmonic of order `harmonic`.

        harmonic = 2 -> [0, 180)
        harmonic = 1 -> [0, 360)
        harmonic = 3 -> [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=(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,
    ):
        """
        Build the profile from paper-style magnitudes and angles.
        Angles are in degrees, anticlockwise from +x.
        """
        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)

    @aa.grid_dec.to_array
    def convergence_2d_from(self, grid, xp=np, **kwargs):
        return xp.zeros(shape=grid.shape[0])

    @aa.grid_dec.to_array
    def potential_2d_from(self, grid, xp=np, **kwargs):
        r, theta = self.radial_and_angle_grid_from(
            grid=grid, 
            centre=self.centre, 
            xp=xp
        )
        gamma_term = 0.5 * r**2 * (
            self.gamma_1 * xp.cos(2.0 * theta) + self.gamma_2 * xp.sin(2.0 * theta)
        )
        tau_term = 0.25 * r**3 * (
            self.tau_1 * xp.cos(theta) + self.tau_2 * xp.sin(theta)
        )
        delta_term = (1.0 / 6.0) * r**3 * (
            self.delta_1 * xp.cos(3.0 * theta) + self.delta_2 * xp.sin(3.0 * theta)
        )

        return gamma_term + tau_term + delta_term

    @aa.grid_dec.to_vector_yx
    @aa.grid_dec.transform
    def deflections_yx_2d_from(self, grid, xp=np, **kwargs):
        r, theta = self.radial_and_angle_grid_from(
            grid=grid, 
            centre=self.centre, 
            xp=xp
        )
        alpha_r = (
            r * (
                self.gamma_1 * xp.cos(2.0 * theta)
                + self.gamma_2 * xp.sin(2.0 * theta)
            )
            + 0.75 * r**2 * (
                self.tau_1 * xp.cos(theta) + self.tau_2 * xp.sin(theta)
            )
            + 0.5 * r**2 * (
                self.delta_1 * xp.cos(3.0 * theta)
                + self.delta_2 * xp.sin(3.0 * theta)
            )
        )
        alpha_theta = (
            r * (
                -self.gamma_1 * xp.sin(2.0 * theta)
                + self.gamma_2 * xp.cos(2.0 * theta)
            )
            + 0.25 * r**2 * (
                -self.tau_1 * xp.sin(theta) + self.tau_2 * xp.cos(theta)
            )
            + 0.5 * r**2 * (
                -self.delta_1 * xp.sin(3.0 * theta)
                + self.delta_2 * xp.cos(3.0 * theta)
            )
        )
        alpha_y = xp.sin(theta) * alpha_r + xp.cos(theta) * alpha_theta
        alpha_x = xp.cos(theta) * alpha_r - xp.sin(theta) * alpha_theta

        return xp.vstack((alpha_y, alpha_x)).T

Original Prompt

Click to expand starting prompt

[Jammy2211]

Hi @Sketos — thanks for the report. I'm Jammy's CLI assistant; I'll take a look and post back shortly with a plan, or follow-up questions if I need more detail.

[Jammy2211]

Here's what I'm planning to do — let me know if anything looks off.

The new ExternalPotential profile generalises ExternalShear by adding the two Powell 2022 (Eq 4) radial-gradient terms — τ (gradient of the surface mass density) and δ (gradient of the external shear). The γ-only subset matches the existing ExternalShear math exactly. I'll add it as a sibling class to ExternalShear in autogalaxy/profiles/mass/sheets/, with the same xp/decorator conventions and a unit-test file mirroring test_external_shear.py.

Plan

• Add ag.mp.ExternalPotential(centre, gamma_1, gamma_2, tau_1, tau_2, delta_1, delta_2) — six free parameters; centre is free (unlike ExternalShear's fixed (0, 0)) because the τ and δ deflections are radial.
• Implement convergence_2d_from, potential_2d_from, deflections_yx_2d_from following Eq 4 of Powell 2022, reusing the existing ExternalShear decorator/xp pattern.
• Add magnitude/angle convenience accessors (gamma_magnitude, tau_magnitude, delta_magnitude, plus their angle counterparts) and a from_magnitudes_and_angles classmethod matching your prototype.
• Wire exports through sheets/__init__.py → mass/__init__.py → ag.mp.*.
• Add tests in test_autogalaxy/profiles/mass/sheets/test_external_potential.py: parity with ExternalShear when τ=δ=0, separate τ-only and δ-only checks, and a non-zero-centre case.
• No workspace demo in this PR — happy to follow up with one in autogalaxy_workspace/ if you'd like.

Detailed implementation plan

Affected repositories

• PyAutoGalaxy (primary, library)

Work classification

Library work — purely additive (new class). No workspace impact expected.

Branch survey

Repository	Current Branch	Dirty?
PyAutoGalaxy	main	clean

Suggested branch: feature/external-potential
Worktree root: ~/Code/PyAutoLabs-wt/external-potential/ (created later by /start_library)

Implementation steps

1. Create autogalaxy/profiles/mass/sheets/external_potential.py — six-parameter sibling of ExternalShear, polar-form implementations of potential and deflections per Powell 2022 Eq 4, plus from_magnitudes_and_angles classmethod and per-term magnitude/angle accessors.
2. Match the actual codebase decorator imports (@aa.decorators.to_array, @aa.decorators.to_vector_yx, @aa.decorators.transform) rather than the @aa.grid_dec.* form in the prototype — the latter is in the CLAUDE.md style guide but the working ExternalShear uses aa.decorators.*.
3. Fix a subtle double-shift in the prototype: @aa.decorators.transform already shifts the grid by centre, so the body should compute r, θ from grid.array[:, 0/1] directly without re-subtracting self.centre.
4. Thread xp=np consistently across all methods and helpers.
5. Wire exports through sheets/__init__.py, mass/__init__.py, and confirm ag.mp.ExternalPotential resolves.
6. Tests: γ-only parity vs ExternalShear (potential + deflections), τ-only with non-zero centre, δ-only, zero-convergence-everywhere, from_magnitudes_and_angles round-trip.
7. Run pytest test_autogalaxy/profiles/mass/sheets/ and the full suite before opening the PR.

Key files

• autogalaxy/profiles/mass/sheets/external_potential.py — new
• autogalaxy/profiles/mass/sheets/__init__.py — add export
• autogalaxy/profiles/mass/__init__.py — add export
• test_autogalaxy/profiles/mass/sheets/test_external_potential.py — new

Decisions and trade-offs

• Sibling, not subclass of ExternalShear. Your prototype is a sibling and that reads cleaner: ExternalShear hardcodes centre=(0, 0) and uses a magnitude/angle parameterisation, while ExternalPotential needs free centre and (γ₁, γ₂)-direct inputs. We doc the γ-only equivalence in the docstring.
• Three terms, one PR. You suggested phasing γ first then τ/δ. The implementation is small enough that one PR keeps the math reviewable as a unit; happy to split into commits inside that PR if it helps review.
• Workspace demo deferred. Class is purely additive. Open a follow-up issue if you'd like a worked example.

I'll post back when the work is in progress and again when there's a PR. If anything in the plan looks off, just reply and I'll adjust.

[Jammy2211]

Smoke tests are all green across every workspace — 42 scripts ran without a single regression. We're good to merge.

Workspace	Passed	Skipped	Failed	Duration
autofit_workspace	6	1	0	4.9s
autogalaxy_workspace	6	0	0	25.4s
autolens_workspace	7	0	0	28.9s
autolens_workspace_test	11	1	0	167.0s
euclid_strong_lens_modeling_pipeline	6	0	0	151.5s
HowToLens	6	0	0	7.9s
Total	42	2	0	—

The two skips are both pre-existing entries in no_run.yaml unrelated to this PR (autofit_workspace/searches/mcmc.py — zeus section breaks test mode; autolens_workspace_test/database/scrape/general.py — NEEDS_FIX since 2026-04-27).

The new ag.mp.ExternalPotential class was validated against the worktree PyAutoGalaxy throughout all runs. No shared code paths were disturbed. Proceeding to merge.

[Jammy2211]

Shipped — 2026-05-18

PR

• feat(mass): add ExternalPotential mass profile (Powell 2022 Eq 4) #422 (merged)

Summary

ag.mp.ExternalPotential is live on main. The class implements all three Powell 2022 (Eq 4) terms in one go — γ (matches ExternalShear exactly when τ=δ=0), τ (linear surface-mass-density gradient), and δ (spin-3 generalised shear) — with a free centre and convenience accessors. Pip-installing the next PyAutoGalaxy release will give you the class; for now you can use it by running off main.

Session Notes

• One math correction vs the prototype: convergence has a non-zero contribution from τ (κ = τ₁·x + τ₂·y) — the prototype's convergence_2d_from = zeros would have been wrong for users who plot κ or compute mass within an aperture. γ and δ remain κ=0 (harmonic). If Powell 2022 uses a normalisation where the τ-convergence vanishes, please reply on this issue and we'll adjust.
• All 42 smoke-test scripts across the six PyAutoGalaxy/PyAutoLens workspaces still pass with the worktree active — no downstream regressions.
• No workspace demo in this PR — happy to open a follow-up if you'd like a worked example (autogalaxy_workspace/scripts/imaging/features/external_potential.py).

Thanks again for the prototype — it was tight enough that the implementation could land in a single PR. Let us know if you spot anything in the merged code.

[Jammy2211]

Small follow-up landed: PyAutoGalaxy#423 added the library-default prior YAML for ExternalPotential, so you can now build a model directly without a config error:

import autofit as af
import autogalaxy as ag

shear_potential = af.Model(ag.mp.ExternalPotential)
# 8 free params: centre + (gamma_1, gamma_2, tau_1, tau_2, delta_1, delta_2)

Defaults: centre Gaussian(0, 0.1); the six γ/τ/δ components Uniform(-0.3, 0.3) with Absolute width modifier 0.05. Tighter ranges are easy to override in your workspace config/priors/mass/sheets/external_potential.yaml. Also added a JAX JIT parity test in autolens_workspace_test#101 so the class is now covered by our integration suite.

Both shipped; nothing else needed from you. Reopen this issue if anything's off.