diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 23c2de6db6e1089b61138a192fdd7a771b9f8938..dbf93a6c1ec77b283839062c92b40a81dc6607ca 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -1032,6 +1032,8 @@ class ProportionalIntegralLoop():
     
     The ProportionalIntegralLoop should be initialized only after generator 
     parameters are set.
+    The method init_Ig2Vg_matrix must be called after each parameter change of 
+    cav_res.
     
     The loop must be added to the CavityResonator object using:
         cav_res.feedback.append(loop)
@@ -1068,19 +1070,11 @@ class ProportionalIntegralLoop():
         True is recommended to prevent a Vc drop in the beginning of the tracking.
         In case of small Pgain (QL ~ 1e4), Vc drop may cause baem loss due to static Robinson.
         Default is True.
-    matrix : bool, optional
-        If True, Ig2Vg_matrix is used in Ig2Vg.
-        Ig2Vg_matrix is faster but init_Ig2Vg_matrix must be called after each 
-        parameter change of cav_res.
-        If False, Ig2Vg_formula is used in Ig2Vg.
-        Default is False.
         
     Methods
     -------
     track()
         Tracking method for the Cavity PI control feedback.
-    Ig2Vg_formula()
-        Return Vg from Ig.
     init_Ig2Vg_matrix()
         Initialize matrix for Ig2Vg_matrix.
     Ig2Vg_matrix()
@@ -1125,7 +1119,7 @@ class ProportionalIntegralLoop():
     
     """
     def __init__(self, ring, cav_res, gain, sample_num, every, delay, 
-                 Vref=None, FF=True, matrix=False):
+                 Vref=None, FF=True):
         self.ring = ring
         self.cav_res = cav_res
         self.Ig2Vg_mat = np.zeros((self.ring.h, self.ring.h), dtype=complex)
@@ -1133,7 +1127,6 @@ class ProportionalIntegralLoop():
         self.FF = FF
         if not self.FF:
             self.FFconst = 0
-        self.matrix = matrix
 
         if delay > 0:
             self.delay = int(delay)
@@ -1163,8 +1156,7 @@ class ProportionalIntegralLoop():
         self.sample_list = range(0, self.ring.h, self.every)
 
         # Pre caclulation for Ig2Vg
-        if self.matrix:
-            self.init_Ig2Vg_matrix()
+        self.init_Ig2Vg_matrix()
 
     def track(self, apply_changes=True):
         """
@@ -1221,19 +1213,7 @@ class ProportionalIntegralLoop():
         self.VrefRot = self._Vref*self.Rot
         if self.FF:
             self.FFconst = np.mean(self.ig_phasor)
-            
-    def Ig2Vg_formula(self):
-        """
-        Return Vg from Ig.
-
-        Eq.26 of ref [2].
-        """
-        generator_phasor_record = np.zeros_like(self.cav_res.generator_phasor_record, dtype=complex)
-        for index in range(self.ring.h):
-            generator_phasor_record[index] = (self.cav_res.generator_phasor_record[index-1]
-                    *np.exp(-1/self.cav_res.filling_time*(1 - 1j*np.tan(self.cav_res.psi))*self.ring.T1)
-                    + self.ig_phasor_record[index]*self.cav_res.loss_factor*self.ring.T1)
-        return generator_phasor_record
+        self.init_Ig2Vg_matrix()
             
     def init_Ig2Vg_matrix(self):
         """
@@ -1251,7 +1231,7 @@ class ProportionalIntegralLoop():
     def Ig2Vg_matrix(self):
         """
         Return Vg from Ig using matrix formalism.
-        Faster but self.init_Ig2Vg should be called after each CavityResonator 
+        Warning: self.init_Ig2Vg should be called after each CavityResonator 
         parameter change.
         """
         generator_phasor_record = (self.Ig2Vg_vec*self.cav_res.generator_phasor_record[-1] +
@@ -1265,10 +1245,7 @@ class ProportionalIntegralLoop():
         Apply new values to cav_res.generator_phasor_record, cav_res.Vg and
         cav_res.theta_g from ig_phasor_record.
         """
-        if self.matrix:
-            self.cav_res.generator_phasor_record = self.Ig2Vg_matrix()
-        else:
-            self.cav_res.generator_phasor_record = self.Ig2Vg_formula()
+        self.cav_res.generator_phasor_record = self.Ig2Vg_matrix()
         self.cav_res.Vg = np.mean(np.abs(self.cav_res.generator_phasor_record))
         self.cav_res.theta_g = np.mean(np.angle(self.cav_res.generator_phasor_record))