From d72181aaa7a6114f6c5e0a2fbc12cb6374de012e Mon Sep 17 00:00:00 2001 From: robjmcgibbon Date: Wed, 22 Apr 2026 17:00:27 +0100 Subject: [PATCH] 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,