From d72181aaa7a6114f6c5e0a2fbc12cb6374de012e Mon Sep 17 00:00:00 2001
From: robjmcgibbon <robjmcgibbon@gmail.com>
Date: Wed, 22 Apr 2026 17:00:27 +0100
Subject: [PATCH 1/7] Implement functions

---
 .../particle_selection/aperture_properties.py | 120 +++++++++++++++++-
 SOAP/property_table.py                        |  26 ++++
 2 files changed, 145 insertions(+), 1 deletion(-)

diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index a14c9be1..e1058db3 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -137,6 +137,7 @@ def MetalFracStar(self):
 import numpy as np
 from numpy.typing import NDArray
 from typing import Dict, List, Tuple
+import healpy as hp
 import unyt
 
 from .halo_properties import HaloProperty, SearchRadiusTooSmallError
@@ -3586,7 +3587,7 @@ def stellar_inertia_tensor(self, **kwargs) -> unyt.unyt_array:
         # aperture
         mass = self.get_dataset("PartType4/Masses")[self.star_mask_all]
         position = (
-            self.get_dataset("PartType4/Coordinates")[self.star_mask_all] - self.centre
+            self.get_dataset("PartType4/Coordinates")[self.star_mask_all] - self.centre - self.shrinking_sphere_centre
         )
 
         return get_inertia_tensor_mass_weighted(
@@ -3698,6 +3699,121 @@ def StellarInertiaTensorReducedNoniterativeLuminosityWeighted(
             reduced=True, max_iterations=1
         )
 
+    @lazy_property
+    def shrinking_sphere_centre(
+            self,
+            min_particles=200,
+            shrink_factor=0.83,
+            max_iter = 256
+        ) -> unyt.unyt_array:
+        """
+        Estimate the galaxy center (center of mass) with iterative shrinking aperture.
+
+        Procedure:
+        1) Initialize center as the COM of all input particles.
+        2) Initialize radius as max(aperture_radius, farthest particle distance).
+        3) Recompute COM inside current aperture, recenter, then shrink radius by shrink_factor.
+        4) Stop when enclosed particle count is < min_particles, and return the last valid center.
+
+        Parameters
+        ----------
+        min_particles : int, default=200
+            Stop iteration when enclosed particle count is below this threshold.
+        shrink_factor : float, default=0.83
+            Radius scaling factor applied at each iteration (0 < shrink_factor < 1).
+        max_iter : int, default=256
+            Maximum number of iterations
+        """
+
+        if self.Mbaryons == 0:
+            return None
+
+        centre = (self.baryon_mass_fraction[:, None] * self.pos_baryons).sum(axis=0)
+        dist = np.linalg.norm(self.pos_baryons - centre, axis=1)
+        radius = max(self.aperture_radius, dist.max())
+
+        for i in range(max_iter):
+            mask = dist <= radius
+            n_in = np.sum(mask)
+            if n_in < min_particles:
+                break
+
+            centre = (self.baryon_mass_fraction[mask, None] * self.pos_baryons[mask]).sum(axis=0)
+            dist = np.linalg.norm(self.pos_baryons - centre, axis=1)
+            radius *= shrink_factor
+
+        return centre
+
+    @lazy_property
+    def ShrinkingSphereCentre(self) -> unyt.unyt_array:
+        """
+        Centre computed by applying shrinking spheres method to baryons
+        """
+        return (self.shrinking_sphere_centre + self.centre) % self.boxsize
+
+    @lazy_property
+    def StellarAsymmetry(self, npix=12):
+        """
+        Compute stellar asymmetry following https://arxiv.org/abs/1805.03210
+
+        TODO: Is equation (3) incorrect?
+
+        Parameters
+        ----------
+        npix : int, default=12
+            Number of equal-area angular regions (HEALPix pixels) for directional
+            partitioning. Must satisfy npix = 12 * nside^2 where nside is a
+            power of 2 (e.g. 12, 48, 192, 768, ...).
+
+        """
+
+        if self.Mstar == 0:
+            return None
+
+        # Check we are using a valid value for npix
+        nside = int(round(np.sqrt(npix // 12)))
+        assert npix == 12 * nside ** 2
+        assert nside.bit_count() == 1
+
+        # Centre using the shrinking sphere
+        pos = self.pos_star - self.shrinking_sphere_centre
+        r = np.linalg.norm(pos, axis=1)
+
+        # Remove particles close to the centre, they are symmetric
+        mask = r.to_value('kpc') < 0.1
+        if np.sum(mask):
+            pos = pos[np.logical_not(mask)]
+            r = r[np.logical_not(mask)]
+            mass_star = self.mass_star[np.logical_not(mask)]
+            if r.shape[0] == 0:
+                return np.float32(0)
+        else:
+            mass_star = self.mass_star
+
+        # Compute the mass in each pixel
+        vecs = (pos / r[:, None]).value
+        idx = hp.vec2pix(nside, vecs[:, 0], vecs[:, 1], vecs[:, 2])
+        # np.bincount will not return a unyt array
+        region_mass_msun = np.bincount(
+            idx,
+            weights=mass_star.to_value('Msun'),
+            minlength=npix,
+        )
+
+        # Create antipodal mapping
+        # Find the center vector of every pixel
+        vecs = hp.pix2vec(nside, np.arange(npix))
+        # Negate vectors to find antipodal points
+        anti_vecs = -np.array(vecs)
+        # Map those points back to pixel IDs
+        anti_indices = hp.vec2pix(nside, anti_vecs[0], anti_vecs[1], anti_vecs[2])
+
+        # Calculate asymmetry
+        mass_diff = np.abs(region_mass_msun - region_mass_msun[anti_indices])
+        asymmetry = np.sum(mass_diff) / (2.0 * self.Mstar.to_value('Msun'))
+
+        return asymmetry
+
 
 class ApertureProperties(HaloProperty):
     """
@@ -3870,6 +3986,8 @@ class ApertureProperties(HaloProperty):
         "StellarInertiaTensorReducedLuminosityWeighted": True,
         "StellarInertiaTensorNoniterativeLuminosityWeighted": False,
         "StellarInertiaTensorReducedNoniterativeLuminosityWeighted": False,
+        "ShrinkingSphereCentre": False,
+        "StellarAsymmetry": False,
     }
 
     property_list = {
diff --git a/SOAP/property_table.py b/SOAP/property_table.py
index d86fdfe0..412e6ea3 100644
--- a/SOAP/property_table.py
+++ b/SOAP/property_table.py
@@ -3686,6 +3686,32 @@ class PropertyTable:
             output_physical=False,
             a_scale_exponent=1,
         ),
+        "ShrinkingSphereCentre": Property(
+            name="ShrinkingSphereCentre",
+            shape=3,
+            dtype=np.float64,
+            unit="snap_length",
+            # TODO: Add description
+            description="Centre of mass of stars.",
+            lossy_compression_filter="DScale6",
+            dmo_property=False,
+            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=False,
+            a_scale_exponent=1,
+        ),
+        "StellarAsymmetry": Property(
+            name="StellarAsymmetry",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            # TODO: Add description
+            description="Centre of mass of stars.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
         "compY": Property(
             name="ComptonY",
             shape=1,

From 5bea8159421afd3cee118d641530d2ce0851c370 Mon Sep 17 00:00:00 2001
From: robjmcgibbon <robjmcgibbon@gmail.com>
Date: Tue, 26 May 2026 09:39:54 +0100
Subject: [PATCH 2/7] Vary npix

---
 .../particle_selection/aperture_properties.py | 19 +++++++++++++-
 SOAP/property_table.py                        | 26 +++++++++++++++++++
 2 files changed, 44 insertions(+), 1 deletion(-)

diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index e1058db3..f4dad5dd 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -3587,6 +3587,7 @@ def stellar_inertia_tensor(self, **kwargs) -> unyt.unyt_array:
         # aperture
         mass = self.get_dataset("PartType4/Masses")[self.star_mask_all]
         position = (
+            # TODO: Remove shrinking sphere
             self.get_dataset("PartType4/Coordinates")[self.star_mask_all] - self.centre - self.shrinking_sphere_centre
         )
 
@@ -3725,6 +3726,8 @@ def shrinking_sphere_centre(
             Maximum number of iterations
         """
 
+        min_particles = min(200, 0.1*self.Nbaryon)
+
         if self.Mbaryons == 0:
             return None
 
@@ -3751,8 +3754,20 @@ def ShrinkingSphereCentre(self) -> unyt.unyt_array:
         """
         return (self.shrinking_sphere_centre + self.centre) % self.boxsize
 
+
+    @lazy_property
+    def StellarAsymmetry(self):
+        return self.stellar_asymmetry()
+
     @lazy_property
-    def StellarAsymmetry(self, npix=12):
+    def StellarAsymmetry48(self):
+        return self.stellar_asymmetry(npix=48)
+
+    @lazy_property
+    def StellarAsymmetry192(self):
+        return self.stellar_asymmetry(npix=192)
+
+    def stellar_asymmetry(self, npix=12):
         """
         Compute stellar asymmetry following https://arxiv.org/abs/1805.03210
 
@@ -3988,6 +4003,8 @@ class ApertureProperties(HaloProperty):
         "StellarInertiaTensorReducedNoniterativeLuminosityWeighted": False,
         "ShrinkingSphereCentre": False,
         "StellarAsymmetry": False,
+        "StellarAsymmetry48": False,
+        "StellarAsymmetry192": False,
     }
 
     property_list = {
diff --git a/SOAP/property_table.py b/SOAP/property_table.py
index 412e6ea3..6e489d03 100644
--- a/SOAP/property_table.py
+++ b/SOAP/property_table.py
@@ -3712,6 +3712,32 @@ class PropertyTable:
             output_physical=True,
             a_scale_exponent=0,
         ),
+        "StellarAsymmetry48": Property(
+            name="StellarAsymmetry48",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            # TODO: Add description
+            description="Centre of mass of stars.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
+        "StellarAsymmetry192": Property(
+            name="StellarAsymmetry192",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            # TODO: Add description
+            description="Centre of mass of stars.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
         "compY": Property(
             name="ComptonY",
             shape=1,

From e388cd238d7e44aece5616046c776057bf481e73 Mon Sep 17 00:00:00 2001
From: robjmcgibbon <robjmcgibbon@gmail.com>
Date: Wed, 3 Jun 2026 10:31:10 +0100
Subject: [PATCH 3/7] Squash bug found by Qinglin

---
 SOAP/particle_selection/aperture_properties.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index f4dad5dd..a005a359 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -3719,21 +3719,21 @@ def shrinking_sphere_centre(
         Parameters
         ----------
         min_particles : int, default=200
-            Stop iteration when enclosed particle count is below this threshold.
+            Stop iteration when enclosed particle count is below this threshold, or 10% of particle count.
         shrink_factor : float, default=0.83
             Radius scaling factor applied at each iteration (0 < shrink_factor < 1).
         max_iter : int, default=256
             Maximum number of iterations
         """
 
-        min_particles = min(200, 0.1*self.Nbaryon)
+        min_particles = min(min_particles, 0.1*self.Nbaryon)
 
         if self.Mbaryons == 0:
             return None
 
         centre = (self.baryon_mass_fraction[:, None] * self.pos_baryons).sum(axis=0)
         dist = np.linalg.norm(self.pos_baryons - centre, axis=1)
-        radius = max(self.aperture_radius, dist.max())
+        radius = np.max(dist) * 1.01
 
         for i in range(max_iter):
             mask = dist <= radius
@@ -3741,7 +3741,8 @@ def shrinking_sphere_centre(
             if n_in < min_particles:
                 break
 
-            centre = (self.baryon_mass_fraction[mask, None] * self.pos_baryons[mask]).sum(axis=0)
+            centre = (self.mass_baryons[mask, None] * self.pos_baryons[mask]).sum(axis=0)
+            centre /= np.sum(self.mass_baryons[mask])
             dist = np.linalg.norm(self.pos_baryons - centre, axis=1)
             radius *= shrink_factor
 
@@ -3752,6 +3753,9 @@ def ShrinkingSphereCentre(self) -> unyt.unyt_array:
         """
         Centre computed by applying shrinking spheres method to baryons
         """
+        if self.Mbaryons == 0:
+            return None
+
         return (self.shrinking_sphere_centre + self.centre) % self.boxsize
 
 

From 9849c0550c740c21f2997c60386b5466e5e23527 Mon Sep 17 00:00:00 2001
From: robjmcgibbon <robjmcgibbon@gmail.com>
Date: Mon, 15 Jun 2026 17:23:02 +0100
Subject: [PATCH 4/7] Add subsample

---
 .../particle_selection/aperture_properties.py | 45 ++++++++++++++++---
 SOAP/property_table.py                        | 26 +++++++++++
 2 files changed, 65 insertions(+), 6 deletions(-)

diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index a005a359..0bceb73d 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -3763,6 +3763,14 @@ def ShrinkingSphereCentre(self) -> unyt.unyt_array:
     def StellarAsymmetry(self):
         return self.stellar_asymmetry()
 
+    @lazy_property
+    def StellarAsymmetryShrink(self):
+        return self.stellar_asymmetry(shrink_centre=True)
+
+    @lazy_property
+    def StellarAsymmetrySubsample(self):
+        return self.stellar_asymmetry(N=8)
+
     @lazy_property
     def StellarAsymmetry48(self):
         return self.stellar_asymmetry(npix=48)
@@ -3771,7 +3779,7 @@ def StellarAsymmetry48(self):
     def StellarAsymmetry192(self):
         return self.stellar_asymmetry(npix=192)
 
-    def stellar_asymmetry(self, npix=12):
+    def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
         """
         Compute stellar asymmetry following https://arxiv.org/abs/1805.03210
 
@@ -3783,6 +3791,10 @@ def stellar_asymmetry(self, npix=12):
             Number of equal-area angular regions (HEALPix pixels) for directional
             partitioning. Must satisfy npix = 12 * nside^2 where nside is a
             power of 2 (e.g. 12, 48, 192, 768, ...).
+        N : int, optional
+            Randomly select 1/N of the original stellar particle set before
+            computing the asymmetry. If the total number of particles is smaller
+            than or equal to N, all particles are used.
 
         """
 
@@ -3794,8 +3806,26 @@ def stellar_asymmetry(self, npix=12):
         assert npix == 12 * nside ** 2
         assert nside.bit_count() == 1
 
+        if N is None:
+            indices = slice(None)
+        else:
+            N = int(N)
+            if N <= 0:
+                raise ValueError("N must be a positive integer")
+            if self.Nstar <= N:
+                indices = slice(None)
+            else:
+                indices = np.random.choice(
+                    self.Nstar, size=self.Nstar // N, replace=False
+                )
+
+        mass_star = self.mass_star[indices]
+
         # Centre using the shrinking sphere
-        pos = self.pos_star - self.shrinking_sphere_centre
+        if shrink_centre:
+            pos = self.pos_star[indices] - self.shrinking_sphere_centre
+        else:
+            pos = self.pos_star[indices]
         r = np.linalg.norm(pos, axis=1)
 
         # Remove particles close to the centre, they are symmetric
@@ -3803,11 +3833,12 @@ def stellar_asymmetry(self, npix=12):
         if np.sum(mask):
             pos = pos[np.logical_not(mask)]
             r = r[np.logical_not(mask)]
-            mass_star = self.mass_star[np.logical_not(mask)]
+            mass_star = mass_star[np.logical_not(mask)]
             if r.shape[0] == 0:
                 return np.float32(0)
-        else:
-            mass_star = self.mass_star
+
+        if mass_star.shape[0] < 3:
+            return None
 
         # Compute the mass in each pixel
         vecs = (pos / r[:, None]).value
@@ -3829,7 +3860,7 @@ def stellar_asymmetry(self, npix=12):
 
         # Calculate asymmetry
         mass_diff = np.abs(region_mass_msun - region_mass_msun[anti_indices])
-        asymmetry = np.sum(mass_diff) / (2.0 * self.Mstar.to_value('Msun'))
+        asymmetry = np.sum(mass_diff) / (2.0 * mass_star.sum().to_value('Msun'))
 
         return asymmetry
 
@@ -4007,6 +4038,8 @@ class ApertureProperties(HaloProperty):
         "StellarInertiaTensorReducedNoniterativeLuminosityWeighted": False,
         "ShrinkingSphereCentre": False,
         "StellarAsymmetry": False,
+        "StellarAsymmetryShrink": False,
+        "StellarAsymmetrySubsample": False,
         "StellarAsymmetry48": False,
         "StellarAsymmetry192": False,
     }
diff --git a/SOAP/property_table.py b/SOAP/property_table.py
index 6e489d03..0d766e00 100644
--- a/SOAP/property_table.py
+++ b/SOAP/property_table.py
@@ -3712,6 +3712,32 @@ class PropertyTable:
             output_physical=True,
             a_scale_exponent=0,
         ),
+        "StellarAsymmetryShrink": Property(
+            name="StellarAsymmetryShrink",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            # TODO: Add description
+            description="Centre of mass of stars.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
+        "StellarAsymmetrySubsample": Property(
+            name="StellarAsymmetrySubsample",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            # TODO: Add description
+            description="Centre of mass of stars.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
         "StellarAsymmetry48": Property(
             name="StellarAsymmetry48",
             shape=1,

From 9c17e33b9195956ca76ac16149b509be40d4a80b Mon Sep 17 00:00:00 2001
From: robjmcgibbon <robjmcgibbon@gmail.com>
Date: Fri, 26 Jun 2026 17:53:27 +0100
Subject: [PATCH 5/7] Shrink centre on stars

---
 .../particle_selection/aperture_properties.py | 21 +++++++++----------
 SOAP/property_table.py                        |  4 ++--
 2 files changed, 12 insertions(+), 13 deletions(-)

diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index 0bceb73d..42145188 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -3587,8 +3587,7 @@ def stellar_inertia_tensor(self, **kwargs) -> unyt.unyt_array:
         # aperture
         mass = self.get_dataset("PartType4/Masses")[self.star_mask_all]
         position = (
-            # TODO: Remove shrinking sphere
-            self.get_dataset("PartType4/Coordinates")[self.star_mask_all] - self.centre - self.shrinking_sphere_centre
+            self.get_dataset("PartType4/Coordinates")[self.star_mask_all] - self.centre
         )
 
         return get_inertia_tensor_mass_weighted(
@@ -3726,13 +3725,13 @@ def shrinking_sphere_centre(
             Maximum number of iterations
         """
 
-        min_particles = min(min_particles, 0.1*self.Nbaryon)
+        min_particles = min(min_particles, 0.1*self.Nstar)
 
-        if self.Mbaryons == 0:
+        if self.Mstar == 0:
             return None
 
-        centre = (self.baryon_mass_fraction[:, None] * self.pos_baryons).sum(axis=0)
-        dist = np.linalg.norm(self.pos_baryons - centre, axis=1)
+        centre = (self.star_mass_fraction[:, None] * self.pos_star).sum(axis=0)
+        dist = np.linalg.norm(self.pos_star - centre, axis=1)
         radius = np.max(dist) * 1.01
 
         for i in range(max_iter):
@@ -3741,9 +3740,9 @@ def shrinking_sphere_centre(
             if n_in < min_particles:
                 break
 
-            centre = (self.mass_baryons[mask, None] * self.pos_baryons[mask]).sum(axis=0)
-            centre /= np.sum(self.mass_baryons[mask])
-            dist = np.linalg.norm(self.pos_baryons - centre, axis=1)
+            centre = (self.mass_star[mask, None] * self.pos_star[mask]).sum(axis=0)
+            centre /= np.sum(self.mass_star[mask])
+            dist = np.linalg.norm(self.pos_star - centre, axis=1)
             radius *= shrink_factor
 
         return centre
@@ -3751,9 +3750,9 @@ def shrinking_sphere_centre(
     @lazy_property
     def ShrinkingSphereCentre(self) -> unyt.unyt_array:
         """
-        Centre computed by applying shrinking spheres method to baryons
+        Centre computed by applying shrinking spheres method to stars
         """
-        if self.Mbaryons == 0:
+        if self.Mstar == 0:
             return None
 
         return (self.shrinking_sphere_centre + self.centre) % self.boxsize
diff --git a/SOAP/property_table.py b/SOAP/property_table.py
index 0d766e00..e7e8c9fc 100644
--- a/SOAP/property_table.py
+++ b/SOAP/property_table.py
@@ -3692,10 +3692,10 @@ class PropertyTable:
             dtype=np.float64,
             unit="snap_length",
             # TODO: Add description
-            description="Centre of mass of stars.",
+            description="Shrinking sphere centre computed using stars.",
             lossy_compression_filter="DScale6",
             dmo_property=False,
-            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
             output_physical=False,
             a_scale_exponent=1,
         ),

From 8b18a4d60247611950a25914f8fa4927c7d848e3 Mon Sep 17 00:00:00 2001
From: robjmcgibbon <robjmcgibbon@gmail.com>
Date: Fri, 25 Sep 2026 12:57:11 +0100
Subject: [PATCH 6/7] Update to match changes to master

---
 README.md                                     | 32 +++++++++++++--
 .../particle_selection/aperture_properties.py | 23 +++++------
 SOAP/property_table.py                        | 39 +++++++++++--------
 documentation/footnote_asymmetry.tex          | 18 +++++++++
 pyproject.toml                                |  5 +++
 5 files changed, 83 insertions(+), 34 deletions(-)
 create mode 100644 documentation/footnote_asymmetry.tex

diff --git a/README.md b/README.md
index dbc66807..494f37c7 100644
--- a/README.md
+++ b/README.md
@@ -18,11 +18,20 @@ The code is written in python and uses mpi4py for parallelism. Whichever
 install method you use, you will need MPI available so that mpi4py can be
 installed.
 
+We recommend cloning the repository, rather than installing SOAP directly from
+GitHub with pip. The example parameter files, test scripts, and the scripts used
+to generate the documentation are only available in the repository, and the
+commands in [Running SOAP](#running-soap) are run from the root of the
+repository. If you only want to import SOAP as a library, then it can be
+installed with `pip install git+https://github.com/SWIFTSIM/SOAP.git`.
+
 ### Quick install (serial HDF5)
 
-SOAP and its dependencies can be installed directly using the command
+SOAP and its dependencies can be installed using the commands
 ```
-pip install git+https://github.com/SWIFTSIM/SOAP.git
+git clone https://github.com/SWIFTSIM/SOAP.git
+cd SOAP
+pip install .
 ```
 This will usually install a serial version of h5py. SOAP will run correctly,
 but for large simulations the I/O will be slower.
@@ -36,11 +45,22 @@ source against that library before installing SOAP:
 ```
 pip install mpi4py
 export HDF5_MPI="ON"; export CC=mpicc; pip install --no-binary=h5py h5py
-pip install git+https://github.com/SWIFTSIM/SOAP.git
+git clone https://github.com/SWIFTSIM/SOAP.git
+cd SOAP
+pip install .
 ```
 If SOAP (and therefore serial h5py) is already installed, then add the flags
 `--no-cache-dir` and `--force-reinstall` when reinstalling h5py.
 
+### Optional dependencies
+
+Some properties require additional python packages. These properties are
+not computed unless they are explicitly enabled in the parameter file. To
+install the packages needed for all of these properties, run
+```
+pip install ".[extra_properties]"
+```
+
 ### Installation on COSMA
 
 If you are using the [COSMA system](https://cosma.readthedocs.io/en/latest/),
@@ -49,6 +69,9 @@ you can install an SOAP virtual environment by running
 
 ## Running SOAP
 
+The commands in this section should be run from the root of the SOAP
+repository.
+
 The command `./tests/run_small_volume.sh` will download a small example
 simulation, run the group membership and halo properties scripts on it.
 This uses the parameter file at `./tests/small_volume.yml`, and the
@@ -238,7 +261,8 @@ If your property is expensive to compute, or requires additional dependencies,
 set `opt_in_reason` in its `SOAP/property_table.py` entry to a short (20
 characters or less) reason. It will then only be calculated if it is explicitly
 enabled in the parameter file. Any additional dependencies must be imported
-within the `@lazy_property`, not at the top of the file.
+within the `@lazy_property`, not at the top of the file, and added to the
+`extra_properties` section of `pyproject.toml`.
 
 If SOAP does crash while evaluating your new property it will try to
 output the ID of the halo it was processing when it crashed. Then you
diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index 8cd35b49..7972068f 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -137,7 +137,6 @@ def MetalFracStar(self):
 import numpy as np
 from numpy.typing import NDArray
 from typing import Dict, List, Tuple
-import healpy as hp
 import unyt
 
 from .halo_properties import HaloProperty, SearchRadiusTooSmallError
@@ -3652,11 +3651,8 @@ def StellarInertiaTensorReducedNoniterativeLuminosityWeighted(
 
     @lazy_property
     def shrinking_sphere_centre(
-            self,
-            min_particles=200,
-            shrink_factor=0.83,
-            max_iter = 256
-        ) -> unyt.unyt_array:
+        self, min_particles=200, shrink_factor=0.83, max_iter=256
+    ) -> unyt.unyt_array:
         """
         Estimate the galaxy center (center of mass) with iterative shrinking aperture.
 
@@ -3676,7 +3672,7 @@ def shrinking_sphere_centre(
             Maximum number of iterations
         """
 
-        min_particles = min(min_particles, 0.1*self.Nstar)
+        min_particles = min(min_particles, 0.1 * self.Nstar)
 
         if self.Mstar == 0:
             return None
@@ -3708,7 +3704,6 @@ def ShrinkingSphereCentre(self) -> unyt.unyt_array:
 
         return (self.shrinking_sphere_centre + self.centre) % self.boxsize
 
-
     @lazy_property
     def StellarAsymmetry(self):
         return self.stellar_asymmetry()
@@ -3733,8 +3728,6 @@ def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
         """
         Compute stellar asymmetry following https://arxiv.org/abs/1805.03210
 
-        TODO: Is equation (3) incorrect?
-
         Parameters
         ----------
         npix : int, default=12
@@ -3748,12 +3741,14 @@ def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
 
         """
 
+        import healpy as hp
+
         if self.Mstar == 0:
             return None
 
         # Check we are using a valid value for npix
         nside = int(round(np.sqrt(npix // 12)))
-        assert npix == 12 * nside ** 2
+        assert npix == 12 * nside**2
         assert nside.bit_count() == 1
 
         if N is None:
@@ -3779,7 +3774,7 @@ def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
         r = np.linalg.norm(pos, axis=1)
 
         # Remove particles close to the centre, they are symmetric
-        mask = r.to_value('kpc') < 0.1
+        mask = r.to_value("kpc") < 0.1
         if np.sum(mask):
             pos = pos[np.logical_not(mask)]
             r = r[np.logical_not(mask)]
@@ -3796,7 +3791,7 @@ def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
         # np.bincount will not return a unyt array
         region_mass_msun = np.bincount(
             idx,
-            weights=mass_star.to_value('Msun'),
+            weights=mass_star.to_value("Msun"),
             minlength=npix,
         )
 
@@ -3810,7 +3805,7 @@ def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
 
         # Calculate asymmetry
         mass_diff = np.abs(region_mass_msun - region_mass_msun[anti_indices])
-        asymmetry = np.sum(mass_diff) / (2.0 * mass_star.sum().to_value('Msun'))
+        asymmetry = np.sum(mass_diff) / (2.0 * mass_star.sum().to_value("Msun"))
 
         return asymmetry
 
diff --git a/SOAP/property_table.py b/SOAP/property_table.py
index e714470f..ad9aad7b 100644
--- a/SOAP/property_table.py
+++ b/SOAP/property_table.py
@@ -222,6 +222,13 @@ class PropertyTable:
         ],
         "footnote_compY.tex": ["compY", "compY_no_agn"],
         "footnote_dopplerB.tex": ["DopplerB"],
+        "footnote_asymmetry.tex": [
+            "StellarAsymmetry",
+            "StellarAsymmetryShrink",
+            "StellarAsymmetrySubsample",
+            "StellarAsymmetry48",
+            "StellarAsymmetry192",
+        ],
         "footnote_coreexcision.tex": [
             "Tgas_cy_weighted_core_excision",
             "Tgas_cy_weighted_core_excision_no_agn",
@@ -3759,78 +3766,78 @@ class PropertyTable:
             shape=3,
             dtype=np.float64,
             unit="snap_length",
-            # TODO: Add description
             description="Shrinking sphere centre computed using stars.",
             lossy_compression_filter="DScale6",
             dmo_property=False,
             particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
             output_physical=False,
             a_scale_exponent=1,
+            opt_in_reason="Expensive to compute",
         ),
         "StellarAsymmetry": Property(
             name="StellarAsymmetry",
             shape=1,
             dtype=np.float32,
             unit="dimensionless",
-            # TODO: Add description
-            description="Centre of mass of stars.",
+            description="Asymmetry of the stellar mass distribution around the halo centre, computed using 12 HEALPix pixels.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
-            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
             output_physical=True,
             a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
         ),
         "StellarAsymmetryShrink": Property(
             name="StellarAsymmetryShrink",
             shape=1,
             dtype=np.float32,
             unit="dimensionless",
-            # TODO: Add description
-            description="Centre of mass of stars.",
+            description="As StellarAsymmetry, but centred on ShrinkingSphereCentre rather than the halo centre.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
-            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
             output_physical=True,
             a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
         ),
         "StellarAsymmetrySubsample": Property(
             name="StellarAsymmetrySubsample",
             shape=1,
             dtype=np.float32,
             unit="dimensionless",
-            # TODO: Add description
-            description="Centre of mass of stars.",
+            description="As StellarAsymmetry, but computed using a random subsample of 1/8 of the star particles. All star particles are used if there are 8 or fewer.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
-            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
             output_physical=True,
             a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
         ),
         "StellarAsymmetry48": Property(
             name="StellarAsymmetry48",
             shape=1,
             dtype=np.float32,
             unit="dimensionless",
-            # TODO: Add description
-            description="Centre of mass of stars.",
+            description="As StellarAsymmetry, but using 48 HEALPix pixels.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
-            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
             output_physical=True,
             a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
         ),
         "StellarAsymmetry192": Property(
             name="StellarAsymmetry192",
             shape=1,
             dtype=np.float32,
             unit="dimensionless",
-            # TODO: Add description
-            description="Centre of mass of stars.",
+            description="As StellarAsymmetry, but using 192 HEALPix pixels.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
-            particle_properties=["PartType0/Coordinates", "PartType0/Masses", "PartType4/Coordinates", "PartType4/Masses"],
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
             output_physical=True,
             a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
         ),
         "compY": Property(
             name="ComptonY",
diff --git a/documentation/footnote_asymmetry.tex b/documentation/footnote_asymmetry.tex
new file mode 100644
index 00000000..bec8af86
--- /dev/null
+++ b/documentation/footnote_asymmetry.tex
@@ -0,0 +1,18 @@
+\paragraph{$^{$FOOTNOTE_NUMBER$}$The stellar asymmetry}\label{footnote:$FOOTNOTE_NUMBER$} is computed following
+\href{https://arxiv.org/abs/1805.03210}{arXiv:1805.03210}. The directions around the centre are divided into
+$N_{\rm{}pix}$ equal-area HEALPix pixels, and the asymmetry is
+
+\begin{equation}
+    A = \frac{1}{2 M_*} \sum_{i=1}^{N_{\rm{}pix}} \left| M_i - M_{\bar{i}} \right|,
+\end{equation}
+
+where $M_i$ is the stellar mass in pixel $i$, $M_{\bar{i}}$ is the stellar mass in the pixel in the opposite
+direction, and $M_*$ is the total stellar mass. Each pair of opposite pixels appears twice in the sum, so $A$
+ranges from 0 for a perfectly symmetric distribution to 1 if all the stellar mass is on one side of the centre.
+Star particles within 0.1 kpc of the centre are excluded, since their direction from the centre is poorly
+defined. The asymmetry is not calculated if there are fewer than 3 star particles.
+
+\verb+StellarAsymmetry+ uses $N_{\rm{}pix}=12$ and is centred on the halo centre. The other variants differ
+as follows: \verb+StellarAsymmetryShrink+ is centred on \verb+ShrinkingSphereCentre+,
+\verb+StellarAsymmetrySubsample+ uses a random subsample of 1/8 of the star particles, and
+\verb+StellarAsymmetry48+ and \verb+StellarAsymmetry192+ use $N_{\rm{}pix}=48$ and $N_{\rm{}pix}=192$.
diff --git a/pyproject.toml b/pyproject.toml
index 44f13618..3847a191 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -32,9 +32,14 @@ dependencies = [
 ]
 
 [project.optional-dependencies]
+# Dependencies which are only required by some opt-in properties
+extra_properties = [
+    "healpy",
+]
 test = [
     "pytest-mpi",
     "numba",
+    "SOAP[extra_properties]",
 ]
 
 [project.urls]

From 3dca05858f18e0ea6ef9f7d834dcfd1827e0746f Mon Sep 17 00:00:00 2001
From: robjmcgibbon <robjmcgibbon@gmail.com>
Date: Mon, 28 Sep 2026 10:23:15 +0100
Subject: [PATCH 7/7] Use fix seed for subsampling

---
 SOAP/particle_selection/aperture_properties.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index 7972068f..54029599 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -3760,9 +3760,9 @@ def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
             if self.Nstar <= N:
                 indices = slice(None)
             else:
-                indices = np.random.choice(
-                    self.Nstar, size=self.Nstar // N, replace=False
-                )
+                # Seed with the halo index so results are reproducible
+                rng = np.random.default_rng(int(self.index))
+                indices = rng.choice(self.Nstar, size=self.Nstar // N, replace=False)
 
         mass_star = self.mass_star[indices]