diff --git a/mbtrack2/tracking/__init__.py b/mbtrack2/tracking/__init__.py
index 38a0f8d4fc6c3e4dceae81ae5a760f3829667432..6f72151ebe7cf8ce7a8e72a12395385fc6d08fc1 100644
--- a/mbtrack2/tracking/__init__.py
+++ b/mbtrack2/tracking/__init__.py
@@ -6,7 +6,11 @@ from mbtrack2.tracking.particles import (Electron,
Particle)
from mbtrack2.tracking.synchrotron import Synchrotron
from mbtrack2.tracking.rf import (RFCavity,
- CavityResonator)
+ CavityResonator,
+ ProportionalLoop,
+ TunerLoop,
+ ProportionalIntegralLoop,
+ DirectFeedback)
from mbtrack2.tracking.parallel import Mpi
from mbtrack2.tracking.element import (Element,
LongitudinalMap,
diff --git a/mbtrack2/tracking/monitors/monitors.py b/mbtrack2/tracking/monitors/monitors.py
index dc6aecf89aa344e47f953ee201ac1101d05bd874..5d55db84cc2681d26e246f02ae3a687766559eba 100644
--- a/mbtrack2/tracking/monitors/monitors.py
+++ b/mbtrack2/tracking/monitors/monitors.py
@@ -1398,6 +1398,8 @@ class CavityMonitor(Monitor):
group_name = cavity_name
dict_buffer = {"cavity_phasor_record":(ring.h, buffer_size,),
"beam_phasor_record":(ring.h, buffer_size,),
+ "generator_phasor_record":(ring.h, buffer_size,),
+ "ig_phasor_record":(ring.h, buffer_size,),
"detune":(buffer_size,),
"psi":(buffer_size,),
"Vg":(buffer_size,),
@@ -1408,6 +1410,8 @@ class CavityMonitor(Monitor):
"QL":(buffer_size,)}
dict_file = {"cavity_phasor_record":(ring.h, total_size,),
"beam_phasor_record":(ring.h, total_size,),
+ "generator_phasor_record":(ring.h, total_size,),
+ "ig_phasor_record":(ring.h, total_size,),
"detune":(total_size,),
"psi":(total_size,),
"Vg":(total_size,),
@@ -1418,6 +1422,8 @@ class CavityMonitor(Monitor):
"QL":(total_size,)}
dict_dtype = {"cavity_phasor_record":complex,
"beam_phasor_record":complex,
+ "generator_phasor_record":complex,
+ "ig_phasor_record":complex,
"detune":float,
"psi":float,
"Vg":float,
diff --git a/mbtrack2/tracking/monitors/plotting.py b/mbtrack2/tracking/monitors/plotting.py
index 6bc498ea640a339d5944e7b3cc93251aa742b0fa..8afe5783d8cb976ddbeacf9f5fbd4e071e6e15fd 100644
--- a/mbtrack2/tracking/monitors/plotting.py
+++ b/mbtrack2/tracking/monitors/plotting.py
@@ -1087,15 +1087,16 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
cavity_name : str
Name of the CavityResonator object.
phasor : str, optional
- Type of the phasor to plot. Can be "beam" or "cavity".
+ Type of the phasor to plot. Can be "beam", "cavity", "generator" or
+ "ig".
plot_type : str, optional
Type of plot:
- "bunch" plots the phasor voltage and angle versus time for a
given bunch.
- "turn" plots the phasor voltage and ange versus bunch index for
a given turn.
- - "streak_volt" plots the phasor voltage versus bunch index and
- time.
+ - "streak_amplitude" plots the phasor amplitude versus bunch index
+ and time.
- "streak_angle" plots the phasor angle versus bunch index and
time.
- "detune" or "psi" plots the detuning or tuning angle versus time.
@@ -1114,26 +1115,29 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
Figure object with the plot on it.
"""
-
file = hp.File(filename, "r")
cavity_data = file[cavity_name]
time = np.array(cavity_data["time"])
- ph = {"cavity":0, "beam":1}
- labels = ["Cavity", "Beam"]
+ ph = {"cavity":0, "beam":1, "generator":2, "ig":3}
+ labels = ["Cavity", "Beam", "Generator", "Generator"]
+ units = [" voltage [MV]", " voltage [MV]", " voltage [MV]", " current [A]"]
+ units_val = [1e-6, 1e-6, 1e-6, 1]
if plot_type == "bunch":
data = [cavity_data["cavity_phasor_record"][bunch_number,:],
- cavity_data["beam_phasor_record"][bunch_number,:]]
+ cavity_data["beam_phasor_record"][bunch_number,:],
+ cavity_data["generator_phasor_record"][bunch_number,:],
+ cavity_data["ig_phasor_record"][bunch_number,:]]
- ylabel1 = labels[ph[phasor]] + " voltage [MV]"
+ ylabel1 = labels[ph[phasor]] + units[ph[phasor]]
ylabel2 = labels[ph[phasor]] + " phase [rad]"
fig, ax = plt.subplots()
twin = ax.twinx()
- p1, = ax.plot(time, np.abs(data[ph[phasor]])*1e-6, color="r",label=ylabel1)
+ p1, = ax.plot(time, np.abs(data[ph[phasor]])*units_val[ph[phasor]], color="r",label=ylabel1)
p2, = twin.plot(time, np.angle(data[ph[phasor]]), color="b", label=ylabel2)
ax.set_xlabel("Turn number")
ax.set_ylabel(ylabel1)
@@ -1150,20 +1154,20 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
index = np.array(time) == turn
if (index.size == 0):
raise ValueError("Turn is not valid.")
- ph = {"cavity":0, "beam":1}
- data = [np.array(cavity_data["cavity_phasor_record"])[:,index],
- np.array(cavity_data["beam_phasor_record"])[:,index]]
- labels = ["Cavity", "Beam"]
+ data = [cavity_data["cavity_phasor_record"][:,index],
+ cavity_data["beam_phasor_record"][:,index],
+ cavity_data["generator_phasor_record"][:,index],
+ cavity_data["ig_phasor_record"][:,index]]
h=len(data[0])
x=np.arange(h)
- ylabel1 = labels[ph[phasor]] + " voltage [MV]"
+ ylabel1 = labels[ph[phasor]] + units[ph[phasor]]
ylabel2 = labels[ph[phasor]] + " phase [rad]"
fig, ax = plt.subplots()
twin = ax.twinx()
- p1, = ax.plot(x, np.abs(data[ph[phasor]])*1e-6, color="r",label=ylabel1)
+ p1, = ax.plot(x, np.abs(data[ph[phasor]])*units_val[ph[phasor]], color="r",label=ylabel1)
p2, = twin.plot(x, np.angle(data[ph[phasor]]), color="b", label=ylabel2)
ax.set_xlabel("Bunch index")
ax.set_ylabel(ylabel1)
@@ -1175,14 +1179,19 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
ax.yaxis.label.set_color("r")
twin.yaxis.label.set_color("b")
- if plot_type == "streak_volt" or plot_type == "streak_phase":
+ if plot_type == "streak_amplitude" or plot_type == "streak_phase":
+
+ data = [cavity_data["cavity_phasor_record"][:,:],
+ cavity_data["beam_phasor_record"][:,:],
+ cavity_data["generator_phasor_record"][:,:],
+ cavity_data["ig_phasor_record"][:,:]]
- if plot_type == "streak_volt":
- data = np.transpose(np.abs(cavity_data["cavity_phasor_record"][:,:])*1e-6)
- ylabel = labels[ph[phasor]] + " voltage [MV]"
+ if plot_type == "streak_amplitude":
+ data = np.transpose(np.abs(data[ph[phasor]])*units_val[ph[phasor]])
+ ylabel = labels[ph[phasor]] + units[ph[phasor]]
cmap = mpl.cm.coolwarm # diverging
elif plot_type == "streak_phase":
- data = np.transpose(np.angle(cavity_data["cavity_phasor_record"][:,:]))
+ data = np.transpose(np.angle(data[ph[phasor]]))
ylabel = labels[ph[phasor]] + " phase [rad]"
cmap = mpl.cm.coolwarm # diverging
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 69a1c20f28f74503f9a11ab24ccc68a187418086..23c2de6db6e1089b61138a192fdd7a771b9f8938 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -58,6 +58,9 @@ class CavityResonator():
If used with mpi, beam.mpi.share_distributions must be called before the
track method call.
+ Different kind of RF feeback and loops can be added using:
+ cavity_resonator.feedback.append(loop)
+
Parameters
----------
ring : Synchrotron object
@@ -93,10 +96,15 @@ class CavityResonator():
Last beam phasor value of each bunch in [V].
generator_phasor : complex
Generator phasor in [V].
+ generator_phasor_record : array of complex
+ Last generator phasor value of each bunch in [V].
cavity_phasor : complex
Cavity phasor in [V].
cavity_phasor_record : array of complex
Last cavity phasor value of each bunch in [V].
+ ig_phasor_record : array of complex
+ Last current generator phasor of each bunch in [A].
+ Only used for some feedback types.
cavity_voltage : float
Cavity total voltage in [V].
cavity_phase : float
@@ -182,7 +190,7 @@ class CavityResonator():
factories. In Frontiers of Particle Beams: Factories with e+ e-Rings
(pp. 293-311). Springer, Berlin, Heidelberg.
- [2] Yamamoto, Naoto, Alexis Gamelin, and Ryutaro Nagaoka. "Investigation
+ [2] Naoto Yamamoto, Alexis Gamelin, and Ryutaro Nagaoka. "Investigation
of Longitudinal Beam Dynamics With Harmonic Cavities by Using the Code
Mbtrack." IPAC’19, Melbourne, Australia, 2019.
"""
@@ -202,6 +210,7 @@ class CavityResonator():
self.theta = theta
self.beam_phasor = np.zeros(1, dtype=complex)
self.beam_phasor_record = np.zeros((self.ring.h), dtype=complex)
+ self.generator_phasor_record = np.zeros((self.ring.h), dtype=complex)
self.tracking = False
self.Vg = 0
self.theta_g = 0
@@ -209,6 +218,7 @@ class CavityResonator():
self.theta_gr = 0
self.Pg = 0
self.n_bin = int(n_bin)
+ self.feedback = []
def init_tracking(self, beam):
"""
@@ -295,7 +305,7 @@ class CavityResonator():
ind = (sorted_index == i)
phase = self.m * self.ring.omega1 * (center0 + self.ring.T1* (index + self.ring.h * self.nturn))
- Vgene = self.Vg*np.cos(phase + self.theta_g)
+ Vgene = np.real(self.generator_phasor_record[index]*np.exp(1j*phase))
Vbeam = np.real(self.beam_phasor)
Vtot = Vgene + Vbeam - charge_per_mp*self.loss_factor*mp_per_bin
energy_change[ind] = Vtot / self.ring.E0
@@ -315,10 +325,13 @@ class CavityResonator():
# phasor decay to be at t=0 of the next bunch
self.phasor_decay(self.ring.T1, ref_frame="beam")
+
+ # apply different kind of RF feedback
+ for fb in self.feedback:
+ fb.track()
self.nturn += 1
-
def init_phasor_track(self, beam):
"""
Initialize the beam phasor for a given beam distribution using a
@@ -515,7 +528,15 @@ class CavityResonator():
@property
def cavity_phasor_record(self):
"""Last cavity phasor value of each bunch in [V]"""
- return self.generator_phasor + self.beam_phasor_record
+ return self.generator_phasor_record + self.beam_phasor_record
+
+ @property
+ def ig_phasor_record(self):
+ """Last current generator phasor of each bunch in [A]"""
+ for FB in self.feedback:
+ if isinstance(FB, (ProportionalIntegralLoop, DirectFeedback)):
+ return FB.ig_phasor_record
+ return np.zeros(self.ring.h)
@property
def cavity_voltage(self):
@@ -786,6 +807,8 @@ class CavityResonator():
self.Vg = self.Vgr*np.cos(self.psi)
# Generator phase [rad]
self.theta_g = self.theta_gr + self.psi
+ # Set generator_phasor_record
+ self.generator_phasor_record = np.ones(self.ring.h)*self.generator_phasor
def plot_phasor(self, I0):
"""
@@ -889,4 +912,607 @@ class CavityResonator():
def deltaVRF(self, z, I0, F = 1, PHI = 0):
"""Return the generator voltage minus beam loading voltage taking into account form factor amplitude F and form factor phase PHI"""
- return -1*self.Vg*(self.ring.k1*self.m)**2*np.cos(self.ring.k1*self.m*z + self.theta_g) - self.Vb(I0)*F*(self.ring.k1*self.m)**2*np.cos(self.ring.k1*self.m*z + self.psi - PHI)
\ No newline at end of file
+ return -1*self.Vg*(self.ring.k1*self.m)**2*np.cos(self.ring.k1*self.m*z + self.theta_g) - self.Vb(I0)*F*(self.ring.k1*self.m)**2*np.cos(self.ring.k1*self.m*z + self.psi - PHI)
+
+class ProportionalLoop():
+ """
+ Proportional feedback loop to control a CavityResonator amplitude and phase.
+
+ Feedback setpoints are cav_res.Vc and cav_res.theta.
+
+ The loop must be added to the CavityResonator object using:
+ cav_res.feedback.append(loop)
+
+ Parameters
+ ----------
+ ring : Synchrotron object
+ cav_res : CavityResonator
+ The CavityResonator which is to be controlled.
+ gain_A : float
+ Amplitude (voltage) gain of the feedback.
+ gain_P : float
+ Phase gain of the feedback.
+ delay : int
+ Feedback delay in unit of turns.
+
+ """
+ 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)*self.cav_res.Vc
+ self.phase_delay = np.ones(self.delay)*self.cav_res.theta
+
+ def track(self):
+ """
+ Tracking method for the amplitude and phase loop.
+
+ Returns
+ -------
+ None.
+
+ """
+ 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] = self.cav_res.cavity_voltage
+ self.phase_delay[0] = self.cav_res.cavity_phase
+
+class TunerLoop():
+ """
+ Cavity tuner loop used to control a CavityResonator tuning (psi or detune)
+ in order to keep the phase between cavity and generator current constant.
+
+ Only a proportional controller is assumed.
+
+ The loop must be added to the CavityResonator object using:
+ cav_res.feedback.append(loop)
+
+ Parameters
+ ----------
+ ring : Synchrotron object
+ cav_res : CavityResonator
+ The CavityResonator which is to be controlled.
+ gain : float
+ Proportional gain of the tuner loop.
+ If not specified, 0.01 is used.
+ avering_period:
+ Period during which the phase difference is monitored and averaged.
+ Then the feedback correction is applied every avering_period turn.
+ Unit is turn number.
+ A value longer than one synchrotron period (1/fs) is recommended.
+ If not specified, 2-synchrotron period (2/fs) is used, although it is
+ too fast compared to the actual situation.
+ offset : float
+ Tuning offset in [rad].
+
+ """
+ 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(self.cav_res.Vc)*self.ring.f1/self.ring.h
+ avering_period = 2/fs/self.ring.T0
+
+ self.Pgain = gain
+ self.offset = offset
+ self.avering_period = int(avering_period)
+ self.diff = 0
+ self.count = 0
+
+ def track(self):
+ """
+ Tracking method for the tuner loop.
+
+ Returns
+ -------
+ None.
+
+ """
+ if self.count == self.avering_period:
+ diff = self.diff/self.avering_period - self.offset
+ self.cav_res.psi -= diff * self.Pgain
+ self.count = 0
+ self.diff = 0
+ else:
+ 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