diff --git a/autogalaxy/config/notation.yaml b/autogalaxy/config/notation.yaml
index e4481b45..50d6bdde 100644
--- a/autogalaxy/config/notation.yaml
+++ b/autogalaxy/config/notation.yaml
@@ -69,6 +69,9 @@ label:
     zeroth_signal_scale: V
   superscript:
     ExternalShear: ext
+    GaussianRandomField: grf
+    InputDeflections: input
+    InputPotential: input
     Mesh: mesh
     Point: point
     SMBH: smbh
diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py
index 4468fb5a..765dfe2b 100644
--- a/autogalaxy/profiles/mass/__init__.py
+++ b/autogalaxy/profiles/mass/__init__.py
@@ -66,3 +66,9 @@
     ChameleonSph,
 )
 from .sheets import ExternalPotential, ExternalShear, MassSheet
+from .input import (
+    InputDeflections,
+    InputPotential,
+    GaussianRandomField,
+    LinearNDInterpolatorExt,
+)
diff --git a/autogalaxy/profiles/mass/input/__init__.py b/autogalaxy/profiles/mass/input/__init__.py
new file mode 100644
index 00000000..b69ee35f
--- /dev/null
+++ b/autogalaxy/profiles/mass/input/__init__.py
@@ -0,0 +1,4 @@
+from .interp import LinearNDInterpolatorExt
+from .input_deflections import InputDeflections
+from .input_potential import InputPotential
+from .gaussian_random_field import GaussianRandomField
diff --git a/autogalaxy/profiles/mass/input/gaussian_random_field.py b/autogalaxy/profiles/mass/input/gaussian_random_field.py
new file mode 100644
index 00000000..a6003c06
--- /dev/null
+++ b/autogalaxy/profiles/mass/input/gaussian_random_field.py
@@ -0,0 +1,127 @@
+from functools import cached_property
+
+import numpy as np
+
+import autoarray as aa
+
+from autogalaxy.profiles.mass.abstract.abstract import MassProfile
+from autogalaxy.profiles.mass.input.input_potential import InputPotential
+
+
+def gaussian_random_field_from(
+    shape_native, pixel_scale: float, power_amplitude: float, power_slope: float, seed: int
+) -> np.ndarray:
+    """
+    A real Gaussian random field realization with isotropic power-law power
+    spectrum P(k) = power_amplitude * k^(-power_slope).
+
+    The field is generated by filtering white Gaussian noise in Fourier space
+    (multiplying its transform by sqrt(P(k)) and inverse transforming), which
+    guarantees a real field whose power spectrum follows P(k); the k = 0 mode
+    is zeroed so the field has zero mean. Note this normalization convention
+    differs from the ``powerbox`` package used by the original
+    implementation — the spectrum's shape and seed-reproducibility are what
+    potential-correction validation relies on.
+
+    Parameters
+    ----------
+    shape_native
+        The 2D shape of the realization.
+    pixel_scale
+        The pixel size, setting the physical frequencies k.
+    power_amplitude
+        The amplitude of the power spectrum.
+    power_slope
+        The (positive) slope of the power-law spectrum, P(k) ~ k^-slope.
+    seed
+        The random seed, making the realization reproducible.
+    """
+    rng = np.random.default_rng(seed)
+    white_noise = rng.normal(size=shape_native)
+
+    ky = np.fft.fftfreq(shape_native[0], d=pixel_scale) * 2.0 * np.pi
+    kx = np.fft.fftfreq(shape_native[1], d=pixel_scale) * 2.0 * np.pi
+    k_grid = np.sqrt(ky[:, None] ** 2 + kx[None, :] ** 2)
+
+    power = np.zeros_like(k_grid)
+    nonzero = k_grid > 0
+    power[nonzero] = power_amplitude * k_grid[nonzero] ** (-power_slope)
+
+    field_ft = np.fft.fft2(white_noise) * np.sqrt(power)
+    return np.fft.ifft2(field_ft).real
+
+
+class GaussianRandomField(MassProfile):
+    def __init__(
+        self,
+        mask: aa.Mask2D,
+        power_amplitude: float = 1.0,
+        power_slope: float = 1.0,
+        seed: int = 1,
+    ):
+        """
+        A mass profile whose lensing potential is a Gaussian random field
+        realization with power spectrum P(k) = power_amplitude * k^(-power_slope),
+        used to simulate extended perturbations of a smooth lens-mass model
+        (e.g. for validating potential-correction reconstructions).
+
+        The realization is evaluated on the unmasked pixels of the input mask
+        and wrapped in an ``InputPotential``, from which the deflection angles
+        and convergence are derived via the mask's sparse derivative
+        operators.
+
+        Ported from the ``potential_correction`` package of Cao et al. 2025
+        (https://github.com/caoxiaoyue/lensing_potential_correction). If you
+        use this profile in your research, please cite Cao et al. 2025;
+        citation materials are provided at
+        https://github.com/caoxiaoyue/potential_correction_paper. The
+        realization here uses a plain-numpy Fourier filter rather than the
+        original's ``powerbox`` dependency (see
+        ``gaussian_random_field_from``).
+
+        Parameters
+        ----------
+        mask
+            The cleaned 2D mask (an ``aa.Mask2D`` carrying the pixel scale)
+            on whose unmasked pixels the potential is defined (see
+            ``aa.util.derivative.cleaned_mask_from``).
+        power_amplitude
+            The amplitude of the potential's power spectrum.
+        power_slope
+            The (positive) slope of the power-law spectrum, P(k) ~ k^-slope.
+        seed
+            The random seed, making the realization reproducible.
+        """
+        self.mask = mask
+        self.power_amplitude = power_amplitude
+        self.power_slope = power_slope
+        self.seed = seed
+        super().__init__()
+
+    @cached_property
+    def lensing_potential_native(self) -> np.ndarray:
+        return gaussian_random_field_from(
+            shape_native=self.mask.shape_native,
+            pixel_scale=self.mask.pixel_scale,
+            power_amplitude=self.power_amplitude,
+            power_slope=self.power_slope,
+            seed=self.seed,
+        )
+
+    @cached_property
+    def input_potential(self) -> InputPotential:
+        grid = aa.Grid2D.from_mask(mask=self.mask)
+        return InputPotential(
+            lensing_potential=self.lensing_potential_native[~np.asarray(self.mask)],
+            image_plane_grid=np.asarray(grid),
+            mask=self.mask,
+        )
+
+    def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        return self.input_potential.convergence_2d_from(grid=grid, xp=xp, **kwargs)
+
+    def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        return self.input_potential.potential_2d_from(grid=grid, xp=xp, **kwargs)
+
+    def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        return self.input_potential.deflections_yx_2d_from(grid=grid, xp=xp, **kwargs)
diff --git a/autogalaxy/profiles/mass/input/input_deflections.py b/autogalaxy/profiles/mass/input/input_deflections.py
new file mode 100644
index 00000000..703d58e0
--- /dev/null
+++ b/autogalaxy/profiles/mass/input/input_deflections.py
@@ -0,0 +1,116 @@
+from typing import Optional
+
+import numpy as np
+from scipy.sparse import spmatrix
+from scipy.spatial import Delaunay
+
+import autoarray as aa
+from autoarray.operators import derivative_util
+
+from autogalaxy.profiles.mass.abstract.abstract import MassProfile
+from autogalaxy.profiles.mass.input.interp import LinearNDInterpolatorExt
+
+
+class InputDeflections(MassProfile):
+    def __init__(
+        self,
+        deflections_y: np.ndarray,
+        deflections_x: np.ndarray,
+        image_plane_grid: aa.type.Grid2DLike,
+        mask: aa.type.Mask2D,
+        Hy: Optional[spmatrix] = None,
+        Hx: Optional[spmatrix] = None,
+    ):
+        """
+        A pixelized mass model defined by known deflection angles on the
+        unmasked pixels of an image-plane grid (e.g. from a previous lens
+        model, a particle simulation, or a potential-correction
+        reconstruction).
+
+        Deflections at arbitrary positions are evaluated by linear
+        interpolation over a Delaunay triangulation of the unmasked pixels
+        (nearest-neighbour fallback outside the convex hull). The convergence
+        is derived from the input deflections via the sparse first-derivative
+        operators of the mask, kappa = 0.5 * (dalpha_y/dy + dalpha_x/dx). The
+        lensing potential is not derivable from deflections alone and returns
+        zeros.
+
+        Ported from the ``potential_correction`` package of Cao et al. 2025
+        (https://github.com/caoxiaoyue/lensing_potential_correction). If you
+        use this profile in your research, please cite Cao et al. 2025;
+        citation materials are provided at
+        https://github.com/caoxiaoyue/potential_correction_paper.
+
+        Parameters
+        ----------
+        deflections_y
+            The 1D (slim) array of the y components of the deflection angles
+            on the unmasked pixels.
+        deflections_x
+            The 1D (slim) array of the x components of the deflection angles
+            on the unmasked pixels.
+        image_plane_grid
+            The [n_unmasked, 2] (y, x) grid of the unmasked pixels the
+            deflection angles are defined on.
+        mask
+            The cleaned 2D mask defining the unmasked pixels (see
+            ``aa.util.derivative.cleaned_mask_from``); its ``pixel_scale``
+            sets the finite-difference step of the derived convergence.
+        Hy
+            The sparse first-derivative operator along y of the mask; built
+            from the mask if not input.
+        Hx
+            The sparse first-derivative operator along x of the mask; built
+            from the mask if not input.
+        """
+        super().__init__()
+
+        self.deflections_y = np.asarray(deflections_y)
+        self.deflections_x = np.asarray(deflections_x)
+        self.image_plane_grid = np.asarray(image_plane_grid)
+        self.mask = mask
+        self.Hy = Hy
+        self.Hx = Hx
+
+        self._build_interpolators()
+
+    def _build_interpolators(self):
+        if self.Hy is None or self.Hx is None:
+            self.Hy, self.Hx = derivative_util.derivative_1st_operators_from(
+                np.asarray(self.mask), pixel_scale=self.mask.pixel_scale
+            )
+
+        self.convergence_slim = (
+            self.Hy @ self.deflections_y + self.Hx @ self.deflections_x
+        ) * 0.5
+
+        self.tri = Delaunay(np.fliplr(self.image_plane_grid))
+        self.interp_defl_y = LinearNDInterpolatorExt(self.tri, self.deflections_y)
+        self.interp_defl_x = LinearNDInterpolatorExt(self.tri, self.deflections_x)
+        self.interp_kappa = LinearNDInterpolatorExt(self.tri, self.convergence_slim)
+
+    @aa.decorators.to_array
+    def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        grid = np.asarray(grid)
+        return self.interp_kappa(grid[:, 1], grid[:, 0])
+
+    @aa.decorators.to_array
+    def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        return np.zeros(shape=np.asarray(grid).shape[0])
+
+    @aa.decorators.to_vector_yx
+    def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        """
+        Calculate the deflection angles at a given set of arc-second gridded
+        coordinates, by interpolating the input deflection angles.
+
+        Parameters
+        ----------
+        grid
+            The grid of (y,x) arc-second coordinates the deflection angles
+            are computed on.
+        """
+        grid = np.asarray(grid)
+        deflections_y = self.interp_defl_y(grid[:, 1], grid[:, 0])
+        deflections_x = self.interp_defl_x(grid[:, 1], grid[:, 0])
+        return np.stack((deflections_y, deflections_x), axis=-1)
diff --git a/autogalaxy/profiles/mass/input/input_potential.py b/autogalaxy/profiles/mass/input/input_potential.py
new file mode 100644
index 00000000..c8fcfbea
--- /dev/null
+++ b/autogalaxy/profiles/mass/input/input_potential.py
@@ -0,0 +1,130 @@
+from typing import Optional
+
+import numpy as np
+from scipy.sparse import spmatrix
+from scipy.spatial import Delaunay
+
+import autoarray as aa
+from autoarray.operators import derivative_util
+
+from autogalaxy.profiles.mass.abstract.abstract import MassProfile
+from autogalaxy.profiles.mass.input.interp import LinearNDInterpolatorExt
+
+
+class InputPotential(MassProfile):
+    def __init__(
+        self,
+        lensing_potential: np.ndarray,
+        image_plane_grid: aa.type.Grid2DLike,
+        mask: aa.type.Mask2D,
+        Hy: Optional[spmatrix] = None,
+        Hx: Optional[spmatrix] = None,
+        Hyy: Optional[spmatrix] = None,
+        Hxx: Optional[spmatrix] = None,
+    ):
+        """
+        A pixelized mass model defined by known lensing-potential values on
+        the unmasked pixels of an image-plane grid (e.g. a Gaussian random
+        field realization, or a potential-correction reconstruction).
+
+        The deflection angles are derived from the input potential via the
+        sparse first-derivative operators of the mask (alpha = grad psi) and
+        the convergence via the second-derivative operators
+        (kappa = 0.5 * laplacian psi). The potential, deflections and
+        convergence at arbitrary positions are evaluated by linear
+        interpolation over a Delaunay triangulation of the unmasked pixels
+        (nearest-neighbour fallback outside the convex hull).
+
+        Ported from the ``potential_correction`` package of Cao et al. 2025
+        (https://github.com/caoxiaoyue/lensing_potential_correction). If you
+        use this profile in your research, please cite Cao et al. 2025;
+        citation materials are provided at
+        https://github.com/caoxiaoyue/potential_correction_paper.
+
+        Parameters
+        ----------
+        lensing_potential
+            The 1D (slim) array of the lensing potential on the unmasked
+            pixels.
+        image_plane_grid
+            The [n_unmasked, 2] (y, x) grid of the unmasked pixels the
+            potential is defined on.
+        mask
+            The cleaned 2D mask defining the unmasked pixels (see
+            ``aa.util.derivative.cleaned_mask_from``); its ``pixel_scale``
+            sets the finite-difference step of the derived deflections and
+            convergence.
+        Hy
+            The sparse first-derivative operator along y of the mask; built
+            from the mask if not input.
+        Hx
+            The sparse first-derivative operator along x of the mask; built
+            from the mask if not input.
+        Hyy
+            The sparse second-derivative operator along y of the mask; built
+            from the mask if not input.
+        Hxx
+            The sparse second-derivative operator along x of the mask; built
+            from the mask if not input.
+        """
+        super().__init__()
+
+        self.lensing_potential = np.asarray(lensing_potential)
+        self.image_plane_grid = np.asarray(image_plane_grid)
+        self.mask = mask
+        self.Hy = Hy
+        self.Hx = Hx
+        self.Hyy = Hyy
+        self.Hxx = Hxx
+
+        self._build_interpolators()
+
+    def _build_interpolators(self):
+        if self.Hy is None or self.Hx is None:
+            self.Hy, self.Hx = derivative_util.derivative_1st_operators_from(
+                np.asarray(self.mask), pixel_scale=self.mask.pixel_scale
+            )
+        if self.Hyy is None or self.Hxx is None:
+            self.Hyy, self.Hxx = derivative_util.derivative_2nd_operators_from(
+                np.asarray(self.mask), pixel_scale=self.mask.pixel_scale
+            )
+
+        self.deflections_y = self.Hy @ self.lensing_potential
+        self.deflections_x = self.Hx @ self.lensing_potential
+        self.convergence_slim = (
+            self.Hyy @ self.lensing_potential + self.Hxx @ self.lensing_potential
+        ) * 0.5
+
+        self.tri = Delaunay(np.fliplr(self.image_plane_grid))
+        self.interp_psi = LinearNDInterpolatorExt(self.tri, self.lensing_potential)
+        self.interp_defl_y = LinearNDInterpolatorExt(self.tri, self.deflections_y)
+        self.interp_defl_x = LinearNDInterpolatorExt(self.tri, self.deflections_x)
+        self.interp_kappa = LinearNDInterpolatorExt(self.tri, self.convergence_slim)
+
+    @aa.decorators.to_array
+    def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        grid = np.asarray(grid)
+        return self.interp_kappa(grid[:, 1], grid[:, 0])
+
+    @aa.decorators.to_array
+    def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        grid = np.asarray(grid)
+        return self.interp_psi(grid[:, 1], grid[:, 0])
+
+    @aa.decorators.to_vector_yx
+    def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        """
+        Calculate the deflection angles at a given set of arc-second gridded
+        coordinates, by interpolating the deflections derived from the input
+        lensing potential.
+
+        Parameters
+        ----------
+        grid
+            The grid of (y,x) arc-second coordinates the deflection angles
+            are computed on.
+        """
+        grid = np.asarray(grid)
+        deflections_y = self.interp_defl_y(grid[:, 1], grid[:, 0])
+        deflections_x = self.interp_defl_x(grid[:, 1], grid[:, 0])
+        return np.stack((deflections_y, deflections_x), axis=-1)
diff --git a/autogalaxy/profiles/mass/input/interp.py b/autogalaxy/profiles/mass/input/interp.py
new file mode 100644
index 00000000..16d8b260
--- /dev/null
+++ b/autogalaxy/profiles/mass/input/interp.py
@@ -0,0 +1,34 @@
+import numpy as np
+from scipy.interpolate import LinearNDInterpolator
+from scipy.interpolate import NearestNDInterpolator
+
+
+class LinearNDInterpolatorExt:
+    def __init__(self, points, values):
+        """
+        Linear interpolation over a Delaunay triangulation of scattered 2D
+        points, falling back to nearest-neighbour interpolation outside the
+        convex hull so extrapolated values are never NaN.
+
+        Ported from the ``potential_correction`` package of Cao et al. 2025
+        (https://github.com/caoxiaoyue/lensing_potential_correction; cite via
+        https://github.com/caoxiaoyue/potential_correction_paper).
+
+        Parameters
+        ----------
+        points
+            The (x, y) coordinates the values are defined at, as an
+            [n_points, 2] array or a pre-built ``scipy.spatial.Delaunay``
+            triangulation of them.
+        values
+            The values interpolated.
+        """
+        self.funcinterp = LinearNDInterpolator(points, values)
+        self.funcnearest = NearestNDInterpolator(points, values)
+
+    def __call__(self, *args):
+        z = self.funcinterp(*args)
+        chk = np.isnan(z)
+        if chk.any():
+            return np.where(chk, self.funcnearest(*args), z)
+        return z
diff --git a/test_autogalaxy/config/notation.yaml b/test_autogalaxy/config/notation.yaml
index f88e205d..5d16619e 100644
--- a/test_autogalaxy/config/notation.yaml
+++ b/test_autogalaxy/config/notation.yaml
@@ -1,82 +1,82 @@
-label:
-  label:
-    alpha: \alpha
-    angle_binary: \theta
-    beta: \beta
-    break_radius: \theta_{\rm B}
-    centre_0: y
-    centre_1: x
-    coefficient: \lambda
-    core_radius: C_{\rm r}
-    core_radius_0: C_{rm r0}
-    core_radius_1: C_{\rm r1}
-    effective_radius: R_{\rm eff}
-    einstein_radius: \theta_{\rm Ein}
-    ell_comps_0: \epsilon_{\rm 1}
-    ell_comps_1: \epsilon_{\rm 2}
-    multipole_comps_0: M_{\rm 1}
-    multipole_comps_1: M_{\rm 2}
-    flux: F
-    gamma: \gamma
-    gamma_1: \gamma
-    gamma_2: \gamma
-    inner_coefficient: \lambda_{\rm 1}
-    inner_slope: t_{\rm 1}
-    intensity: I_{\rm b}
-    kappa: \kappa
-    kappa_s: \kappa_{\rm s}
-    log10m_vir: log_{\rm 10}(m_{vir})
-    m: m
-    mass: M
-    mass_at_200: M_{\rm 200}
-    mass_ratio: M_{\rm ratio}
-    mass_to_light_gradient: \Gamma
-    mass_to_light_ratio: \Psi
-    mass_to_light_ratio_base: \Psi_{\rm base}
-    mass_to_light_radius: R_{\rm ref}
-    noise_factor: \omega_{\rm 1}
-    noise_power: \omega{\rm 2}
-    noise_scale: \sigma_{\rm 1}
-    normalization_scale: n
-    outer_coefficient: \lambda_{\rm 2}
-    outer_slope: t_{\rm 2}
-    overdens: \Delta_{\rm vir}
-    pixels: N_{\rm pix}
-    radius_break: R_{\rm b}
-    redshift: z
-    redshift_object: z_{\rm obj}
-    redshift_source: z_{\rm src}
-    scale_radius: R_{\rm s}
-    scatter: \sigma
-    separation: s
-    sersic_index: n
-    shape_0: y_{\rm pix}
-    shape_1: x_{\rm pix}
-    sigma: \sigma
-    signal_scale: V
-    sky_scale: \sigma_{\rm 0}
-    slope: \gamma
-    truncation_radius: R_{\rm t}
-    weight_floor: W_{\rm f}
-    weight_power: W_{\rm p}
-  superscript:
-    externalshear: ext
-    pixelization: pix
-    point: point
-    redshift: ''
-    regularization: reg
-label_format:
-  format:
-    angular_diameter_distance_to_earth: '{:.2f}'
-    concentration: '{:.2f}'
-    einstein_mass: '{:.4e}'
-    einstein_radius: '{:.2f}'
-    kpc_per_arcsec: '{:.2f}'
-    luminosity: '{:.4e}'
-    m: '{:.1f}'
-    mass: '{:.4e}'
-    mass_at_truncation_radius: '{:.4e}'
-    radius: '{:.2f}'
-    redshift: '{:.2f}'
-    rho: '{:.2f}'
-    sersic_luminosity: '{:.4e}'
+label:
+  label:
+    alpha: \alpha
+    angle_binary: \theta
+    beta: \beta
+    break_radius: \theta_{\rm B}
+    centre_0: y
+    centre_1: x
+    coefficient: \lambda
+    core_radius: C_{\rm r}
+    core_radius_0: C_{rm r0}
+    core_radius_1: C_{\rm r1}
+    effective_radius: R_{\rm eff}
+    einstein_radius: \theta_{\rm Ein}
+    ell_comps_0: \epsilon_{\rm 1}
+    ell_comps_1: \epsilon_{\rm 2}
+    multipole_comps_0: M_{\rm 1}
+    multipole_comps_1: M_{\rm 2}
+    flux: F
+    gamma: \gamma
+    gamma_1: \gamma
+    gamma_2: \gamma
+    inner_coefficient: \lambda_{\rm 1}
+    inner_slope: t_{\rm 1}
+    intensity: I_{\rm b}
+    kappa: \kappa
+    kappa_s: \kappa_{\rm s}
+    log10m_vir: log_{\rm 10}(m_{vir})
+    m: m
+    mass: M
+    mass_at_200: M_{\rm 200}
+    mass_ratio: M_{\rm ratio}
+    mass_to_light_gradient: \Gamma
+    mass_to_light_ratio: \Psi
+    mass_to_light_ratio_base: \Psi_{\rm base}
+    mass_to_light_radius: R_{\rm ref}
+    noise_factor: \omega_{\rm 1}
+    noise_power: \omega{\rm 2}
+    noise_scale: \sigma_{\rm 1}
+    normalization_scale: n
+    outer_coefficient: \lambda_{\rm 2}
+    outer_slope: t_{\rm 2}
+    overdens: \Delta_{\rm vir}
+    pixels: N_{\rm pix}
+    radius_break: R_{\rm b}
+    redshift: z
+    redshift_object: z_{\rm obj}
+    redshift_source: z_{\rm src}
+    scale_radius: R_{\rm s}
+    scatter: \sigma
+    separation: s
+    sersic_index: n
+    shape_0: y_{\rm pix}
+    shape_1: x_{\rm pix}
+    sigma: \sigma
+    signal_scale: V
+    sky_scale: \sigma_{\rm 0}
+    slope: \gamma
+    truncation_radius: R_{\rm t}
+    weight_floor: W_{\rm f}
+    weight_power: W_{\rm p}
+  superscript:
+    externalshear: ext
+    pixelization: pix
+    point: point
+    redshift: ''
+    regularization: reg
+label_format:
+  format:
+    angular_diameter_distance_to_earth: '{:.2f}'
+    concentration: '{:.2f}'
+    einstein_mass: '{:.4e}'
+    einstein_radius: '{:.2f}'
+    kpc_per_arcsec: '{:.2f}'
+    luminosity: '{:.4e}'
+    m: '{:.1f}'
+    mass: '{:.4e}'
+    mass_at_truncation_radius: '{:.4e}'
+    radius: '{:.2f}'
+    redshift: '{:.2f}'
+    rho: '{:.2f}'
+    sersic_luminosity: '{:.4e}'
diff --git a/test_autogalaxy/profiles/mass/input/__init__.py b/test_autogalaxy/profiles/mass/input/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/test_autogalaxy/profiles/mass/input/test_gaussian_random_field.py b/test_autogalaxy/profiles/mass/input/test_gaussian_random_field.py
new file mode 100644
index 00000000..f0dd1ba6
--- /dev/null
+++ b/test_autogalaxy/profiles/mass/input/test_gaussian_random_field.py
@@ -0,0 +1,81 @@
+import numpy as np
+import pytest
+
+import autoarray as aa
+import autogalaxy as ag
+from autogalaxy.profiles.mass.input.gaussian_random_field import (
+    gaussian_random_field_from,
+)
+
+
+def cleaned_circular_mask(shape=(16, 16), pixel_scale=0.5, radius=3.4):
+    mask = aa.Mask2D.circular(
+        shape_native=shape, pixel_scales=pixel_scale, radius=radius
+    )
+    cleaned, _ = aa.util.derivative.cleaned_mask_from(np.asarray(mask))
+    return aa.Mask2D(mask=cleaned, pixel_scales=pixel_scale)
+
+
+def test__realization_is_reproducible_and_zero_mean():
+    field_a = gaussian_random_field_from(
+        shape_native=(32, 32),
+        pixel_scale=0.5,
+        power_amplitude=1.0,
+        power_slope=2.0,
+        seed=3,
+    )
+    field_b = gaussian_random_field_from(
+        shape_native=(32, 32),
+        pixel_scale=0.5,
+        power_amplitude=1.0,
+        power_slope=2.0,
+        seed=3,
+    )
+    field_c = gaussian_random_field_from(
+        shape_native=(32, 32),
+        pixel_scale=0.5,
+        power_amplitude=1.0,
+        power_slope=2.0,
+        seed=4,
+    )
+
+    assert field_a == pytest.approx(field_b, abs=0.0)
+    assert not np.allclose(field_a, field_c)
+    assert field_a.mean() == pytest.approx(0.0, abs=1.0e-10)
+
+
+def test__amplitude_scales_field():
+    field_1 = gaussian_random_field_from(
+        shape_native=(32, 32),
+        pixel_scale=0.5,
+        power_amplitude=1.0,
+        power_slope=2.0,
+        seed=3,
+    )
+    field_4 = gaussian_random_field_from(
+        shape_native=(32, 32),
+        pixel_scale=0.5,
+        power_amplitude=4.0,
+        power_slope=2.0,
+        seed=3,
+    )
+
+    assert field_4 == pytest.approx(2.0 * field_1, abs=1.0e-10)
+
+
+def test__profile_potential_matches_realization_on_unmasked_pixels():
+    mask = cleaned_circular_mask()
+    grid = aa.Grid2D.from_mask(mask=mask)
+
+    profile = ag.mp.GaussianRandomField(
+        mask=mask, power_amplitude=1.0, power_slope=1.0, seed=2
+    )
+
+    potential = np.asarray(profile.potential_2d_from(grid=grid))
+    expected = profile.lensing_potential_native[~np.asarray(mask)]
+
+    assert potential == pytest.approx(expected, abs=1.0e-10)
+
+    deflections = np.asarray(profile.deflections_yx_2d_from(grid=grid))
+    assert deflections.shape == (expected.shape[0], 2)
+    assert np.isfinite(deflections).all()
diff --git a/test_autogalaxy/profiles/mass/input/test_input_deflections.py b/test_autogalaxy/profiles/mass/input/test_input_deflections.py
new file mode 100644
index 00000000..8328b7d6
--- /dev/null
+++ b/test_autogalaxy/profiles/mass/input/test_input_deflections.py
@@ -0,0 +1,80 @@
+import numpy as np
+import pytest
+
+import autoarray as aa
+import autogalaxy as ag
+
+
+def masked_setup(shape=(16, 16), pixel_scale=0.5, radius=3.4):
+    mask = aa.Mask2D.circular(
+        shape_native=shape, pixel_scales=pixel_scale, radius=radius
+    )
+    cleaned, _ = aa.util.derivative.cleaned_mask_from(np.asarray(mask))
+    mask = aa.Mask2D(mask=cleaned, pixel_scales=pixel_scale)
+    grid = aa.Grid2D.from_mask(mask=mask)
+    return mask, grid
+
+
+def test__linear_deflections__interp_exact_at_nodes_and_convergence_exact():
+    mask, grid = masked_setup()
+    grid_arr = np.asarray(grid)
+
+    # alpha_y = a y, alpha_x = b x (the deflections of psi = a y^2/2 + b x^2/2):
+    # kappa = 0.5 (a + b), exact for every first-difference scheme since the
+    # deflections are linear in the coordinates.
+    a, b = 1.5, 0.5
+    deflections_y = a * grid_arr[:, 0]
+    deflections_x = b * grid_arr[:, 1]
+
+    profile = ag.mp.InputDeflections(
+        deflections_y=deflections_y,
+        deflections_x=deflections_x,
+        image_plane_grid=grid_arr,
+        mask=mask,
+    )
+
+    deflections = np.asarray(profile.deflections_yx_2d_from(grid=grid))
+    assert deflections[:, 0] == pytest.approx(deflections_y, abs=1.0e-8)
+    assert deflections[:, 1] == pytest.approx(deflections_x, abs=1.0e-8)
+
+    convergence = np.asarray(profile.convergence_2d_from(grid=grid))
+    assert convergence == pytest.approx(0.5 * (a + b), abs=1.0e-8)
+
+
+def test__potential_returns_zeros():
+    mask, grid = masked_setup()
+    grid_arr = np.asarray(grid)
+
+    profile = ag.mp.InputDeflections(
+        deflections_y=np.ones(grid_arr.shape[0]),
+        deflections_x=np.ones(grid_arr.shape[0]),
+        image_plane_grid=grid_arr,
+        mask=mask,
+    )
+
+    potential = np.asarray(profile.potential_2d_from(grid=grid))
+    assert potential == pytest.approx(0.0, abs=1.0e-12)
+
+
+def test__interpolation_off_nodes_is_linear():
+    mask, grid = masked_setup()
+    grid_arr = np.asarray(grid)
+
+    a, b = 1.5, 0.5
+    profile = ag.mp.InputDeflections(
+        deflections_y=a * grid_arr[:, 0],
+        deflections_x=b * grid_arr[:, 1],
+        image_plane_grid=grid_arr,
+        mask=mask,
+    )
+
+    # linear interpolation reproduces a linear deflection field exactly at
+    # positions inside the convex hull, including off-node positions.
+    off_node = aa.Grid2DIrregular(values=[(0.13, -0.41), (1.07, 0.66)])
+    deflections = np.asarray(profile.deflections_yx_2d_from(grid=off_node))
+    assert deflections[:, 0] == pytest.approx(
+        a * np.asarray(off_node)[:, 0], abs=1.0e-8
+    )
+    assert deflections[:, 1] == pytest.approx(
+        b * np.asarray(off_node)[:, 1], abs=1.0e-8
+    )
diff --git a/test_autogalaxy/profiles/mass/input/test_input_potential.py b/test_autogalaxy/profiles/mass/input/test_input_potential.py
new file mode 100644
index 00000000..c382d4a0
--- /dev/null
+++ b/test_autogalaxy/profiles/mass/input/test_input_potential.py
@@ -0,0 +1,82 @@
+import numpy as np
+import pytest
+
+import autoarray as aa
+import autogalaxy as ag
+
+
+def masked_setup(shape=(16, 16), pixel_scale=0.5, radius=3.4):
+    mask = aa.Mask2D.circular(
+        shape_native=shape, pixel_scales=pixel_scale, radius=radius
+    )
+    cleaned, _ = aa.util.derivative.cleaned_mask_from(np.asarray(mask))
+    mask = aa.Mask2D(mask=cleaned, pixel_scales=pixel_scale)
+    grid = aa.Grid2D.from_mask(mask=mask)
+    return mask, grid
+
+
+def test__linear_potential__deflections_exact_and_convergence_zero():
+    mask, grid = masked_setup()
+    grid_arr = np.asarray(grid)
+
+    # psi = 2 y + 3 x: alpha = (2, 3) everywhere, kappa = 0. All the mask's
+    # finite-difference schemes are exact on a linear potential.
+    potential = 2.0 * grid_arr[:, 0] + 3.0 * grid_arr[:, 1]
+
+    profile = ag.mp.InputPotential(
+        lensing_potential=potential, image_plane_grid=grid_arr, mask=mask
+    )
+
+    deflections = np.asarray(profile.deflections_yx_2d_from(grid=grid))
+    assert deflections[:, 0] == pytest.approx(2.0, abs=1.0e-8)
+    assert deflections[:, 1] == pytest.approx(3.0, abs=1.0e-8)
+
+    convergence = np.asarray(profile.convergence_2d_from(grid=grid))
+    assert convergence == pytest.approx(0.0, abs=1.0e-8)
+
+
+def test__quadratic_potential__convergence_exact_and_potential_interp_at_nodes():
+    mask, grid = masked_setup()
+    grid_arr = np.asarray(grid)
+
+    # psi = a y^2 + b x^2: kappa = 0.5 (2a + 2b) = a + b. Second-difference
+    # schemes (central and one-sided) are exact on quadratics.
+    a, b = 0.75, 0.25
+    potential = a * grid_arr[:, 0] ** 2 + b * grid_arr[:, 1] ** 2
+
+    profile = ag.mp.InputPotential(
+        lensing_potential=potential, image_plane_grid=grid_arr, mask=mask
+    )
+
+    convergence = np.asarray(profile.convergence_2d_from(grid=grid))
+    assert convergence == pytest.approx(a + b, abs=1.0e-8)
+
+    # linear interpolation at the triangulation nodes returns the inputs
+    potential_interp = np.asarray(profile.potential_2d_from(grid=grid))
+    assert potential_interp == pytest.approx(potential, abs=1.0e-8)
+
+
+def test__operators_can_be_preloaded():
+    mask, grid = masked_setup()
+    grid_arr = np.asarray(grid)
+    potential = np.asarray(grid_arr[:, 0])
+
+    Hy, Hx = aa.util.derivative.derivative_1st_operators_from(
+        np.asarray(mask), pixel_scale=mask.pixel_scale
+    )
+    Hyy, Hxx = aa.util.derivative.derivative_2nd_operators_from(
+        np.asarray(mask), pixel_scale=mask.pixel_scale
+    )
+
+    profile = ag.mp.InputPotential(
+        lensing_potential=potential,
+        image_plane_grid=grid_arr,
+        mask=mask,
+        Hy=Hy,
+        Hx=Hx,
+        Hyy=Hyy,
+        Hxx=Hxx,
+    )
+
+    deflections = np.asarray(profile.deflections_yx_2d_from(grid=grid))
+    assert deflections[:, 0] == pytest.approx(1.0, abs=1.0e-8)