diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py
index 6a07f6cb..5d0cbe4b 100644
--- a/SOAP/particle_selection/aperture_properties.py
+++ b/SOAP/particle_selection/aperture_properties.py
@@ -147,6 +147,10 @@ def MetalFracStar(self):
 )
 from SOAP.property_calculation.kinematic_properties import (
     get_velocity_dispersion_matrix,
+    get_cylindrical_velocity_dispersion_vector_mass_weighted,
+    get_cylindrical_velocity_dispersion_vector_luminosity_weighted,
+    get_rotation_velocity_mass_weighted,
+    get_rotation_velocity_luminosity_weighted,
     get_angular_momentum,
     get_angular_momentum_and_kappa_corot_mass_weighted,
     get_angular_momentum_and_kappa_corot_luminosity_weighted,
@@ -156,6 +160,9 @@ def MetalFracStar(self):
     get_inertia_tensor_mass_weighted,
     get_inertia_tensor_luminosity_weighted,
 )
+from SOAP.property_calculation.cylindrical_coordinates import (
+    calculate_cylindrical_velocities,
+)
 from SOAP.core.swift_cells import SWIFTCellGrid
 from SOAP.property_calculation.stellar_age_calculator import StellarAgeCalculator
 from SOAP.particle_filter.cold_dense_gas_filter import ColdDenseGasFilter
@@ -1448,6 +1455,162 @@ def veldisp_matrix_star(self) -> unyt.unyt_array:
             self.star_mass_fraction, self.vel_star, self.vcom_star
         )
 
+    @lazy_property
+    def star_cylindrical_velocities(self) -> unyt.unyt_array:
+        """
+        Calculate the velocities of the star particles in cyclindrical
+        coordinates, where the axes are centred on the stellar CoM,
+        and the z axis is aligned with the stellar angular momentum.
+        """
+
+        # We need at least 2 particles to have an angular momentum vector
+        if self.Nstar < 2:
+            return None
+
+        # This can happen if we have particles on top of each other
+        # or with the same velocity
+        if np.sum(self.Lstar) == 0:
+            return None
+
+        # Get velocities in cylindrical coordinates
+        return calculate_cylindrical_velocities(
+            self.pos_star,
+            self.vel_star,
+            self.Lstar,
+            reference_velocity=self.vcom_star,
+        )
+
+    @lazy_property
+    def StellarRotationalVelocity(self) -> unyt.unyt_array:
+        """
+        Mass-weighted average azimuthal velocity of stars.
+        """
+        if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
+            return None
+        return get_rotation_velocity_mass_weighted(
+            self.mass_star, self.star_cylindrical_velocities[:, 1]
+        )
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionVector(self) -> unyt.unyt_array:
+        if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
+            return None
+        return get_cylindrical_velocity_dispersion_vector_mass_weighted(
+            self.mass_star, self.star_cylindrical_velocities
+        )
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersion(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVector is None:
+            return None
+        return np.sqrt((self.StellarCylindricalVelocityDispersionVector**2).sum() / 3)
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionVertical(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVector is None:
+            return None
+        return self.StellarCylindricalVelocityDispersionVector[2]
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionDiscPlane(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVector is None:
+            return None
+        return np.sqrt((self.StellarCylindricalVelocityDispersionVector[:2] ** 2).sum())
+
+    @lazy_property
+    def star_cylindrical_velocities_luminosity_weighted(self) -> unyt.unyt_array:
+        """
+        Calculate the velocities of the star particles in cylindrical
+        coordinates, where the origin and reference frame are centred on the
+        stellar centre of light in each band. The z axis is aligned with the
+        stellar angular momentum obtained from each band.
+        """
+
+        # We need at least 2 particles to have an angular momentum vector
+        if self.Nstar < 2:
+            return None
+
+        # This can happen if we have particles on top of each other
+        # or with the same velocity
+        if np.sum(self.Lstar_luminosity_weighted) == 0:
+            return None
+
+        # We iterate over bands to use their own reference vector and luminosity-
+        # weighted centre of mass phase space coordinates.
+        cylindrical_velocities = (
+            np.zeros(
+                (
+                    self.stellar_luminosities.shape[1],
+                    self.stellar_luminosities.shape[0],
+                    3,
+                )
+            )
+            * self.vel_star.units
+        )
+        for i_band, particle_luminosities_i_band in enumerate(
+            self.stellar_luminosities.T
+        ):
+            cylindrical_velocities[i_band] = calculate_cylindrical_velocities(
+                self.pos_star,
+                self.vel_star,
+                self.Lstar_luminosity_weighted[i_band * 3 : (1 + i_band) * 3],
+                reference_velocity=self.vcom_star,
+            )
+
+        return cylindrical_velocities
+
+    @lazy_property
+    def StellarRotationalVelocityLuminosityWeighted(self) -> unyt.unyt_array:
+        if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
+            return None
+        return get_rotation_velocity_luminosity_weighted(
+            self.stellar_luminosities,
+            self.star_cylindrical_velocities_luminosity_weighted[:, :, 1],
+        )
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionVectorLuminosityWeighted(
+        self,
+    ) -> unyt.unyt_array:
+        if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
+            return None
+        return get_cylindrical_velocity_dispersion_vector_luminosity_weighted(
+            self.stellar_luminosities,
+            self.star_cylindrical_velocities_luminosity_weighted,
+        )
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionLuminosityWeighted(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None:
+            return None
+        return np.sqrt(
+            (
+                self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted**2
+            ).sum(axis=1)
+            / 3
+        )
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionVerticalLuminosityWeighted(
+        self,
+    ) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None:
+            return None
+        return self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, 2]
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted(
+        self,
+    ) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None:
+            return None
+        return np.sqrt(
+            (
+                self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, :2]
+                ** 2
+            ).sum(axis=1)
+        )
+
     @lazy_property
     def KineticEnergyStars(self) -> unyt.unyt_quantity:
         """
@@ -3583,6 +3746,14 @@ class ApertureProperties(HaloProperty):
         "veldisp_matrix_gas": False,
         "veldisp_matrix_dm": False,
         "veldisp_matrix_star": False,
+        "StellarRotationalVelocity": False,
+        "StellarCylindricalVelocityDispersion": False,
+        "StellarCylindricalVelocityDispersionVertical": False,
+        "StellarCylindricalVelocityDispersionDiscPlane": False,
+        "StellarRotationalVelocityLuminosityWeighted": False,
+        "StellarCylindricalVelocityDispersionLuminosityWeighted": False,
+        "StellarCylindricalVelocityDispersionVerticalLuminosityWeighted": False,
+        "StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted": False,
         "KineticEnergyGas": False,
         "KineticEnergyStars": False,
         "Mgas_SF": False,
diff --git a/SOAP/particle_selection/subhalo_properties.py b/SOAP/particle_selection/subhalo_properties.py
index a11231f5..bd2c94cd 100644
--- a/SOAP/particle_selection/subhalo_properties.py
+++ b/SOAP/particle_selection/subhalo_properties.py
@@ -33,6 +33,10 @@
     get_angular_momentum_and_kappa_corot_luminosity_weighted,
     get_vmax,
     get_velocity_dispersion_matrix,
+    get_cylindrical_velocity_dispersion_vector_mass_weighted,
+    get_cylindrical_velocity_dispersion_vector_luminosity_weighted,
+    get_rotation_velocity_mass_weighted,
+    get_rotation_velocity_luminosity_weighted,
 )
 from SOAP.property_calculation.inertia_tensors import (
     get_inertia_tensor_mass_weighted,
@@ -1389,7 +1393,7 @@ def star_mass_fraction(self) -> unyt.unyt_array:
     @lazy_property
     def star_cylindrical_velocities(self) -> unyt.unyt_array:
         """
-        Calculate the velocities of the star particles in cyclindrical
+        Calculate the velocities of the star particles in cylindrical
         coordinates, where the axes are centred on the stellar CoM,
         and the z axis is aligned with the stellar angular momentum.
         """
@@ -1403,63 +1407,141 @@ def star_cylindrical_velocities(self) -> unyt.unyt_array:
         if np.sum(self.Lstar) == 0:
             return None
 
-        # Calculate the position of the stars relative to their CoM
-        pos = self.pos_star - (self.star_mass_fraction[:, None] * self.pos_star).sum(
-            axis=0
-        )
-
-        # Calculate relative velocity of stars
-        vrel = self.vel_star - self.vcom[None, :]
-
         # Get velocities in cylindrical coordinates
         return calculate_cylindrical_velocities(
-            pos,
-            vrel,
+            self.pos_star,
+            self.vel_star,
             self.Lstar,
+            reference_velocity=self.vcom_star,
         )
 
     @lazy_property
     def StellarRotationalVelocity(self) -> unyt.unyt_array:
         if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
             return None
-        v_cylindrical = self.star_cylindrical_velocities
-        v_phi = v_cylindrical[:, 1]
-        return (self.star_mass_fraction * v_phi).sum()
+        return get_rotation_velocity_mass_weighted(
+            self.mass_star, self.star_cylindrical_velocities[:, 1]
+        )
 
     @lazy_property
-    def stellar_cylindrical_squared_velocity_dispersion_vector(self) -> unyt.unyt_array:
+    def StellarCylindricalVelocityDispersionVector(self) -> unyt.unyt_array:
         if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
             return None
-        v_cylindrical = self.star_cylindrical_velocities
+        return get_cylindrical_velocity_dispersion_vector_mass_weighted(
+            self.mass_star, self.star_cylindrical_velocities
+        )
 
-        # This implementation of standard deviation is more numerically stable than using <x^2> - <x>^2
-        mean_velocity = (self.star_mass_fraction[:, None] * v_cylindrical).sum(axis=0)
-        squared_velocity_dispersion = (
-            self.star_mass_fraction[:, None] * (v_cylindrical - mean_velocity) ** 2
-        ).sum(axis=0)
+    @lazy_property
+    def StellarCylindricalVelocityDispersion(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVector is None:
+            return None
+        return np.sqrt((self.StellarCylindricalVelocityDispersionVector**2).sum() / 3)
 
-        return squared_velocity_dispersion
+    @lazy_property
+    def StellarCylindricalVelocityDispersionVertical(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVector is None:
+            return None
+        return self.StellarCylindricalVelocityDispersionVector[2]
 
     @lazy_property
-    def StellarCylindricalVelocityDispersion(self) -> unyt.unyt_array:
-        if self.stellar_cylindrical_squared_velocity_dispersion_vector is None:
+    def StellarCylindricalVelocityDispersionDiscPlane(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVector is None:
+            return None
+        return np.sqrt((self.StellarCylindricalVelocityDispersionVector[:2] ** 2).sum())
+
+    @lazy_property
+    def star_cylindrical_velocities_luminosity_weighted(self) -> unyt.unyt_array:
+        """
+        Calculate the velocities of the star particles in cylindrical
+        coordinates, where the origin and reference frame are centred on the
+        stellar centre of light in each band. The z axis is aligned with the
+        stellar angular momentum obtained from each band.
+        """
+
+        # We need at least 2 particles to have an angular momentum vector
+        if self.Nstar < 2:
+            return None
+
+        # This can happen if we have particles on top of each other
+        # or with the same velocity
+        if np.sum(self.Lstar_luminosity_weighted) == 0:
+            return None
+
+        # We iterate over bands to use their own reference vector and luminosity-
+        # weighted centre of mass phase space coordinates.
+        cylindrical_velocities = (
+            np.zeros(
+                (
+                    self.stellar_luminosities.shape[1],
+                    self.stellar_luminosities.shape[0],
+                    3,
+                )
+            )
+            * self.vel_star.units
+        )
+        for i_band, particle_luminosities_i_band in enumerate(
+            self.stellar_luminosities.T
+        ):
+            cylindrical_velocities[i_band] = calculate_cylindrical_velocities(
+                self.pos_star,
+                self.vel_star,
+                self.Lstar_luminosity_weighted[i_band * 3 : (1 + i_band) * 3],
+                reference_velocity=self.vcom_star,
+            )
+
+        return cylindrical_velocities
+
+    @lazy_property
+    def StellarRotationalVelocityLuminosityWeighted(self) -> unyt.unyt_array:
+        if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
+            return None
+
+        return get_rotation_velocity_luminosity_weighted(
+            self.stellar_luminosities,
+            self.star_cylindrical_velocities_luminosity_weighted[:, :, 1],
+        )
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionVectorLuminosityWeighted(
+        self,
+    ) -> unyt.unyt_array:
+        if (self.Nstar < 2) or (np.sum(self.Lstar) == 0):
+            return None
+        return get_cylindrical_velocity_dispersion_vector_luminosity_weighted(
+            self.stellar_luminosities,
+            self.star_cylindrical_velocities_luminosity_weighted,
+        )
+
+    @lazy_property
+    def StellarCylindricalVelocityDispersionLuminosityWeighted(self) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None:
             return None
         return np.sqrt(
-            self.stellar_cylindrical_squared_velocity_dispersion_vector.sum() / 3
+            (
+                self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted**2
+            ).sum(axis=1)
+            / 3
         )
 
     @lazy_property
-    def StellarCylindricalVelocityDispersionVertical(self) -> unyt.unyt_array:
-        if self.stellar_cylindrical_squared_velocity_dispersion_vector is None:
+    def StellarCylindricalVelocityDispersionVerticalLuminosityWeighted(
+        self,
+    ) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None:
             return None
-        return np.sqrt(self.stellar_cylindrical_squared_velocity_dispersion_vector[2])
+        return self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, 2]
 
     @lazy_property
-    def StellarCylindricalVelocityDispersionDiscPlane(self) -> unyt.unyt_array:
-        if self.stellar_cylindrical_squared_velocity_dispersion_vector is None:
+    def StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted(
+        self,
+    ) -> unyt.unyt_array:
+        if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None:
             return None
         return np.sqrt(
-            self.stellar_cylindrical_squared_velocity_dispersion_vector[:2].sum()
+            (
+                self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, :2]
+                ** 2
+            ).sum(axis=1)
         )
 
     @lazy_property
@@ -2307,10 +2389,14 @@ class SubhaloProperties(HaloProperty):
             "Lgas",
             "Ldm",
             "Lstar",
+            "StellarRotationalVelocity",
             "StellarCylindricalVelocityDispersion",
             "StellarCylindricalVelocityDispersionVertical",
             "StellarCylindricalVelocityDispersionDiscPlane",
-            "StellarRotationalVelocity",
+            "StellarRotationalVelocityLuminosityWeighted",
+            "StellarCylindricalVelocityDispersionLuminosityWeighted",
+            "StellarCylindricalVelocityDispersionVerticalLuminosityWeighted",
+            "StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted",
             "kappa_corot_gas",
             "kappa_corot_star",
             "Lbaryons",
diff --git a/SOAP/property_calculation/cylindrical_coordinates.py b/SOAP/property_calculation/cylindrical_coordinates.py
index 223ac266..410b01c4 100644
--- a/SOAP/property_calculation/cylindrical_coordinates.py
+++ b/SOAP/property_calculation/cylindrical_coordinates.py
@@ -42,7 +42,9 @@ def build_rotation_matrix(z_target):
     return R
 
 
-def calculate_cylindrical_velocities(positions, velocities, z_target):
+def calculate_cylindrical_velocities(
+    positions, velocities, z_target, reference_position=None, reference_velocity=None
+):
     """
     Convert 3D Cartesian velocities to cylindrical coordinates (v_r, v_phi, v_z),
     after rotating the system such that the z-axis aligns with `z_target`.
@@ -51,16 +53,30 @@ def calculate_cylindrical_velocities(positions, velocities, z_target):
         positions: (N, 3) array of particle positions in the original Cartesian frame.
         velocities: (N, 3) array of particle velocities in the original Cartesian frame.
         z_target: A 3-element vector indicating the new z-axis direction.
+        reference_position: (3,) array with a reference position on which to centre the Cartesian coordinate system.
+        reference_velocity: (3,) array with a reference velocity on which to centre the Cartesian reference frame.
 
     Returns:
         cyl_velocities: (N, 3) array of velocities in cylindrical coordinates:
               [v_r, v_phi, v_z] for each particle.
     """
+
+    # We need to declare a relative pos/vel array to not overwrite the original
+    # values.
+    if reference_position is None:
+        prel = positions
+    else:
+        prel = positions - reference_position
+    if reference_velocity is None:
+        vrel = velocities
+    else:
+        vrel = velocities - reference_velocity
+
     R = build_rotation_matrix(z_target)
 
     # Rotate positions and velocities into new frame
-    positions_rot = positions @ R.T
-    velocities_rot = velocities @ R.T
+    positions_rot = prel @ R.T
+    velocities_rot = vrel @ R.T
 
     x = positions_rot[:, 0]
     y = positions_rot[:, 1]
diff --git a/SOAP/property_calculation/half_mass_radius.py b/SOAP/property_calculation/half_mass_radius.py
index 1daecf1d..13896b16 100644
--- a/SOAP/property_calculation/half_mass_radius.py
+++ b/SOAP/property_calculation/half_mass_radius.py
@@ -52,7 +52,7 @@ def get_half_weight_radius(
     # Consistency check.
     # np.sum() and np.cumsum() use different orders, so we have to allow for
     # some small difference.
-    if cumulative_weights[-1] < 0.999 * total_weight:
+    if cumulative_weights[-1] < 0.998 * total_weight:
         raise RuntimeError(
             "Weights sum up to less than the given total weight value:"
             f" cumulative_weights[-1] = {cumulative_weights[-1]},"
diff --git a/SOAP/property_calculation/kinematic_properties.py b/SOAP/property_calculation/kinematic_properties.py
index 14aeab68..1e5c900f 100644
--- a/SOAP/property_calculation/kinematic_properties.py
+++ b/SOAP/property_calculation/kinematic_properties.py
@@ -14,6 +14,80 @@
 import unyt
 
 
+def get_weighted_rotation_velocity(
+    particle_weights: unyt.unyt_array, particle_azimuthal_velocities: unyt.unyt_array
+) -> unyt.unyt_quantity:
+    """
+    Get the weighted average azimuthal velocity of a particle distribution.
+
+    Parameters:
+     - particle_weights: unyt.unyt_array
+       Weight assigned to each particle.
+     - particle_azimuthal_velocities: unyt.unyt_array
+       Azimuthal velocity of each particle.
+
+    Returns:
+     - Weighted average of the azimuthal velocity of particles.
+    """
+    return (particle_weights * particle_azimuthal_velocities).sum()
+
+
+def get_rotation_velocity_mass_weighted(
+    particle_masses, particle_azimuthal_velocities
+) -> unyt.unyt_quantity:
+    """
+    Return the mass-weighted average azimuthal velocity of a particle distribution.
+
+    Parameters:
+     - particle_masses: unyt.unyt_array
+       Mass of particle.
+     - particle_azimuthal_velocities: unyt.unyt_array
+       Azimuthal velocity of each particle.
+
+    Returns:
+     - Mass-weighted average of the azimuthal velocity of particles.
+    """
+    mass_weights = particle_masses / particle_masses.sum()
+    return get_weighted_rotation_velocity(mass_weights, particle_azimuthal_velocities)
+
+
+def get_rotation_velocity_luminosity_weighted(
+    particle_luminosities: unyt.unyt_array,
+    particle_azimuthal_velocities: unyt.unyt_array,
+) -> unyt.unyt_array:
+    """
+    Return the luminosity-weighted average azimuthal velocity of a particle distribution, for each
+    provided luminosity band.
+
+    Parameters:
+     - particle_luminosities: unyt.unyt_array
+       Luminosity of each particle in the provided bands.
+     - particle_azimuthal_velocities: unyt.unyt_array
+       Azimuthal velocity of each particle, pre-computed for each luminosity band.
+
+    Returns:
+     - Luminosity-weighted average of the azimuthal velocity of particles, with a value for each band.
+    """
+
+    number_luminosity_bands = particle_luminosities.shape[1]
+    rotation_velocities = (
+        np.zeros(number_luminosity_bands) * particle_azimuthal_velocities.units
+    )
+
+    for i_band, (
+        particle_luminosities_i_band,
+        particle_azimuthal_velocities_i_band,
+    ) in enumerate(zip(particle_luminosities.T, particle_azimuthal_velocities)):
+        luminosity_weights = (
+            particle_luminosities_i_band / particle_luminosities_i_band.sum()
+        )
+        rotation_velocities[i_band] = get_weighted_rotation_velocity(
+            luminosity_weights, particle_azimuthal_velocities_i_band
+        )
+
+    return rotation_velocities
+
+
 def get_velocity_dispersion_matrix(
     mass_fraction: unyt.unyt_array,
     velocity: unyt.unyt_array,
@@ -53,6 +127,98 @@ def get_velocity_dispersion_matrix(
     return result
 
 
+def get_weighted_cylindrical_velocity_dispersion_vector(
+    particle_weights: unyt.unyt_array,
+    particle_cylindrical_velocities: unyt.unyt_array,
+) -> unyt.unyt_array:
+    """
+    Compute the velocity dispersion along the radial, azimuthal and vertical
+    directions for the input particles and their specified weights.
+
+    Parameters:
+     - particle_weights: unyt.unyt_array
+       Weight assigned to each particle.
+     - particle_cylindrical_velocities: unyt.unyt_array
+       Velocity of the particles in a cylindrical coordinate system.
+
+    Returns a 3 element vector containing [sigma_r, sigma_phi, sigma_z].
+    """
+
+    # This implementation of standard deviation is more numerically stable than using <x^2> - <x>^2
+    mean_velocity = (particle_weights[:, None] * particle_cylindrical_velocities).sum(
+        axis=0
+    )
+    squared_velocity_dispersion = (
+        particle_weights[:, None]
+        * (particle_cylindrical_velocities - mean_velocity) ** 2
+    ).sum(axis=0)
+
+    return np.sqrt(squared_velocity_dispersion)
+
+
+def get_cylindrical_velocity_dispersion_vector_mass_weighted(
+    particle_masses: unyt.unyt_array,
+    particle_cylindrical_velocities: unyt.unyt_array,
+) -> unyt.unyt_array:
+    """
+    Compute the mass-weighted velocity dispersion along the radial, azimuthal and vertical
+    directions for the input particles.
+
+    Parameters:
+     - particle_masses: unyt.unyt_array
+       Mass of each particle.
+     - particle_cylindrical_velocities: unyt.unyt_array
+       Velocity of each particle in a cylindrical coordinate system.
+
+    Returns a 3 element vector containing [sigma_r, sigma_phi, sigma_z].
+    """
+    mass_weights = particle_masses / particle_masses.sum()
+    return get_weighted_cylindrical_velocity_dispersion_vector(
+        mass_weights, particle_cylindrical_velocities
+    )
+
+
+def get_cylindrical_velocity_dispersion_vector_luminosity_weighted(
+    particle_luminosities: unyt.unyt_array,
+    particle_cylindrical_velocities: unyt.unyt_array,
+) -> unyt.unyt_array:
+    """
+    Compute the luminosity-weighted velocity dispersion along the radial, azimuthal and vertical
+    directions for the input particles. We return a vector for each luminosity band.
+
+    Parameters:
+     - particle_luminosities: unyt.unyt_array
+       Luminosity of each particle in different luminosity bands.
+     - particle_cylindrical_velocities: unyt.unyt_array
+       Velocity of each particle in a cylindrical coordinate system, which varies between different
+       luminosity bands.
+
+    Returns a 3 element vector for each luminosity band, which contains [sigma_r, sigma_phi, sigma_z].
+    The velocity dispersion vectors for each band are appended to the same vector, hence the shape is
+    (number_luminosity_bands, 3).
+    """
+
+    number_luminosity_bands = particle_luminosities.shape[1]
+    velocity_dispersion_vectors = (
+        np.zeros((number_luminosity_bands, 3)) * particle_cylindrical_velocities.units
+    )
+
+    for i_band, (
+        particle_luminosities_i_band,
+        particle_cylindrical_velocities_i_band,
+    ) in enumerate(zip(particle_luminosities.T, particle_cylindrical_velocities)):
+        luminosity_weights = (
+            particle_luminosities_i_band / particle_luminosities_i_band.sum()
+        )
+        velocity_dispersion_vectors[i_band] = (
+            get_weighted_cylindrical_velocity_dispersion_vector(
+                luminosity_weights, particle_cylindrical_velocities_i_band
+            )
+        )
+
+    return velocity_dispersion_vectors
+
+
 def get_angular_momentum(
     mass: unyt.unyt_array,
     position: unyt.unyt_array,
diff --git a/SOAP/property_table.py b/SOAP/property_table.py
index 9666b4ac..f10b4f99 100644
--- a/SOAP/property_table.py
+++ b/SOAP/property_table.py
@@ -1657,7 +1657,7 @@ class PropertyTable:
             shape=1,
             dtype=np.float32,
             unit="snap_length/snap_time",
-            description="Mass weighted mean rotational velocity of the stars, in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.",
+            description="Mass-weighted mean rotational velocity of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -1673,7 +1673,7 @@ class PropertyTable:
             shape=1,
             dtype=np.float32,
             unit="snap_length/snap_time",
-            description="One-dimensional velocity dispersion of the stars computed in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.",
+            description="One-dimensional mass-weighted velocity dispersion of the star, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -1689,7 +1689,7 @@ class PropertyTable:
             shape=1,
             dtype=np.float32,
             unit="snap_length/snap_time",
-            description="Velocity dispersion perpendicular to the orbital plane of the stars, computed in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.",
+            description="Mass-weighted velocity dispersion perpendicular to the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -1705,7 +1705,7 @@ class PropertyTable:
             shape=1,
             dtype=np.float32,
             unit="snap_length/snap_time",
-            description="Total velocity dispersion in the orbital plane of the stars, computed in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.",
+            description="Mass-weighted total velocity dispersion in the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -1716,6 +1716,70 @@ class PropertyTable:
             output_physical=True,
             a_scale_exponent=0,
         ),
+        "StellarRotationalVelocityLuminosityWeighted": Property(
+            name="StellarRotationalVelocityLuminosityWeighted",
+            shape=9,  # GAMA bands
+            dtype=np.float32,
+            unit="snap_length/snap_time",
+            description="Luminosity-weighted mean rotational velocity of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=[
+                "PartType4/Coordinates",
+                "PartType4/Luminosities",
+                "PartType4/Velocities",
+            ],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
+        "StellarCylindricalVelocityDispersionLuminosityWeighted": Property(
+            name="StellarCylindricalVelocityDispersionLuminosityWeighted",
+            shape=9,  # GAMA bands
+            dtype=np.float32,
+            unit="snap_length/snap_time",
+            description="One-dimensional luminosity-weighted velocity dispersion of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=[
+                "PartType4/Coordinates",
+                "PartType4/Luminosities",
+                "PartType4/Velocities",
+            ],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
+        "StellarCylindricalVelocityDispersionVerticalLuminosityWeighted": Property(
+            name="StellarCylindricalVelocityDispersionVerticalLuminosityWeighted",
+            shape=9,  # GAMA bands
+            dtype=np.float32,
+            unit="snap_length/snap_time",
+            description="Luminosity-weighted velocity dispersion perpendicular to the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=[
+                "PartType4/Coordinates",
+                "PartType4/Luminosities",
+                "PartType4/Velocities",
+            ],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
+        "StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted": Property(
+            name="StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted",
+            shape=9,  # GAMA bands
+            dtype=np.float32,
+            unit="snap_length/snap_time",
+            description="Luminosity-weighted total velocity dispersion in the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=[
+                "PartType4/Coordinates",
+                "PartType4/Luminosities",
+                "PartType4/Velocities",
+            ],
+            output_physical=True,
+            a_scale_exponent=0,
+        ),
         "MaximumStellarBirthDensity": Property(
             name="MaximumStellarBirthDensity",
             shape=1,
diff --git a/parameter_files/COLIBRE_HYBRID.yml b/parameter_files/COLIBRE_HYBRID.yml
index a8de3712..0f9b3f6e 100644
--- a/parameter_files/COLIBRE_HYBRID.yml
+++ b/parameter_files/COLIBRE_HYBRID.yml
@@ -236,6 +236,14 @@ ApertureProperties:
     StellarMassFractionInMetals: true
     StellarMassFractionInOxygen: true
     StellarVelocityDispersionMatrix: general
+    StellarRotationalVelocity: false
+    StellarCylindricalVelocityDispersion: false
+    StellarCylindricalVelocityDispersionVertical: false
+    StellarCylindricalVelocityDispersionDiscPlane: false
+    StellarRotationalVelocityLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false
     TotalMass: true
     TotalSNIaRate: true
     GasMassInColdDenseDiffuseMetals:
@@ -730,6 +738,10 @@ SubhaloProperties:
     StellarCylindricalVelocityDispersion: general
     StellarCylindricalVelocityDispersionVertical: false
     StellarCylindricalVelocityDispersionDiscPlane: false
+    StellarRotationalVelocityLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false
 aliases:
   PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent
   PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent
diff --git a/parameter_files/COLIBRE_THERMAL.yml b/parameter_files/COLIBRE_THERMAL.yml
index 8b2ddb4e..403267a5 100644
--- a/parameter_files/COLIBRE_THERMAL.yml
+++ b/parameter_files/COLIBRE_THERMAL.yml
@@ -236,6 +236,14 @@ ApertureProperties:
     StellarMassFractionInMetals: true
     StellarMassFractionInOxygen: true
     StellarVelocityDispersionMatrix: general
+    StellarCylindricalVelocityDispersion: false
+    StellarCylindricalVelocityDispersionVertical: false
+    StellarCylindricalVelocityDispersionDiscPlane: false
+    StellarCylindricalVelocityDispersionLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false
+    StellarRotationalVelocity: false
+    StellarRotationalVelocityLuminosityWeighted: false
     TotalMass: true
     TotalSNIaRate: true
     GasMassInColdDenseDiffuseMetals:
@@ -730,6 +738,10 @@ SubhaloProperties:
     StellarCylindricalVelocityDispersion: general
     StellarCylindricalVelocityDispersionVertical: false
     StellarCylindricalVelocityDispersionDiscPlane: false
+    StellarRotationalVelocityLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false
+    StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false
 aliases:
   PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent
   PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent