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/monitors/monitors.py b/mbtrack2/tracking/monitors/monitors.py
index 69c62ae91f832551be3a31c8dacdebea076151be..01794ebea901beaf33c956b94d084616e26d55f4 100644
--- a/mbtrack2/tracking/monitors/monitors.py
+++ b/mbtrack2/tracking/monitors/monitors.py
@@ -1567,6 +1567,10 @@ class CavityMonitor(Monitor):
ring.h,
buffer_size,
),
+ "DFB_ig_phasor":(
+ ring.h,
+ buffer_size,
+ ),
"detune": (buffer_size, ),
"psi": (buffer_size, ),
"Vg": (buffer_size, ),
@@ -1595,6 +1599,10 @@ class CavityMonitor(Monitor):
ring.h,
total_size,
),
+ "DFB_ig_phasor":(
+ ring.h,
+ total_size,
+ ),
"detune": (total_size, ),
"psi": (total_size, ),
"Vg": (total_size, ),
@@ -1611,6 +1619,7 @@ class CavityMonitor(Monitor):
"beam_phasor_record": complex,
"generator_phasor_record": complex,
"ig_phasor_record": complex,
+ "DFB_ig_phasor":complex,
"detune": float,
"psi": float,
"Vg": float,
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 6e731c5c8eb1aa8b4699e5aa26b52dfbf3966865..c01e63c3e98142e30a186c8007d68351dbfa0981 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):
@@ -579,6 +579,14 @@ class CavityResonator():
if isinstance(FB, (ProportionalIntegralLoop, DirectFeedback)):
return FB.ig_phasor_record
return np.zeros(self.ring.h)
+
+ @property
+ def DFB_ig_phasor(self):
+ """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
+ return np.zeros(self.ring.h)
@property
def cavity_voltage(self):
@@ -930,6 +938,70 @@ class CavityResonator():
return fig
+ def is_CBI_stable(self,
+ I0,
+ synchrotron_tune=None,
+ modes=None,
+ bool_return=False):
+ """
+ Check Coupled-Bunch-Instability stability,
+ Effect of Direct RF feedback is not included.
+
+ 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].
+ 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
+ -------
+ 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:
+ if growth_rate > 1/self.ring.tau[2]:
+ return True
+ else:
+ return False
+ else:
+ return growth_rate, mu
+
def is_DC_Robinson_stable(self, I0):
"""
Check DC Robinson stability - Eq. (6.1.1) [1]
@@ -1212,16 +1284,23 @@ class ProportionalIntegralLoop():
In case of super conducting cavity (QL > 1e6), the Pgain of ~100
can be used.
In a "bad" parameter set, unstable oscillation of the cavity voltage
- can be caused. So, a parameter scan of the gain should be made.
+ can be caused. So, a parameter scan of the gain should be made.
+ Igain * ring.T1 / dtau is Ki defined as a the coefficients for
+ integral part in of [1], where dtau is a clock period of the PI controller.
sample_num : int
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
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
@@ -1244,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
-----
@@ -1260,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
--------
@@ -1270,6 +1355,8 @@ class ProportionalIntegralLoop():
QL=11800, fs0=23kHz
==> gain=[0.5,1e4], sample_num=8, every=7(13ns), delay=500(1us)
is one reasonable parameter set.
+ The practical clock period is 13ns.
+ ==> Igain_PF = Igain_mbtrack * Trf / 13ns = Igain_mbtrack * 0.153
References
----------
@@ -1279,10 +1366,19 @@ class ProportionalIntegralLoop():
of synchrotron light sources. PRAB, 21(1), 012001.
"""
- def __init__(self, ring, cav_res, gain, sample_num, every, delay, 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)
+ self.ig_modulation_signal = np.zeros(self.ring.h, dtype=complex)
self.gain = gain
self.FF = FF
@@ -1310,6 +1406,7 @@ class ProportionalIntegralLoop():
self.sample_list = range(0, self.ring.h, self.every)
+ self.IIR_init(IIR_cutoff)
self.init_FFconst()
# Pre caclulation for Ig2Vg
@@ -1340,13 +1437,14 @@ class ProportionalIntegralLoop():
# 1) recording diff as a first item of the list
mean_vc = np.mean(vc_list[index:self.sample_num + index]) * np.exp(
-1j * self.cav_res.theta)
- self.diff_record[0] = self.cav_res.Vc - mean_vc
+ self.diff_record[0] = self.cav_res.Vc - self.IIR(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 = self.ig_phasor + self.ig_modulation_signal
self.ig_phasor_record = self.ig_phasor
if apply_changes:
@@ -1407,6 +1505,40 @@ class ProportionalIntegralLoop():
"""
return Vg * (1 - 1j * np.tan(self.cav_res.psi)) / self.cav_res.RL
+ def IIR_init(self, cutoff):
+ """
+ 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
+ else:
+ omega = 2.0 * np.pi * cutoff
+ T = self.ring.T1 * self.every
+ alpha = np.cos(omega*T)-1
+ tmp = alpha * alpha - 2 * alpha
+ 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):
+ """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 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))
class DirectFeedback(ProportionalIntegralLoop):
"""