diff --git a/mbtrack2/instability/instabilities.py b/mbtrack2/instability/instabilities.py
index 12e2237a9a4855689b3fb3b093a7adcb5333fcf3..75988ebd269336066f2562f12aad3873b4354b52 100644
--- a/mbtrack2/instability/instabilities.py
+++ b/mbtrack2/instability/instabilities.py
@@ -96,9 +96,10 @@ def cbi_threshold(ring, I, Vrf, f, beta, Ncav=1):
def lcbi_growth_rate_mode(ring,
I,
- Vrf,
M,
mu,
+ synchrotron_tune=None,
+ Vrf=None,
fr=None,
RL=None,
QL=None,
@@ -114,12 +115,14 @@ def lcbi_growth_rate_mode(ring,
ring : Synchrotron object
I : float
Total beam current in [A].
- Vrf : float
- Total RF voltage in [V].
M : int
Nomber of bunches in the beam.
mu : int
Coupled bunch mode number (= 0, ..., M-1).
+ synchrotron_tune : float, optional
+ Synchrotron tune.
+ Vrf : float, optinal
+ Total RF voltage in [V] used to compute synchrotron tune if not given.
fr : float or list, optional
Frequency of the resonator in [Hz].
RL : float or list, optional
@@ -132,7 +135,7 @@ def lcbi_growth_rate_mode(ring,
Returns
-------
float
- Coupled bunch instability growth rate for the mode mu.
+ Coupled bunch instability growth rate for the mode mu in [s-1].
References
----------
@@ -142,7 +145,12 @@ def lcbi_growth_rate_mode(ring,
"""
- nu_s = ring.synchrotron_tune(Vrf)
+ if synchrotron_tune is None and Vrf is None:
+ raise ValueError("Either synchrotron_tune or Vrf is needed.")
+ if synchrotron_tune is None:
+ nu_s = ring.synchrotron_tune(Vrf)
+ else:
+ nu_s = synchrotron_tune
factor = ring.eta() * I / (4 * np.pi * ring.E0 * nu_s)
if isinstance(fr, (float, int)):
@@ -186,7 +194,15 @@ def lcbi_growth_rate_mode(ring,
return factor * sum_val
-def lcbi_growth_rate(ring, I, Vrf, M, fr=None, RL=None, QL=None, Z=None):
+def lcbi_growth_rate(ring,
+ I,
+ M,
+ synchrotron_tune=None,
+ Vrf=None,
+ fr=None,
+ RL=None,
+ QL=None,
+ Z=None):
"""
Compute the maximum growth rate for longitudinal coupled bunch instability
driven an impedance [1].
@@ -198,10 +214,12 @@ def lcbi_growth_rate(ring, I, Vrf, M, fr=None, RL=None, QL=None, Z=None):
ring : Synchrotron object
I : float
Total beam current in [A].
- Vrf : float
- Total RF voltage in [V].
M : int
Nomber of bunches in the beam.
+ synchrotron_tune : float, optional
+ Synchrotron tune.
+ Vrf : float, optinal
+ Total RF voltage in [V] used to compute synchrotron tune if not given.
fr : float or list, optional
Frequency of the HOM in [Hz].
RL : float or list, optional
@@ -214,12 +232,13 @@ def lcbi_growth_rate(ring, I, Vrf, M, fr=None, RL=None, QL=None, Z=None):
Returns
-------
growth_rate : float
- Maximum coupled bunch instability growth rate.
+ Maximum coupled bunch instability growth rate in [s-1].
mu : int
Coupled bunch mode number corresponding to the maximum coupled bunch
instability growth rate.
growth_rates : array
- Coupled bunch instability growth rates for the different mode numbers.
+ Coupled bunch instability growth rates for the different mode numbers
+ in [s-1].
References
----------
@@ -232,9 +251,10 @@ def lcbi_growth_rate(ring, I, Vrf, M, fr=None, RL=None, QL=None, Z=None):
for i in range(M):
growth_rates[i] = lcbi_growth_rate_mode(ring,
I,
- Vrf,
M,
i,
+ synchrotron_tune=synchrotron_tune,
+ Vrf=Vrf,
fr=fr,
RL=RL,
QL=QL,
@@ -246,7 +266,14 @@ def lcbi_growth_rate(ring, I, Vrf, M, fr=None, RL=None, QL=None, Z=None):
return growth_rate, mu, growth_rates
-def lcbi_stability_diagram(ring, I, Vrf, M, modes, cavity_list, detune_range):
+def lcbi_stability_diagram(ring,
+ I,
+ M,
+ modes,
+ cavity_list,
+ detune_range,
+ synchrotron_tune=None,
+ Vrf=None):
"""
Plot longitudinal coupled bunch instability stability diagram for a
arbitrary list of CavityResonator objects around a detuning range.
@@ -260,8 +287,6 @@ def lcbi_stability_diagram(ring, I, Vrf, M, modes, cavity_list, detune_range):
Ring parameters.
I : float
Total beam current in [A].
- Vrf : float
- Total RF voltage in [V].
M : int
Nomber of bunches in the beam.
modes : list
@@ -274,6 +299,10 @@ def lcbi_stability_diagram(ring, I, Vrf, M, modes, cavity_list, detune_range):
- mode dampers
detune_range : array
Detuning range (list of points) of the last CavityResonator object.
+ synchrotron_tune : float, optional
+ Synchrotron tune.
+ Vrf : float, optinal
+ Total RF voltage in [V] used to compute synchrotron tune if not given.
Returns
-------
@@ -290,23 +319,25 @@ def lcbi_stability_diagram(ring, I, Vrf, M, modes, cavity_list, detune_range):
for cav in cavity_list[:-1]:
fixed_gr += lcbi_growth_rate_mode(ring,
I=I,
- Vrf=Vrf,
mu=mu,
fr=cav.fr,
RL=cav.RL,
QL=cav.QL,
- M=M)
+ M=M,
+ synchrotron_tune=synchrotron_tune,
+ Vrf=Vrf)
cav = cavity_list[-1]
for i, det in enumerate(detune_range):
gr = lcbi_growth_rate_mode(ring,
I=I,
- Vrf=Vrf,
mu=mu,
fr=cav.m * ring.f1 + det,
RL=cav.RL,
QL=cav.QL,
- M=M)
+ M=M,
+ synchrotron_tune=synchrotron_tune,
+ Vrf=Vrf)
Rth[i] = (1 / ring.tau[2] - fixed_gr) * cav.Rs / gr
ax.plot(detune_range * 1e-3,
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 3d37de9b07f80a2bb01472f405e282b334ce359c..3defbc1e29c1d472b03810e6d34c0a37360a9e92 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -7,8 +7,8 @@ import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.legend_handler import HandlerPatch
-
from mbtrack2.tracking.element import Element
+from mbtrack2.instability import lcbi_growth_rate
class RFCavity(Element):
@@ -582,7 +582,7 @@ class CavityResonator():
@property
def DFB_ig_phasor(self):
- """Last current generator phasor of each bunch in [A]"""
+ """Last direct feedback current generator phasor of each bunch in [A]"""
for FB in self.feedback:
if isinstance(FB, DirectFeedback):
return FB.DFB_ig_phasor
@@ -938,59 +938,70 @@ class CavityResonator():
return fig
- def is_CBI_stable(self,I0,tune=0,mode=[0],max_freq=5,bool_return=False):
+ def is_CBI_stable(self,
+ I0,
+ synchrotron_tune=None,
+ modes=None,
+ bool_return=False):
"""
- Check Coupled-Bunch-Instability stability
+ Check Coupled-Bunch-Instability stability,
Effect of Direct RF feedback is not included.
- This method caluclates the CBI growth rate from own impedance.
+ This method is a wraper around lcbi_growth_rate to caluclate the CBI
+ growth rate from the CavityResonator own impedance.
+
+ See lcbi_growth_rate for details.
Parameters
----------
I0 : float
Beam current in [A].
- tune : float
- fractional number of longitudinal tune
- mode : float or list of float
- Coupled Bunch Instability mode number
- bool_return : bool
- if True, return bool
+ synchrotron_tune : float, optinal
+ Fractional number of longitudinal tune.
+ If None, synchrotron_tune is computed using self.Vc as total
+ voltage.
+ Default is None.
+ modes : float or list of float, optional
+ Coupled Bunch Instability mode numbers to consider.
+ If not None, return growth_rates of the input coupled bunch modes.
+ Default is None.
+ bool_return : bool, optional
+ If True:
+ - return True if the most unstable mode is stronger than
+ longitudinal damping time.
+ - return False if it is not the case.
+ Default is False.
Returns
-------
- bool
-
- """
- if tune == 0:
- tune = self.ring.synchrotron_tune(self.Vc)
- if tune > 1:
- tune = tune - int(tune)
- cof = self.ring.ac*I0/2.0/tune/self.ring.E0
- if isinstance(mode,list):
- CBImode = mode
- else:
- CBImode = [mode]
-
- gr = np.zeros(len(CBImode))
- gr_bool = np.zeros(len(CBImode),dtype=bool)
- count = 0
- for i in CBImode:
- #fp = np.array([0,1,2,3,4])*self.ring.f1+i*self.ring.f0*(i+tune)
- fp = self.ring.f0*(np.arange(max_freq)*self.ring.h+(i+tune))
- #fm = np.array([1,2,3,4,5])*self.ring.f1-i*self.ring.f0*(i+tune)
- fm = self.ring.f0*(np.arange(1,max_freq+1)*self.ring.h-(i+tune))
- sumZ = np.sum(fp*np.real(self.Z(fp))) - np.sum(fm*np.real(self.Z(fm)))
-
- gr[count] = cof*sumZ
- gr_bool[count] = ( gr[count] > 1/self.ring.tau[2] )
- count +=1
-
+ growth_rate : float
+ Maximum coupled bunch instability growth rate in [s-1].
+ mu : int
+ Coupled bunch mode number corresponding to the maximum coupled bunch
+ instability growth rate.
+
+ """
+ growth_rate, mu, growth_rates = lcbi_growth_rate(self.ring,
+ I0,
+ M=self.ring.h,
+ synchrotron_tune=synchrotron_tune,
+ Vrf=self.Vc,
+ fr=self.fr,
+ RL=self.RL,
+ QL=self.QL)
+
+ if modes is not None:
+ growth_rates = growth_rates[modes]
+ return growth_rates
+
if bool_return:
- return gr_bool
+ if growth_rate > 1/self.ring.tau[2]:
+ return True
+ else:
+ return False
else:
- return gr
-
-
+ return growth_rate, mu
+
def is_DC_Robinson_stable(self, I0):
"""
Check DC Robinson stability - Eq. (6.1.1) [1]
@@ -1280,12 +1291,16 @@ class ProportionalIntegralLoop():
Number of bunch over which the mean cavity voltage is computed.
Units are in bucket numbers.
every : int
- Sampling and Clock period of the feedback controller
+ Sampling and clock period of the feedback controller
Time interval between two cavity voltage monitoring and feedback.
Units are in bucket numbers.
delay : int
Loop delay of the PI feedback system.
Units are in bucket numbers.
+ IIR_cutoff : float, optinal
+ Cutoff frequency of the IIR filter in [Hz].
+ If 0, cutoff frequency is infinity.
+ Default is 0.
FF : bool, optional
Boolean switch to use feedforward constant.
True is recommended to prevent a cavity voltage drop in the beginning
@@ -1308,6 +1323,12 @@ class ProportionalIntegralLoop():
Go from Ig to Vg and apply values.
Vg2Ig(Vg)
Return Ig from Vg (assuming constant Vg).
+ IIR_init(cutoff)
+ Initialization for the IIR filter.
+ IIR(input)
+ Return IIR filter output.
+ IIRcutoff()
+ Return IIR cutoff frequency in [Hz].
Notes
-----
@@ -1324,9 +1345,9 @@ class ProportionalIntegralLoop():
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-->(rot)-->(-)-->(V->I,fac)-->PI-->Ig --> Vg
- |
- Ref
+ Vc-->(rot)-->IIR-->(-)-->(V->I,fac)-->PI-->Ig --> Vg
+ |
+ Ref
Examples
--------
@@ -1345,7 +1366,15 @@ class ProportionalIntegralLoop():
of synchrotron light sources. PRAB, 21(1), 012001.
"""
- def __init__(self, ring, cav_res, gain, sample_num, every, delay, IIR_cutoff=0,FF=True):
+ def __init__(self,
+ ring,
+ cav_res,
+ gain,
+ sample_num,
+ every,
+ delay,
+ IIR_cutoff=0,
+ FF=True):
self.ring = ring
self.cav_res = cav_res
self.Ig2Vg_mat = np.zeros((self.ring.h, self.ring.h), dtype=complex)
@@ -1377,8 +1406,6 @@ class ProportionalIntegralLoop():
self.sample_list = range(0, self.ring.h, self.every)
- # IIR filter
- ## IIR_cutoff ; cutoff frequecny in Hz. If 0, cutoff frequency is infinity.
self.IIR_init(IIR_cutoff)
self.init_FFconst()
@@ -1478,12 +1505,16 @@ class ProportionalIntegralLoop():
"""
return Vg * (1 - 1j * np.tan(self.cav_res.psi)) / self.cav_res.RL
- def IIR_init(self,cutoff):
+ def IIR_init(self, cutoff):
"""
- f = 1/2/pi/T arccos((2-2r-r*r)/2(1-r))
- T: Sample time
- r: IIRcoef
- f: cutoff freqeuncy
+ Initialization for the IIR filter.
+
+ Parameters
+ ----------
+ cutoff : float
+ Cutoff frequency of the IIR filter in [Hz].
+ If 0, cutoff frequency is infinity.
+
"""
if cutoff == 0:
self.IIRcoef = 1.0
@@ -1492,20 +1523,20 @@ class ProportionalIntegralLoop():
T = self.ring.T1 * self.every
alpha = np.cos(omega*T)-1
tmp = alpha * alpha - 2 * alpha
- print(tmp)
if tmp > 0:
self.IIRcoef = alpha + np.sqrt(tmp)
else:
self.IIRcoef = T * cutoff * 2 * np.pi
self.IIRout = self.cav_res.Vc
- def IIR(self,input):
- self.IIRout = (1-self.IIRcoef)*self.IIRout+self.IIRcoef*input
+ def IIR(self, input):
+ """Return IIR filter output."""
+ self.IIRout = (1 - self.IIRcoef)*self.IIRout + self.IIRcoef*input
return self.IIRout
@property
def IIRcutoff(self):
- """Return IIR cutoff frequency."""
+ """Return IIR cutoff frequency in [Hz]."""
T = self.ring.T1 * self.every
return 1.0/2.0/np.pi/T * np.arccos((2-2*self.IIRcoef-self.IIRcoef*self.IIRcoef)/2/(1-self.IIRcoef))