diff --git a/mbtrack2/tracking/monitors/monitors.py b/mbtrack2/tracking/monitors/monitors.py
index 5d55db84cc2681d26e246f02ae3a687766559eba..f8641fe9985aabeaa56348c3e276e0db91bfd8b5 100644
--- a/mbtrack2/tracking/monitors/monitors.py
+++ b/mbtrack2/tracking/monitors/monitors.py
@@ -1407,7 +1407,9 @@ class CavityMonitor(Monitor):
"Pg":(buffer_size,),
"Rs":(buffer_size,),
"Q":(buffer_size,),
- "QL":(buffer_size,)}
+ "QL":(buffer_size,),
+ "Vc":(buffer_size,),
+ "theta":(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,),
@@ -1419,7 +1421,9 @@ class CavityMonitor(Monitor):
"Pg":(total_size,),
"Rs":(total_size,),
"Q":(total_size,),
- "QL":(total_size,)}
+ "QL":(total_size,),
+ "Vc":(total_size,),
+ "theta":(total_size,)}
dict_dtype = {"cavity_phasor_record":complex,
"beam_phasor_record":complex,
"generator_phasor_record":complex,
@@ -1431,7 +1435,9 @@ class CavityMonitor(Monitor):
"Pg":float,
"Rs":float,
"Q":float,
- "QL":float}
+ "QL":float,
+ "Vc":float,
+ "theta":float}
self.monitor_init(group_name, save_every, buffer_size, total_size,
dict_buffer, dict_file, file_name, mpi_mode,
diff --git a/mbtrack2/tracking/monitors/plotting.py b/mbtrack2/tracking/monitors/plotting.py
index 8afe5783d8cb976ddbeacf9f5fbd4e071e6e15fd..8adbae41cd7dfe91631637609645f043cf84df0b 100644
--- a/mbtrack2/tracking/monitors/plotting.py
+++ b/mbtrack2/tracking/monitors/plotting.py
@@ -1076,7 +1076,8 @@ def streak_beamspectrum(filename, dim="tau", f0=None, log_scale=True, fmin=None,
return fig
def plot_cavitydata(filename, cavity_name, phasor="cavity",
- plot_type="bunch", bunch_number=0, turn=None, cm_lim=None):
+ plot_type="mean", bunch_number=0, turn=None, cm_lim=None,
+ show_objective=False):
"""
Plot data recorded by CavityMonitor.
@@ -1091,6 +1092,8 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
"ig".
plot_type : str, optional
Type of plot:
+ - "mean" plots the mean phasor voltage and angle versus time for
+ non empty bunches.
- "bunch" plots the phasor voltage and angle versus time for a
given bunch.
- "turn" plots the phasor voltage and ange versus bunch index for
@@ -1108,6 +1111,10 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
Turn to plot. The default is None.
cm_lim : list [vmin, vmax], optional
Colormap scale for the "streak" plots.
+ show_objective : bool, optional
+ If True, shows the Cavity voltage and phase objective value (Vc and
+ theta) in plot type "mean" and "bunch" for "cavity" phasor.
+ Default is False.
Returns
-------
@@ -1125,13 +1132,23 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
units = [" voltage [MV]", " voltage [MV]", " voltage [MV]", " current [A]"]
units_val = [1e-6, 1e-6, 1e-6, 1]
- if plot_type == "bunch":
+ if plot_type == "bunch" or "mean":
- data = [cavity_data["cavity_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,:]]
-
+ if plot_type == "bunch":
+ data = [cavity_data["cavity_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,:]]
+ elif plot_type == "mean":
+ try:
+ bunch_index = (file["Beam"]["current"][:,0] != 0).nonzero()[0]
+ except:
+ print("Beam monitor is needed to show mean voltage.")
+ data = [np.mean(cavity_data["cavity_phasor_record"][bunch_index,:],0),
+ np.mean(cavity_data["beam_phasor_record"][bunch_index,:],0),
+ np.mean(cavity_data["generator_phasor_record"][bunch_index,:],0),
+ np.mean(cavity_data["ig_phasor_record"][bunch_index,:],0)]
+
ylabel1 = labels[ph[phasor]] + units[ph[phasor]]
ylabel2 = labels[ph[phasor]] + " phase [rad]"
@@ -1139,11 +1156,17 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
twin = ax.twinx()
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)
+ if show_objective:
+ o1, = ax.plot(time, np.array(cavity_data["Vc"])*1e-6, "r--", label="Cavity voltage objective value [MV]")
+ o2, = twin.plot(time, np.array(cavity_data["theta"]), "b--", label="Cavity phase objective value [rad]")
ax.set_xlabel("Turn number")
ax.set_ylabel(ylabel1)
twin.set_ylabel(ylabel2)
- plots = [p1, p2]
+ if show_objective:
+ plots = [p1, o1, p2, o2]
+ else:
+ plots = [p1, p2]
ax.legend(handles=plots, loc="best")
ax.yaxis.label.set_color("r")
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index dbf93a6c1ec77b283839062c92b40a81dc6607ca..7afa2dcb80f39a69d7142c986b6e7d64596e72de 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -79,9 +79,9 @@ class CavityResonator():
Ncav : int, optional
Number of cavities.
Vc : float, optinal
- Total cavity voltage in [V].
+ Total cavity voltage objective value in [V].
theta : float, optional
- Total cavity phase in [rad].
+ Total cavity phase objective value in [rad].
n_bin : int, optional
Number of bins used for the beam loading computation.
Only used if MPI is not used, otherwise n_bin must be specified in the
@@ -183,6 +183,8 @@ class CavityResonator():
Return the second derivative of total RF voltage.
deltaVRF(z, I0)
Return the generator voltage minus beam loading voltage.
+ update_feedback()
+ Force feedback update from current CavityResonator parameters.
References
----------
@@ -197,6 +199,7 @@ class CavityResonator():
def __init__(self, ring, m, Rs, Q, QL, detune, Ncav=1, Vc=0, theta=0,
n_bin=75):
self.ring = ring
+ self.feedback = []
self.m = m
self.Ncav = Ncav
if Ncav != 1:
@@ -218,7 +221,6 @@ class CavityResonator():
self.theta_gr = 0
self.Pg = 0
self.n_bin = int(n_bin)
- self.feedback = []
def init_tracking(self, beam):
"""
@@ -623,6 +625,7 @@ class CavityResonator():
def QL(self, value):
self._QL = value
self._beta = self.Q/self.QL - 1
+ self.update_feedback()
@property
def beta(self):
@@ -645,6 +648,7 @@ class CavityResonator():
self._wr = self.fr*2*np.pi
self._psi = np.arctan(self.QL*(self.fr/(self.m*self.ring.f1) -
(self.m*self.ring.f1)/self.fr))
+ self.update_feedback()
@property
def fr(self):
@@ -913,6 +917,12 @@ 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)
+
+ def update_feedback(self):
+ """Force feedback update from current CavityResonator parameters."""
+ for FB in self.feedback:
+ if isinstance(FB, (ProportionalIntegralLoop, DirectFeedback)):
+ FB.init_Ig2Vg_matrix()
class ProportionalLoop():
"""
@@ -1030,21 +1040,24 @@ class ProportionalIntegralLoop():
Proportional Integral (PI) loop to control a CavityResonator amplitude and
phase via generator current (ig) [1,2].
+ Feedback reference targets (setpoints) are cav_res.Vc and cav_res.theta.
+
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)
+
+ When the reference target is changed, it might be needed to re-initialize
+ the feedforward constant by calling init_FFconst().
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
+ gain : float list
+ Proportional gain (Pgain) and integral gain (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.
@@ -1062,13 +1075,10 @@ class ProportionalIntegralLoop():
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
+ Boolean switch to use feedforward constant.
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.
+ In case of small Pgain (QL ~ 1e4), Vc drop may cause beam loss due to static Robinson.
Default is True.
Methods
@@ -1077,6 +1087,8 @@ class ProportionalIntegralLoop():
Tracking method for the Cavity PI control feedback.
init_Ig2Vg_matrix()
Initialize matrix for Ig2Vg_matrix.
+ init_FFconst()
+ Initialize feedforward constant.
Ig2Vg_matrix()
Return Vg from Ig using matrix formalism.
Ig2Vg()
@@ -1118,15 +1130,12 @@ class ProportionalIntegralLoop():
of synchrotron light sources. PRAB, 21(1), 012001.
"""
- def __init__(self, ring, cav_res, gain, sample_num, every, delay,
- Vref=None, FF=True):
+ def __init__(self, ring, cav_res, gain, sample_num, every, delay, FF=True):
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
if delay > 0:
self.delay = int(delay)
@@ -1148,12 +1157,9 @@ class ProportionalIntegralLoop():
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)
+
+ self.init_FFconst()
# Pre caclulation for Ig2Vg
self.init_Ig2Vg_matrix()
@@ -1179,8 +1185,8 @@ class ProportionalIntegralLoop():
# 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
+ 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
# 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
@@ -1190,30 +1196,6 @@ class ProportionalIntegralLoop():
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)
- self.init_Ig2Vg_matrix()
def init_Ig2Vg_matrix(self):
"""
@@ -1228,6 +1210,13 @@ class ProportionalIntegralLoop():
for idx in np.arange(self.ring.h):
self.Ig2Vg_mat[idx:,idx] = tempV[:self.ring.h-idx]
+ def init_FFconst(self):
+ """Initialize feedforward constant."""
+ if self.FF:
+ self.FFconst = np.mean(self.ig_phasor)
+ else:
+ self.FFconst = 0
+
def Ig2Vg_matrix(self):
"""
Return Vg from Ig using matrix formalism.
@@ -1255,9 +1244,7 @@ class ProportionalIntegralLoop():
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
+ return Vg * ( 1 - 1j*np.tan(self.cav_res.psi) ) / self.cav_res.RL
class DirectFeedback(ProportionalIntegralLoop):
"""
@@ -1271,6 +1258,9 @@ class DirectFeedback(ProportionalIntegralLoop):
the CavityResonator generator parameters should be set before
initialization.
+ The loop must be added to the CavityResonator object using:
+ cav_res.feedback.append(loop)
+
Parameters
----------
DFB_gain : float