diff --git a/mbtrack2/impedance/wakefield.py b/mbtrack2/impedance/wakefield.py
index 31712518626ca947d6faedfee0219d01d1a90876..0976f0c59983bd2297f8eeb23080cf82d43c4f2c 100644
--- a/mbtrack2/impedance/wakefield.py
+++ b/mbtrack2/impedance/wakefield.py
@@ -794,7 +794,7 @@ class WakeField:
         """
         Return an array of the impedance component names for the element.
         """
-        valid = ["Zlong", "Zxdip", "Zydip", "Zxquad", "Zyquad"]
+        valid = ["Zlong", "Zxdip", "Zydip", "Zxquad", "Zyquad", "Zxcst", "Zycst"]
         return np.array([comp for comp in dir(self) if comp in valid])
     
     @property
@@ -802,7 +802,7 @@ class WakeField:
         """
         Return an array of the wake function component names for the element.
         """
-        valid = ["Wlong", "Wxdip", "Wydip", "Wxquad", "Wyquad"]
+        valid = ["Wlong", "Wxdip", "Wydip", "Wxquad", "Wyquad", "Wxcst", "Wycst"]
         return np.array([comp for comp in dir(self) if comp in valid])
     
     @property
@@ -810,8 +810,8 @@ class WakeField:
         """
         Return an array of the component names for the element.
         """
-        valid = ["Wlong", "Wxdip", "Wydip", "Wxquad", "Wyquad", 
-                 "Zlong", "Zxdip", "Zydip", "Zxquad", "Zyquad"]
+        valid = ["Wlong", "Wxdip", "Wydip", "Wxquad", "Wyquad", "Wxcst", "Wycst",
+                 "Zlong", "Zxdip", "Zydip", "Zxquad", "Zyquad", "Zxcst", "Zycst"]
         return np.array([comp for comp in dir(self) if comp in valid])
     
     def drop(self, component):
diff --git a/mbtrack2/tracking/particles.py b/mbtrack2/tracking/particles.py
index e663ee0ac563c871951452098a1d4565a1226237..5c900bac11c0eff3566ce4791269d804a3148d3a 100644
--- a/mbtrack2/tracking/particles.py
+++ b/mbtrack2/tracking/particles.py
@@ -235,12 +235,13 @@ class Bunch:
         """
         Return the bunch emittance for each plane.
         """
-        emitX = (np.mean(self['x']**2)*np.mean(self['xp']**2) - 
-                 np.mean(self['x']*self['xp'])**2)**(0.5)
-        emitY = (np.mean(self['y']**2)*np.mean(self['yp']**2) - 
-                 np.mean(self['y']*self['yp'])**2)**(0.5)
-        emitS = (np.mean(self['tau']**2)*np.mean(self['delta']**2) - 
-                 np.mean(self['tau']*self['delta'])**2)**(0.5)
+        cor = np.squeeze([[self[name] - self[name].mean()] for name in self])
+        emitX = np.sqrt(np.mean(cor[0]**2)*np.mean(cor[1]**2) - 
+                        np.mean(cor[0]*cor[1])**2)
+        emitY = np.sqrt(np.mean(cor[2]**2)*np.mean(cor[3]**2) - 
+                        np.mean(cor[2]*cor[3])**2)
+        emitS = np.sqrt(np.mean(cor[4]**2)*np.mean(cor[5]**2) - 
+                        np.mean(cor[4]*cor[5])**2)
         return np.array([emitX, emitY, emitS])
     
     @property
diff --git a/mbtrack2/tracking/wakepotential.py b/mbtrack2/tracking/wakepotential.py
index 19bee813701ab7a7d27299a7f271e2ad23c99f3b..0d30f44f57964e780279f38a255c8a8cb7e3c500 100644
--- a/mbtrack2/tracking/wakepotential.py
+++ b/mbtrack2/tracking/wakepotential.py
@@ -273,7 +273,7 @@ class WakePotential(Element):
         
         profile0 = np.interp(tau0, self.tau, self.rho, 0, 0)
         
-        if wake_type == "Wlong" or wake_type == "Wxquad" or wake_type == "Wyquad":
+        if wake_type == "Wlong" or wake_type == "Wxquad" or wake_type == "Wyquad" or wake_type == "Wxcst" or wake_type == "Wycst":
             Wp = signal.convolve(profile0, W0*-1, mode='same')*dtau0
         elif wake_type == "Wxdip":
             dipole0 = self.dipole_moment(bunch, "x", tau0)
@@ -318,6 +318,10 @@ class WakePotential(Element):
                 elif wake_type == "Wyquad":
                     bunch["yp"] += (bunch["y"] * Wp_interp * bunch.charge 
                                     / self.ring.E0)
+                elif wake_type == "Wxcst":
+                    bunch["xp"] += Wp_interp * bunch.charge / self.ring.E0
+                elif wake_type == "Wycst":
+                    bunch["yp"] += Wp_interp * bunch.charge / self.ring.E0
                 
     def plot_last_wake(self, wake_type, plot_rho=True, plot_dipole=False, 
                        plot_wake_function=True):
@@ -346,7 +350,9 @@ class WakePotential(Element):
                   "Wxdip" : r"$W_{p,x}^{D} (V/pC)$",
                   "Wydip" : r"$W_{p,y}^{D} (V/pC)$",
                   "Wxquad" : r"$W_{p,x}^{Q} (V/pC/m)$",
-                  "Wyquad" : r"$W_{p,y}^{Q} (V/pC/m)$"}
+                  "Wyquad" : r"$W_{p,y}^{Q} (V/pC/m)$",
+                  "Wxcst" : r"$W_{p,x}^{M}$ (V/pC)",
+                  "Wycst" : r"$W_{p,y}^{M}$ (V/pC)",}
         
         Wp = getattr(self, wake_type)
         tau0 = getattr(self, "tau0_" + wake_type)