Merge pull request #366 from PyAutoLabs/claude/list-admin-prompts-Mid1U

docs(weak-lensing): document shear API across LensCalc, Isothermal, ExternalShear

Author: Jammy2211

autogalaxy/operate/lens_calc.py

@@ -519,33 +519,53 @@ def convergence_2d_via_hessian_from(self, grid, xp=np) -> aa.ArrayIrregular:
         return convergence
 
     def shear_yx_2d_via_hessian_from(self, grid, xp=np) -> ShearYX2DIrregular:
-        """
-        Returns the 2D (y,x) shear vectors of the lensing object, which are computed from the 2D deflection angle map
-        via the Hessian using the expressions (see equation 57 https://inspirehep.net/literature/419263):
+        r"""
+        Returns the 2D weak-lensing shear vector field of the lensing object, computed from the deflection-angle
+        Hessian (see equation 57 of https://inspirehep.net/literature/419263):
 
-        `shear_y = hessian_{1,0} =  hessian_{0,1} = hessian_yx = hessian_xy`
-        `shear_x = 0.5 * (hessian_{0,0} - hessian_{1,1}) = 0.5 * (hessian_xx - hessian_yy)`
+        .. math::
 
-        By going via the Hessian, the shear vectors can be calculated at any (y,x) coordinate, therefore using either a
-        2D uniform or irregular grid.
+            \gamma_1 = \tfrac{1}{2} (\partial^2 \psi / \partial x^2 - \partial^2 \psi / \partial y^2)
+                     = 0.5 \, (H_{xx} - H_{yy})
 
-        This calculation of the shear vectors is independent of analytic calculations defined within `MassProfile`
-        objects and can therefore be used as a cross-check.
+            \gamma_2 = \partial^2 \psi / \partial x \partial y
+                     = H_{xy} = H_{yx}
 
-        The result is returned as a `ShearYX2D` dats structure, which has shape [total_shear_vectors, 2], where
-        entries for [:,0] are the gamma_2 values and entries for [:,1] are the gamma_1 values.
+        where :math:`\psi(\theta)` is the lensing potential and :math:`H_{ij}` are its second partial derivatives,
+        i.e. the components of the Hessian returned by ``hessian_from``.  The two components together encode the
+        full distortion tensor of a background source: :math:`\gamma_1` is the stretch along the x/y axes and
+        :math:`\gamma_2` is the stretch along the diagonals.
 
-        Note therefore that this convention means the FIRST entries in the array are the gamma_2 values and the SECOND
-        entries are the gamma_1 values.
+        Because the Hessian is computed numerically (Richardson-extrapolated finite differences in NumPy, or
+        ``jax.jacfwd`` in JAX), this routine works for **any** ``MassProfile``, ``Galaxy``, ``Galaxies`` or
+        ``Tracer`` that exposes ``deflections_yx_2d_from`` — no per-profile analytic formula is required.  It can
+        therefore be used as a numerical cross-check of analytic shear methods such as
+        ``Isothermal.shear_yx_2d_from`` (see ``test_isothermal.py``).  It also accepts both uniform ``Grid2D`` and
+        irregular ``Grid2DIrregular`` grids, which makes it the natural choice for simulating a weak-lensing
+        shear field on either a regular pixel grid or at the discrete positions of background source galaxies.
+
+        Convention
+        ----------
+        The result is returned as a ``ShearYX2DIrregular`` data structure of shape ``[total_shear_vectors, 2]``,
+        where:
+
+        - ``[:, 0]`` are the :math:`\gamma_2` values
+        - ``[:, 1]`` are the :math:`\gamma_1` values
+
+        Note therefore that the FIRST column is :math:`\gamma_2` and the SECOND column is :math:`\gamma_1`.  This
+        ordering matches the ``[y, x]`` ordering used elsewhere in the project for vector fields, since
+        :math:`\gamma_2` plays the role of the "y" component and :math:`\gamma_1` the "x" component when the shear
+        field is treated as a 2D pseudo-vector.
 
         Parameters
         ----------
-        grids
+        grid
             The 2D grid of (y,x) arc-second coordinates the deflection angles and Hessian are computed on.
         xp
-            The array module to use for the computation (e.g. `numpy` or `jax.numpy`). Passed through to
-            `hessian_from`. When `xp` is not `numpy` (e.g. inside a `jax.jit` trace) the result is returned
-            as a raw array of shape `(N, 2)` rather than a `ShearYX2DIrregular` wrapper.
+            The array module to use for the computation (e.g. ``numpy`` or ``jax.numpy``).  Passed through to
+            ``hessian_from``.  When ``xp`` is not ``numpy`` (e.g. inside a ``jax.jit`` trace) the result is
+            returned as a raw array of shape ``(N, 2)`` rather than a ``ShearYX2DIrregular`` wrapper, because
+            ``ShearYX2DIrregular`` is not a registered JAX pytree.
         """
 
         hessian_yy, hessian_xy, hessian_yx, hessian_xx = self.hessian_from(

autogalaxy/profiles/mass/sheets/external_shear.py

@@ -9,27 +9,51 @@
 
 class ExternalShear(MassProfile):
     def __init__(self, gamma_1: float = 0.0, gamma_2: float = 0.0):
-        """
-        An `ExternalShear` term, to model the line-of-sight contribution of other galaxies / satellites.
+        r"""
+        Constant external shear term used in strong-lens mass models.
+
+        ``ExternalShear`` represents the line-of-sight contribution to the lensing potential from mass that is
+        not part of the primary lens — typically nearby group/cluster members, large-scale structure, or
+        unmodelled satellites.  Because this contribution is approximately uniform across the small angular
+        extent of a strong-lens system, it is parameterised as a constant shear with two components,
+        :math:`\gamma_1` and :math:`\gamma_2`, in the same convention as
+        ``LensCalc.shear_yx_2d_via_hessian_from`` and ``ShearYX2DIrregular``:
+
+        - :math:`\gamma_1` produces stretching along the x/y axes
+        - :math:`\gamma_2` produces stretching along the diagonals
+
+        The associated shear *magnitude* and *position angle* (degrees, anticlockwise from the +x axis) are:
 
-        The shear angle is defined in the direction of stretching of the image. Therefore, if an object located \
-        outside the lens is responsible for the shear, it will be offset 90 degrees from the value of angle.
+        .. math::
+
+            |\gamma| = \sqrt{\gamma_1^2 + \gamma_2^2}, \qquad
+            \phi = \tfrac{1}{2} \, \mathrm{arctan2}(\gamma_2, \gamma_1).
+
+        Note that the shear *position angle* lies in the direction of image stretching.  An external mass
+        located in the *direction of compression* therefore appears at an angle offset by 90 degrees from
+        :math:`\phi`.
 
         Parameters
         ----------
-        gamma
+        gamma_1
+            The :math:`\gamma_1` shear component.
+        gamma_2
+            The :math:`\gamma_2` shear component.
         """
 
         super().__init__(centre=(0.0, 0.0), ell_comps=(0.0, 0.0))
         self.gamma_1 = gamma_1
         self.gamma_2 = gamma_2
 
     def magnitude(self, xp=np):
+        r"""Returns the shear magnitude :math:`|\gamma| = \sqrt{\gamma_1^2 + \gamma_2^2}`."""
         return convert.shear_magnitude_from(
             gamma_1=self.gamma_1, gamma_2=self.gamma_2, xp=xp
         )
 
     def angle(self, xp=np):
+        r"""Returns the shear position angle :math:`\phi = \tfrac{1}{2}\,\mathrm{arctan2}(\gamma_2, \gamma_1)`
+        in degrees, in the [0, 180) convention used elsewhere in the package."""
         return convert.shear_angle_from(
             gamma_1=self.gamma_1, gamma_2=self.gamma_2, xp=xp
         )
@@ -42,13 +66,25 @@ def average_convergence_of_1_radius(self):
 
     @aa.decorators.to_array
     def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        """A pure shear term has zero convergence at every grid point."""
         return xp.zeros(shape=grid.shape[0])
 
     @aa.decorators.to_array
     def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
+        r"""
+        Returns the lensing potential of the constant external shear, given by:
+
+        .. math::
+
+            \psi(\theta) = -\tfrac{1}{2} |\gamma| \, r^2 \, \cos\!\big(2\,(\varphi - \phi_\gamma)\big)
+
+        where :math:`r, \varphi` are the polar coordinates of ``grid`` and :math:`\phi_\gamma` is the shear
+        position angle (offset by ``-90 deg`` to remain consistent with the deflection-angle convention used
+        elsewhere in PyAutoLens).
+        """
         shear_angle = (
             self.angle(xp) - 90
-        )  ##to be onsistent with autolens deflection angle calculation
+        )  # offset by -90 deg to match the deflection-angle convention used elsewhere
         phig = xp.deg2rad(shear_angle)
         shear_amp = self.magnitude(xp=xp)
         phicoord = xp.arctan2(grid.array[:, 0], grid.array[:, 1])
@@ -59,14 +95,28 @@ def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
     @aa.decorators.to_vector_yx
     @aa.decorators.transform(rotate_back=True)
     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.
+        r"""
+        Returns the deflection angles of the constant external shear at each ``(y, x)`` arc-second
+        coordinate.
+
+        In the profile's rotated reference frame (where the shear position angle is aligned with the x-axis)
+        the deflection reduces to:
+
+        .. math::
+
+            \alpha_y(y, x) = -|\gamma| \, y, \qquad \alpha_x(y, x) = +|\gamma| \, x.
+
+        The ``@transform(rotate_back=True)`` decorator rotates the input grid into this aligned frame and
+        rotates the resulting deflection vectors back into the original frame.  Because the shear is a spin-2
+        field, the *position angle* :math:`\phi` of the shear is encoded in the rotation of the grid and not
+        as an extra factor in the deflection formula above.
 
         Parameters
         ----------
         grid
             The grid of (y,x) arc-second coordinates the deflection angles are computed on.
-
+        xp
+            The array module (``numpy`` or ``jax.numpy``).
         """
         deflection_y = -xp.multiply(self.magnitude(xp=xp), grid.array[:, 0])
         deflection_x = xp.multiply(self.magnitude(xp=xp), grid.array[:, 1])

autogalaxy/profiles/mass/total/isothermal.py

@@ -110,20 +110,45 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
     @aa.decorators.to_vector_yx
     @aa.decorators.transform
     def shear_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
-        """
-        Calculate the (gamma_y, gamma_x) shear vector field on a grid of (y,x) arc-second coordinates.
+        r"""
+        Returns the analytic 2D weak-lensing shear vector field :math:`(\gamma_2, \gamma_1)` of the elliptical
+        isothermal mass distribution on a grid of ``(y, x)`` arc-second coordinates.
+
+        For an axis-aligned isothermal profile centred on the origin the shear components reduce to:
+
+        .. math::
 
-        The result is returned as a `ShearYX2D` dats structure, which has shape [total_shear_vectors, 2], where
-        entries for [:,0] are the gamma_2 values and entries for [:,1] are the gamma_1 values.
+            \gamma_1 = -\kappa(\theta) \, \frac{x^2 - y^2}{x^2 + y^2}
 
-        Note therefore that this convention means the FIRST entries in the array are the gamma_2 values and the SECOND
-        entries are the gamma_1 values.
+            \gamma_2 = -2 \, \kappa(\theta) \, \frac{x \, y}{x^2 + y^2}
+
+        where :math:`\kappa(\theta)` is the convergence at the rotated grid coordinate.  After evaluation in the
+        profile's reference frame the shear vector field is rotated back into the original frame using the
+        ``2 * angle`` rotation appropriate for a spin-2 quantity (the shear transforms as a spin-2 field, so a
+        coordinate rotation by ``angle`` rotates the components by ``2 * angle``).
+
+        This analytic path is mathematically equivalent to ``LensCalc.shear_yx_2d_via_hessian_from``, which
+        derives the same shear from finite-difference (or JAX) derivatives of ``deflections_yx_2d_from``.  The
+        cross-check is exercised in
+        ``test_autogalaxy/profiles/mass/total/test_isothermal.py::test__shear_yx_2d_from__matches_via_hessian``.
+
+        Convention
+        ----------
+        The result is returned as a vector-field with shape ``[total_shear_vectors, 2]`` where:
+
+        - ``[:, 0]`` are the :math:`\gamma_2` values
+        - ``[:, 1]`` are the :math:`\gamma_1` values
+
+        i.e. the FIRST column is :math:`\gamma_2` and the SECOND column is :math:`\gamma_1`.  This ordering
+        matches the convention used by ``ShearYX2D`` / ``ShearYX2DIrregular`` and
+        ``LensCalc.shear_yx_2d_via_hessian_from``.
 
         Parameters
         ----------
         grid
-            The grid of (y,x) arc-second coordinates the deflection angles are computed on.
-
+            The grid of (y,x) arc-second coordinates the shear vectors are computed on.
+        xp
+            The array module (``numpy`` or ``jax.numpy``).
         """
 
         convergence = self.convergence_2d_from(grid=grid, xp=xp, **kwargs)