Add ProportionalIntegralLoop and DirectFeedback

Add ProportionalIntegralLoop and DirectFeedback based on a code refactoring of the RF-FB branch of Naoto.
Small changes to ProportionalLoop and TunerLoop to store the cav_res as class attribute to agree with other FBs.
Small change to the CavityResonator.feedback interface.

mbtrack2/tracking/__init__.py

@@ -8,7 +8,9 @@ from mbtrack2.tracking.synchrotron import Synchrotron
 from mbtrack2.tracking.rf import (RFCavity, 
                                   CavityResonator,
                                   ProportionalLoop,
-                                  TunerLoop)
+                                  TunerLoop,
+                                  ProportionalIntegralLoop,
+                                  DirectFeedback)
 from mbtrack2.tracking.parallel import Mpi
 from mbtrack2.tracking.element import (Element, 
                                        LongitudinalMap, 

mbtrack2/tracking/rf.py

@@ -325,7 +325,7 @@ class CavityResonator():
             
         # apply different kind of RF feedback
         for fb in self.feedback:
-            fb.track(self)
+            fb.track()
                 
         self.nturn += 1
                 
@@ -907,15 +907,15 @@ class ProportionalLoop():
     """
     Proportional feedback loop to control a CavityResonator amplitude and phase.
     
-    Feedback setpoints are cavity_resonator.Vc and cavity_resonator.theta.
+    Feedback setpoints are cav_res.Vc and cav_res.theta.
     
     The loop must be added to the CavityResonator object using:
-        cavity_resonator.feedback.append(loop)
+        cav_res.feedback.append(loop)
 
     Parameters
     ----------
     ring : Synchrotron object
-    cavity_resonator : CavityResonator
+    cav_res : CavityResonator
         The CavityResonator which is to be controlled.
     gain_A : float
         Amplitude (voltage) gain of the feedback.
@@ -925,15 +925,16 @@ class ProportionalLoop():
         Feedback delay in unit of turns.
 
     """
-    def __init__(self, ring, cavity_resonator, gain_A, gain_P, delay):
+    def __init__(self, ring, cav_res, gain_A, gain_P, delay):
         self.ring = ring
+        self.cav_res = cav_res
         self.gain_A = gain_A
         self.gain_P = gain_P
         self.delay = int(delay)
-        self.volt_delay = np.ones(self.delay)*cavity_resonator.Vc
-        self.phase_delay = np.ones(self.delay)*cavity_resonator.theta
+        self.volt_delay = np.ones(self.delay)*self.cav_res.Vc
+        self.phase_delay = np.ones(self.delay)*self.cav_res.theta
     
-    def track(self, cavity_resonator):
+    def track(self):
         """
         Tracking method for the amplitude and phase loop.
 
@@ -942,15 +943,15 @@ class ProportionalLoop():
         None.
 
         """
-        diff_A = self.volt_delay[-1] - cavity_resonator.Vc
-        diff_P = self.phase_delay[-1] - cavity_resonator.theta
-        cavity_resonator.Vg -= self.gain_A*diff_A
-        cavity_resonator.theta_g -= self.gain_P*diff_P
-        cavity_resonator.generator_phasor_record = np.ones(self.ring.h)*cavity_resonator.generator_phasor
+        diff_A = self.volt_delay[-1] - self.cav_res.Vc
+        diff_P = self.phase_delay[-1] - self.cav_res.theta
+        self.cav_res.Vg -= self.gain_A*diff_A
+        self.cav_res.theta_g -= self.gain_P*diff_P
+        self.cav_res.generator_phasor_record = np.ones(self.ring.h)*self.cav_res.generator_phasor
         self.volt_delay = np.roll(self.volt_delay, 1)
         self.phase_delay = np.roll(self.phase_delay, 1)
-        self.volt_delay[0] = cavity_resonator.cavity_voltage
-        self.phase_delay[0] = cavity_resonator.cavity_phase
+        self.volt_delay[0] = self.cav_res.cavity_voltage
+        self.phase_delay[0] = self.cav_res.cavity_phase
 
 class TunerLoop():
     """
@@ -960,12 +961,12 @@ class TunerLoop():
     Only a proportional controller is assumed.
     
     The loop must be added to the CavityResonator object using:
-        cavity_resonator.feedback.append(loop)
+        cav_res.feedback.append(loop)
 
     Parameters
     ----------
     ring : Synchrotron object
-    cavity_resonator : CavityResonator
+    cav_res : CavityResonator
         The CavityResonator which is to be controlled.
     gain : float
         Proportional gain of the tuner loop.
@@ -981,11 +982,12 @@ class TunerLoop():
         Tuning offset in [rad].
         
     """
-    def __init__(self, ring, cavity_resonator, gain=0.01, avering_period=0, 
+    def __init__(self, ring, cav_res, gain=0.01, avering_period=0, 
                  offset=0):
         self.ring = ring
+        self.cav_res = cav_res
         if avering_period == 0:
-            fs = self.ring.synchrotron_tune(cavity_resonator.Vc)*self.ring.f1/self.ring.h
+            fs = self.ring.synchrotron_tune(self.cav_res.Vc)*self.ring.f1/self.ring.h
             avering_period = 2/fs/self.ring.T0
     
         self.Pgain = gain
@@ -994,7 +996,7 @@ class TunerLoop():
         self.diff = 0
         self.count = 0
     
-    def track(self, cavity_resonator):
+    def track(self):
         """
         Tracking method for the tuner loop.
     
@@ -1005,9 +1007,501 @@ class TunerLoop():
         """
         if self.count == self.avering_period:
             diff = self.diff/self.avering_period - self.offset
-            cavity_resonator.psi -= diff * self.Pgain
+            self.cav_res.psi -= diff * self.Pgain
             self.count = 0
             self.diff = 0
         else:
-            self.diff += cavity_resonator.cavity_phase - cavity_resonator.theta_g + cavity_resonator.psi
+            self.diff += self.cav_res.cavity_phase - self.cav_res.theta_g + self.cav_res.psi
             self.count += 1
+
+class ProportionalIntegralLoop():
+    """
+    Proportional Integral (PI) loop to control a CavityResonator amplitude and 
+    phase via generator current (ig) [1,2].
+    
+    The ProportionalIntegralLoop should be initialized only after generator 
+    parameters are set.
+    
+    The loop must be added to the CavityResonator object using:
+        cav_res.feedback.append(loop)
+
+    Parameters
+    ----------
+    ring : Synchrotron object
+    cav_res : CavityResonator object
+        CavityResonator in which the loop will be added. 
+    gain : float or float list
+        Pgain and Igain of the feedback system
+        For IQ feedback, same gain set is used for I and Q.
+        In case of normal conducting cavity (QL~1e4),
+            the Pgain of ~1.0 and Igain of ~1e4(5) are usually used.
+        In case of super conducting cavity (QL > 1e6),
+            the Pgain of ~100 can be used.
+        In a "bad" parameter set, unstable oscillation of Vc is caused.
+        So, parameter scan of Vcloop should be made. 
+    sample_num : int
+        Sample number to monitor Vc 
+        The averaged Vc value in sample_num is monitored.
+        Units are in bucket numbers.
+    every : int
+        Interval to monitor and change Vc
+        Units are in bucket numbers.
+    delay : int
+        Loop delay of the PI feedback system. 
+        Units are in bucket numbers.
+    Vref : array or complex, optional
+        Target (set point) value of the feedback.
+        If None, cav_res.Vc and cav_res.theta are set as the reference.
+        Default is None.
+    FF : bool, optional
+        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()
+        Return Vg from Ig using matrix formalism.
+    Ig2Vg()
+        Go from Ig to Vg and apply values.
+    Vg2Ig(Vg)
+        Return Ig from Vg (assuming constant Vg).
+
+    Notes
+    -----
+    Assumption : delay >> every ~ sample_num
+    
+    Adjusting ig(Vg) parameter to keep the Vc constant (to target(Vref) values).
+    The "track" method is
+       1) Monitoring the Vc phasor
+       Mean Vc value between specified bunch number (sample_num) is monitored
+       with specified interval period (every).
+       2) Changing the ig phasor
+       The ig phasor is changed according the difference of the specified
+       reference values (Vref) with specified gain (gain).
+       By using ig instead of Vg, the cavity response can be taken account.
+       3) ig changes are reflected to Vg after the specifed delay (delay) of the system
+
+    Vc,IQ-->(rot)-->(-)-->(V->I,fac)-->PI-->Ig,IQ {--> Vg,IQ}
+                     |
+                     Ref
+
+    Examples
+    --------
+    PF; E0=2.5GeV, C=187m, frf=500MHz
+        QL=11800, fs0=23kHz
+        ==> gain=[0.5,1e4], sample_num=8, every=7(13ns), delay=500(1us)
+                     is one reasonable parameter set.
+     
+    References
+    ----------
+    [1] : https://en.wikipedia.org/wiki/PID_controller
+    [2] : Yamamoto, N., Takahashi, T., & Sakanaka, S. (2018). Reduction and 
+    compensation of the transient beam loading effect in a double rf system 
+    of synchrotron light sources. PRAB, 21(1), 012001.
+    
+    """
+    def __init__(self, ring, cav_res, gain, sample_num, every, delay, 
+                 Vref=None, FF=True, matrix=False):
+        self.ring = ring
+        self.cav_res = cav_res
+        self.Ig2Vg_mat = np.zeros((self.ring.h, self.ring.h), dtype=complex)
+        self.gain = gain
+        self.FF = FF
+        if not self.FF:
+            self.FFconst = 0
+        self.matrix = matrix
+
+        if delay > 0:
+            self.delay = int(delay)
+        else:
+            self.delay = 1
+        if every > 0:
+            self.every = int(every)
+        else:
+            self.every = 1
+        record_size = int(np.ceil(self.delay/self.every))
+        if record_size < 1:
+            raise ValueError("Bad parameter set : delay or sample_every")
+        self.sample_num = int(sample_num)
+
+        # init lists for FB process
+        self.ig_phasor = np.ones(self.ring.h, dtype=complex)*self.Vg2Ig(self.cav_res.generator_phasor)
+        self.ig_phasor_record = self.ig_phasor
+        self.vc_previous = np.ones(self.sample_num)*self.cav_res.cavity_phasor
+        self.diff_record = np.zeros(record_size, dtype=complex)
+        self.I_record = 0+0j
+
+        if Vref is None:
+            self.Vref = self.cav_res.Vc*np.exp(1j*self.cav_res.theta)
+        else: 
+            self.Vref = Vref
+
+        self.sample_list = range(0, self.ring.h, self.every)
+
+        # Pre caclulation for Ig2Vg
+        if self.matrix:
+            self.init_Ig2Vg_matrix()
+
+    def track(self, apply_changes=True):
+        """
+        Tracking method for the Cavity PI control feedback.
+
+        Returns
+        -------
+        None.
+
+        """
+        vc_list = np.concatenate([self.vc_previous, self.cav_res.cavity_phasor_record]) #This line is slowing down the process.
+        self.ig_phasor.fill(self.ig_phasor[-1])
+        
+        for index in self.sample_list: 
+            # 2) updating Ig using last item of the list
+            diff = self.diff_record[-1] - self.FFconst
+            self.I_record += diff/self.ring.f1
+            fb_value = self.gain[0]*diff + self.gain[1]*self.I_record
+            self.ig_phasor[index:] = self.Vg2Ig(fb_value) + self.FFconst
+            # Shift the record
+            self.diff_record = np.roll(self.diff_record, 1)
+            # 1) recording diff as a first item of the list
+            mean_vc = np.mean(vc_list[index:self.sample_num+index])*self.Rot
+            self.diff_record[0] = self.VrefRot - mean_vc
+        # update sample_list for next turn
+        self.sample_list = range(index+self.every-self.ring.h,self.ring.h,self.every)
+        # update vc_previous for next turn
+        self.vc_previous = self.cav_res.cavity_phasor_record[- self.sample_num:]
+        
+        self.ig_phasor_record = self.ig_phasor
+        
+        if apply_changes:
+            self.Ig2Vg()
+        
+    @property
+    def Vref(self):
+        """Return the target (set point) value of the feedback as a complex."""
+        return self._Vref
+
+    @Vref.setter
+    def Vref(self, value):
+        """
+        Set the target (set point) value of the feedback.
+
+        Vref can be specified as an IQ-complex style or array-like [Amp, Phase].
+        If not specified, cav_res.Vc and cav_res.theta are used.
+
+        """
+        if isinstance(value,complex):
+            self._Vref = value
+        else:
+            self._Vref = value[0]*np.exp(1j*value[1])
+        self.Rot = np.exp(-1j*np.angle(self._Vref))
+        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
+            
+    def init_Ig2Vg_matrix(self):
+        """
+        Initialize matrix for Ig2Vg_matrix.
+        
+        Shoud be called before first use of Ig2Vg_matrix and after each cavity 
+        parameter change.
+        """
+        k = np.arange(0, self.ring.h)
+        self.Ig2Vg_vec = np.exp(-1/self.cav_res.filling_time*(1 - 1j*np.tan(self.cav_res.psi))*self.ring.T1*(k+1))
+        tempV = np.exp(-1/self.cav_res.filling_time*self.ring.T1*k*(1 - 1j*np.tan(self.cav_res.psi)))
+        for idx in np.arange(self.ring.h):
+            self.Ig2Vg_mat[idx:,idx] = tempV[:self.ring.h-idx]
+            
+    def Ig2Vg_matrix(self):
+        """
+        Return Vg from Ig using matrix formalism.
+        Faster but self.init_Ig2Vg should be called after each CavityResonator 
+        parameter change.
+        """
+        generator_phasor_record = (self.Ig2Vg_vec*self.cav_res.generator_phasor_record[-1] +
+                    self.Ig2Vg_mat.dot(self.ig_phasor_record)*self.cav_res.loss_factor*self.ring.T1)
+        return generator_phasor_record
+            
+    def Ig2Vg(self):
+        """
+        Go from Ig to Vg.
+        
+        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.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))
+        
+    def Vg2Ig(self, Vg):
+        """
+        Return Ig from Vg (assuming constant Vg).
+
+        Eq.25 of ref [2] assuming the dVg/dt = 0.
+        """
+        fac = 1/self.cav_res.filling_time/self.cav_res.loss_factor # 1/R_L
+        ig = fac*Vg*( 1 - 1j*np.tan(self.cav_res.psi) )
+        return ig
+
+class DirectFeedback(ProportionalIntegralLoop):
+    """
+    Direct RF FB (DFB) via generator current (ig) using PI controller [1].
+    
+    The DirectFeedback inherits from ProportionalIntegralLoop so all 
+    mandatory parameters from ProportionalIntegralLoop should be passed as
+    kwargs when initializing a DirectFeedback object.
+        
+    To avoid cavity-beam unmatching (large synchrotron oscilation of beam),
+    the CavityResonator generator parameters should be set before 
+    initialization.
+    
+    Parameters
+    ----------
+    DFB_gain : float
+        Gain of the DFB.
+    DFB_phase_shift : float
+        Phase shift of the DFB.
+    DFB_sample_num : int, optional
+        Sample number to monitor Vc for the DFB.
+        The averaged Vc value in DFB_sample_num is monitored.
+        Units are in bucket numbers.
+        If None, value from the ProportionalIntegralLoop is used.
+        Default is None.
+    DFB_every : int, optional
+        Interval to monitor and change Vc for the DFB.
+        Units are in bucket numbers.
+        If None, value from the ProportionalIntegralLoop is used.
+        Default is None.
+    DFB_delay : int, optional
+        Loop delay of the DFB.
+        Units are in bucket numbers.
+        If None, value from the ProportionalIntegralLoop is used.
+        Default is None.
+        
+    Attributes
+    ----------
+    phase_shift : float
+        Phase shift applied to DFB, defined as psi - DFB_phase_shift.
+    DFB_psi: float
+        Return detuning angle value with Direct RF feedback in [rad].
+    DFB_alpha : float
+        Return the amplitude ratio alpha of the DFB.
+    DFB_gamma : float
+        Return the gain gamma of the DFB.
+    DFB_Rs : float
+        Return the shunt impedance of the DFB in [ohm].
+    
+    Methods
+    -------
+    DFB_parameter_set(DFB_gain, DFB_phase_shift)
+        Set DFB gain and phase shift parameters.
+    track()
+        Tracking method for the Direct RF feedback.
+    DFB_Vg()
+        Return the generator voltage main and DFB components in [V].
+    DFB_fs()
+        Return the modified synchrotron frequency in [Hz].
+        
+    References
+    ----------
+    [1] : Akai, K. (2022). Stability analysis of rf accelerating mode with 
+    feedback loops under heavy beam loading in SuperKEKB. PRAB, 25(10), 
+    102002.
+
+    """
+    def __init__(self, DFB_gain, DFB_phase_shift, DFB_sample_num=None,
+                 DFB_every=None, DFB_delay=None, **kwargs):
+        super(DirectFeedback, self).__init__(**kwargs)
+        
+        if DFB_delay is not None:
+            self.DFB_delay = int(DFB_delay)
+        else:
+            self.DFB_delay = self.delay
+
+        if DFB_sample_num is not None:
+            self.DFB_sample_num = int(DFB_sample_num)
+        else:
+            self.DFB_sample_num = self.sample_num
+            
+        if DFB_every is not None:
+            self.DFB_every = int(DFB_every)
+        else:
+            self.DFB_every = self.every
+
+        record_size = int(np.ceil(self.DFB_delay/self.DFB_every))
+        if record_size < 1:
+            raise ValueError("Bad parameter set : delay or sample_every")
+
+        self.DFB_parameter_set(DFB_gain, DFB_phase_shift)
+        if np.sum(np.abs(self.cav_res.beam_phasor)) == 0:
+            cavity_phasor = self.cav_res.Vc*np.exp(1j*self.cav_res.theta)
+        else:
+            cavity_phasor = np.mean(self.cav_res.cavity_phasor_record)
+        self.DFB_VcRecord = np.ones(record_size, dtype=complex)*cavity_phasor
+        self.DFB_vc_previous = np.ones(self.DFB_sample_num, dtype=complex)*cavity_phasor 
+
+        self.DFB_sample_list = range(0, self.ring.h, self.DFB_every)
+        
+    @property
+    def DFB_phase_shift(self):
+        """Return DFB phase shift."""
+        return self._DFB_phase_shift
+    
+    @DFB_phase_shift.setter
+    def DFB_phase_shift(self, value):
+        """Set DFB_phase_shift and phase_shift"""
+        self._DFB_phase_shift = value
+        self._phase_shift = self.cav_res.psi - value
+       
+    @property
+    def phase_shift(self):
+        """
+        Return phase shift applied to DFB.
+        Defined as self.cav_res.psi - self.DFB_phase_shift.
+        """
+        return self._phase_shift
+    
+    @property
+    def DFB_psi(self):
+        """
+        Return detuning angle value with Direct RF feedback in [rad].
+    
+        Fig.4 of ref [1].
+        """
+        return (np.angle(np.mean(self.cav_res.cavity_phasor_record)) - 
+                np.angle(np.mean(self.ig_phasor_record)))
+    
+    @property
+    def DFB_alpha(self):
+        """
+        Return the amplitude ratio alpha of the DFB.
+        
+        Near Eq. (13) of [1].
+        """
+        fac = np.abs(np.mean(self.DFB_ig_phasor)/np.mean(self.ig_phasor_record))
+        return 20*np.log10(fac)
+    
+    @property
+    def DFB_gamma(self):
+        """
+        Return the gain gamma of the DFB.
+        
+        Near Eq. (13) of [1].
+        """
+        fac = np.abs(np.mean(self.DFB_ig_phasor)/np.mean(self.ig_phasor_record))
+        return fac/(1-fac)
+    
+    @property
+    def DFB_Rs(self):
+        """
+        Return the shunt impedance of the DFB in [ohm].
+        
+        Eq. (15) of [1].
+        """
+        return self.cav_res.Rs/(1+self.DFB_gamma*np.cos(self.DFB_psi))
+        
+    def DFB_parameter_set(self, DFB_gain, DFB_phase_shift):
+        """
+        Set DFB gain and phase shift parameters.
+
+        Parameters
+        ----------
+        DFB_gain : float
+            Gain of the DFB.
+        DFB_phase_shift : float
+            Phase shift of the DFB.
+
+        """
+        self.DFB_gain = DFB_gain
+        self.DFB_phase_shift = DFB_phase_shift
+
+        if np.sum(np.abs(self.cav_res.beam_phasor)) == 0:
+            vc = np.ones(self.ring.h)*self.cav_res.Vc*np.exp(1j*self.cav_res.theta)
+        else:
+            vc = self.cav_res.cavity_phasor_record
+        vg_drf = self.DFB_gain*vc*np.exp(1j*self.phase_shift)
+        self.DFB_ig_phasor = self.Vg2Ig(vg_drf)
+
+        self.ig_phasor = self.ig_phasor_record - self.DFB_ig_phasor
+        if self.FF:
+            self.FFconst = np.mean(self.ig_phasor)
+
+    def track(self):
+        """
+        Tracking method for the Direct RF feedback.
+
+        Returns
+        -------
+        None.
+
+        """
+        super(DirectFeedback, self).track(False)
+        
+        vc_list = np.concatenate([self.DFB_vc_previous, self.cav_res.cavity_phasor_record])
+        self.DFB_ig_phasor = np.roll(self.DFB_ig_phasor, 1)
+        for index in self.DFB_sample_list:
+            # 2) updating Ig using last item of the list
+            vg_drf = self.DFB_gain*self.DFB_VcRecord[-1]*np.exp(1j*self.phase_shift)
+            self.DFB_ig_phasor[index:] = self.Vg2Ig(vg_drf)
+            # Shift the record
+            self.DFB_VcRecord = np.roll(self.DFB_VcRecord, 1)
+            # 1) recording Vc
+            mean_vc = np.mean(vc_list[index:self.DFB_sample_num+index])
+            self.DFB_VcRecord[0] = mean_vc
+        # update sample_list for next turn
+        self.DFB_sample_list = range(index+self.DFB_every-self.ring.h, self.ring.h, self.DFB_every)
+        # update vc_previous for next turn
+        self.DFB_vc_previous = self.cav_res.cavity_phasor_record[- self.DFB_sample_num:]
+
+        self.ig_phasor_record = self.ig_phasor + self.DFB_ig_phasor
+        
+        self.Ig2Vg()
+    
+    def DFB_Vg(self, vc=-1):
+        """Return the generator voltage main and DFB components in [V]."""
+        if vc == -1:
+            vc = np.mean(self.cav_res.cavity_phasor_record)
+        vg_drf=self.DFB_gain*vc*np.exp(1j*self.phase_shift)
+        vg_main=np.mean(self.cav_res.generator_phasor_record)-vg_drf
+        return vg_main, vg_drf
+
+    def DFB_fs(self, vg_main=-1, vg_drf=-1):
+        """Return the modified synchrotron frequency in [Hz]."""
+        vc = np.mean(self.cav_res.cavity_phasor_record)
+        if vg_drf ==-1:
+            vg_drf=self.DFB_gain*vc*np.exp(1j*self.phase_shift)
+        if vg_main == -1:
+            vg_main=np.mean(self.cav_res.generator_phasor_record)-vg_drf
+        vg_sum = np.abs(vg_main)*np.sin(np.angle(vg_main))+np.abs(vg_drf)*np.sin(np.angle(vg_drf))
+        omega_s = 0
+        if (vg_sum) > 0.0:
+            omega_s=np.sqrt(self.ring.ac*self.ring.omega1*(vg_sum)/self.ring.E0/self.ring.T0)
+        return omega_s/2/np.pi