diff --git a/tracking/element.py b/tracking/element.py
index 34fdc76fb0f49d39c9cbfe24bd65e719cb95e6b1..911817a0c8516cf3057a10e4e3add3914cf5ef84 100644
--- a/tracking/element.py
+++ b/tracking/element.py
@@ -203,8 +203,8 @@ class TransverseMap(Element):
         
 class SkewQuadrupole:
     """
-    Skew quadrupole element used to introduce betatron coupling. The length of
-    the quadrupole is neglected.
+    Thin skew quadrupole element used to introduce betatron coupling (the 
+    length of the quadrupole is neglected).
     
     Parameters
     ----------
@@ -226,16 +226,7 @@ class SkewQuadrupole:
         ----------
         bunch : Bunch or Beam object
         """
-        matrix = np.zeros((4,4,len(bunch)))
-        matrix[0,0,:] = 1
-        matrix[1,1,:] = 1
-        matrix[2,2,:] = 1
-        matrix[3,3,:] = 1
-        matrix[1,2,:] = -1*self.strength
-        matrix[3,0,:] = -1*self.strength
         
-        xp = bunch['xp'] + matrix[1,2,:]*bunch['y']
-        yp = bunch['yp'] + matrix[3,0,:]*bunch['x']
-    
-        bunch['xp'] = xp
-        bunch['yp'] = yp
+        bunch['xp'] = bunch['xp'] - self.strength * bunch['y']
+        bunch['yp'] = bunch['yp'] - self.strength * bunch['x']
+
diff --git a/tracking/monitors/monitors.py b/tracking/monitors/monitors.py
index 88ae25859e1b501e4d648730062eb478d47bf981..32f0af252a5876a107ed43cd3c33435cbc1c0bc9 100644
--- a/tracking/monitors/monitors.py
+++ b/tracking/monitors/monitors.py
@@ -16,7 +16,7 @@ from mbtrack2.tracking.element import Element
 from mbtrack2.tracking.particles import Bunch, Beam
 from scipy.interpolate import interp1d
 from abc import ABCMeta
-#from mpi4py import MPI
+from mpi4py import MPI
 
 class Monitor(Element, metaclass=ABCMeta):
     """
@@ -269,10 +269,10 @@ class BunchMonitor(Monitor):
         group_name = "BunchData_" + str(self.bunch_number)
         dict_buffer = {"mean":(6, buffer_size), "std":(6, buffer_size),
                      "emit":(3, buffer_size), "current":(buffer_size,),
-                     "action":(2, buffer_size)}
+                     "cs_invariant":(2, buffer_size)}
         dict_file = {"mean":(6, total_size), "std":(6, total_size),
                      "emit":(3, total_size), "current":(total_size,),
-                     "action":(2, total_size)}
+                     "cs_invariant":(2, total_size)}
         self.monitor_init(group_name, save_every, buffer_size, total_size,
                           dict_buffer, dict_file, file_name, mpi_mode)
         
diff --git a/tracking/monitors/plotting.py b/tracking/monitors/plotting.py
index bbd39ba21c4f35d367d84351cb858ada379b4094..a457a39865874290f81666bfef337107fb2003a8 100644
--- a/tracking/monitors/plotting.py
+++ b/tracking/monitors/plotting.py
@@ -123,7 +123,7 @@ def plot_bunchdata(filename, bunch_number, dataset, dimension="x", x_var="time")
     bunch_number : int
         Bunch to plot. This has to be identical to 'bunch_number' parameter in 
         'BunchMonitor' object.
-    dataset : {"current", "emit", "mean", "std", "action"}
+    dataset : {"current", "emit", "mean", "std", "cs_invariant"}
         HDF5 file's dataset to be plotted.
     dimension : str, optional
         The dimension of the dataset to plot. Use "None" for "current",
@@ -175,7 +175,7 @@ def plot_bunchdata(filename, bunch_number, dataset, dimension="x", x_var="time")
         
         label = label_list[axis_index]
         
-    elif dataset == "action":
+    elif dataset == "cs_invariant":
         dimension_dict = {"x":0, "y":1}
         axis_index = dimension_dict[dimension]
         y_var = file[group][dataset][axis_index]
diff --git a/tracking/particles.py b/tracking/particles.py
index c9ad5968da65fc520f47e78c46e82057dc199f1b..f662f679d21b92a0b303f25379d19cb305108c94 100644
--- a/tracking/particles.py
+++ b/tracking/particles.py
@@ -232,9 +232,9 @@ class Bunch:
         return np.array([emitX, emitY, emitS])
     
     @property
-    def action(self):
+    def cs_invariant(self):
         """
-        Return the average action (Couran-Snyder invariant) of each plane.
+        Return the average Courant-Snyder invariant of each plane.
 
         """
         Jx = (self.ring.optics.local_gamma[0] * self['x']**2) + \