Compute stellar asymmetry

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 = {

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,