diff --git a/vlasov/beamloading.py b/vlasov/beamloading.py
index 07688a10ae62b3cf31e6cfe0cef6cdf15f6c930d..a91b5b22b2405942463d7d0a191af2e978584731 100644
--- a/vlasov/beamloading.py
+++ b/vlasov/beamloading.py
@@ -88,6 +88,13 @@ class BeamLoadingVlasov():
             delta += cavity.Vb(self.I0) * self.F[i] * np.cos(cavity.psi - self.PHI[i])
             delta -= cavity.Vg * np.cos(cavity.theta_g)
         return delta
+    
+    def center_of_mass(self):
+        """Return center of mass position in [s]"""
+        z0 = np.linspace(self.B1, self.B2, 1000)
+        rho = self.rho(z0)
+        CM = np.average(z0, weights=rho)
+        return CM/c
 
     def u(self, z):
         """Scaled potential u(z)"""