diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 15b6f7c1080a94ce157c03bf707824f19b1c0cf8..9d2d6ea6753dad834c26eefbb1a098b334e92f46 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -56,7 +56,7 @@ class RFCavity(Element):
np.cos(self.m * self.ring.omega1 * val + self.theta))
-class CavityResonator:
+class CavityResonator():
"""Cavity resonator class for active or passive RF cavity with beam
loading or HOM, based on [1,2].
@@ -172,6 +172,8 @@ class CavityResonator:
Plot phasor diagram.
is_DC_Robinson_stable(I0)
Check DC Robinson stability.
+ is_CBI_stable(I0)
+ Check Coupled-Bunch-Instability stability.
plot_DC_Robinson_stability()
Plot DC Robinson stability limit.
init_tracking(beam)
@@ -206,20 +208,17 @@ class CavityResonator:
Mbtrack." IPAC’19, Melbourne, Australia, 2019.
"""
- def __init__(
- self,
- ring,
- m,
- Rs,
- Q,
- QL,
- detune,
- Ncav=1,
- Vc=0,
- theta=0,
- phase_noise=None,
- n_bin=75,
- ):
+ 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
@@ -233,10 +232,7 @@ class CavityResonator:
self.detune = detune
self.Vc = Vc
self.theta = theta
- # =============== add: phase noise generator==================
- self.phase_noise_generator = None
- # =============== END ============================
- self.beam_phasor = np.zeros(1, dtype=complex)
+ self.beam_phasor = 0 + 0j
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
@@ -257,8 +253,7 @@ class CavityResonator:
"""
if beam.mpi_switch:
- # Number of the tracked bunch in this processor
- self.bunch_index = beam.mpi.bunch_num
+ self.bunch_index = beam.mpi.bunch_num # Number of the tracked bunch in this processor
self.distance = beam.distance_between_bunches
self.valid_bunch_index = beam.bunch_index
@@ -283,13 +278,7 @@ class CavityResonator:
if self.tracking is False:
self.init_tracking(beam)
- # =============== add：before every turn, add phase noise to generator
- if self.phase_noise_generator is not None:
- # get a noise from phase noise generator(unit：rad)，add it to
- # theta_g
- noise = self.phase_noise_generator.generate_noise(size=1)
- self.theta_g += noise
- # =============== END =============================
+
for index, bunch in enumerate(beam):
if beam.filling_pattern[index]:
@@ -300,16 +289,16 @@ class CavityResonator:
# mpi -> get shared bunch profile for current bunch
center = beam.mpi.tau_center[rank]
profile = beam.mpi.tau_profile[rank]
- bin_length = center[1] - center[0]
- charge_per_mp = beam.mpi.charge_per_mp_all[rank]
+ bin_length = float(beam.mpi.tau_bin_length[rank][0])
+ charge_per_mp = float(beam.mpi.charge_per_mp_all[rank])
if index == self.bunch_index:
sorted_index = beam.mpi.tau_sorted_index
else:
# no mpi -> get bunch profile for current bunch
- if len(bunch) != 0:
+ if not bunch.is_empty:
(bins, sorted_index, profile,
center) = bunch.binning(n_bin=self.n_bin)
- bin_length = center[1] - center[0]
+ bin_length = bins[1] - bins[0]
charge_per_mp = bunch.charge_per_mp
self.bunch_index = index
else:
@@ -324,8 +313,7 @@ class CavityResonator:
energy_change = bunch["tau"] * 0
- # remove part of beam phasor decay to be at the start of the
- # binning (=bins[0])
+ # remove part of beam phasor decay to be at the start of the binning (=bins[0])
self.phasor_decay(center[0] - bin_length/2, ref_frame="beam")
if index != self.bunch_index:
@@ -336,25 +324,23 @@ class CavityResonator:
else:
# modify beam phasor
for i, center0 in enumerate(center):
- mp_per_bin = profile[i]
+ mp_per_bin = int(profile[i])
if mp_per_bin == 0:
self.phasor_decay(bin_length, ref_frame="beam")
continue
- ind = sorted_index == i
- phase = (self.m * self.ring.omega1 *
- (center0 + self.ring.T1 *
- (index + self.ring.h * self.nturn)))
+ ind = (sorted_index == i)
+ phase = self.m * self.ring.omega1 * (
+ center0 + self.ring.T1 *
+ (index + self.ring.h * self.nturn))
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)
+ Vtot = Vgene + Vbeam - charge_per_mp * self.loss_factor * mp_per_bin
energy_change[ind] = Vtot / self.ring.E0
- self.beam_phasor -= (2 * charge_per_mp *
- self.loss_factor * mp_per_bin)
+ self.beam_phasor -= 2 * charge_per_mp * self.loss_factor * mp_per_bin
self.phasor_decay(bin_length, ref_frame="beam")
# phasor decay to be at t=0 of the current bunch (=-1*bins[-1])
@@ -402,10 +388,11 @@ class CavityResonator:
if beam.mpi_switch:
# get shared bunch profile for current bunch
+ beam.mpi.share_distributions(beam, n_bin=self.n_bin)
center = beam.mpi.tau_center[j]
profile = beam.mpi.tau_profile[j]
- bin_length = center[1] - center[0]
- charge_per_mp = beam.mpi.charge_per_mp_all[j]
+ bin_length = float(beam.mpi.tau_bin_length[j][0])
+ charge_per_mp = float(beam.mpi.charge_per_mp_all[j])
else:
if i == 0:
# get bunch profile for current bunch
@@ -426,7 +413,7 @@ class CavityResonator:
profile = self.profile_save[j, :]
center = self.center_save[j, :]
- bin_length = center[1] - center[0]
+ bin_length = bins[1] - bins[0]
charge_per_mp = bunch.charge_per_mp
self.phasor_decay(center[0] - bin_length/2, ref_frame="rf")
@@ -455,7 +442,7 @@ class CavityResonator:
if ref_frame == "beam":
delta = self.wr
elif ref_frame == "rf":
- delta = self.wr - self.m * self.ring.omega1
+ delta = (self.wr - self.m * self.ring.omega1)
self.beam_phasor = self.beam_phasor * np.exp(
(-1 / self.filling_time + 1j*delta) * time)
@@ -489,7 +476,7 @@ class CavityResonator:
if ref_frame == "beam":
delta = self.wr
elif ref_frame == "rf":
- delta = self.wr - self.m * self.ring.omega1
+ delta = (self.wr - self.m * self.ring.omega1)
n_bin = len(profile)
@@ -526,13 +513,14 @@ class CavityResonator:
if self.tracking is False:
self.init_tracking(beam)
- N = self.n_bin - 1
- delta = self.wr - self.m * self.ring.omega1
+ N = self.n_bin
+ delta = (self.wr - self.m * self.ring.omega1)
n_turn = int(self.filling_time / self.ring.T0 * 10)
T = np.ones(self.ring.h) * self.ring.T1
bin_length = np.zeros(self.ring.h)
charge_per_mp = np.zeros(self.ring.h)
+ bins = np.zeros((N + 1, self.ring.h))
profile = np.zeros((N, self.ring.h))
center = np.zeros((N, self.ring.h))
@@ -543,17 +531,17 @@ class CavityResonator:
beam.mpi.share_distributions(beam, n_bin=self.n_bin)
center[:, index] = beam.mpi.tau_center[j]
profile[:, index] = beam.mpi.tau_profile[j]
- bin_length[index] = center[1, index] - center[0, index]
- charge_per_mp[index] = beam.mpi.charge_per_mp_all[j]
+ bin_length[index] = float(beam.mpi.tau_bin_length[j][0])
+ charge_per_mp[index] = float(beam.mpi.charge_per_mp_all[j])
else:
- (bins, sorted_index, profile[:, index],
- center[:, index]) = (bunch.binning(n_bin=self.n_bin))
- bin_length[index] = center[1, index] - center[0, index]
+ (bins[:, index], sorted_index, profile[:, index],
+ center[:, index]) = bunch.binning(n_bin=self.n_bin)
+ bin_length[index] = bins[1, index] - bins[0, index]
charge_per_mp[index] = bunch.charge_per_mp
- T[index] -= center[-1, index] + bin_length[index] / 2
+ T[index] -= (center[-1, index] + bin_length[index] / 2)
if index != 0:
- T[index - 1] += center[0, index] - bin_length[index] / 2
- T[self.ring.h - 1] += center[0, 0] - bin_length[0] / 2
+ T[index - 1] += (center[0, index] - bin_length[index] / 2)
+ T[self.ring.h - 1] += (center[0, 0] - bin_length[0] / 2)
# Compute matrix coefficients
k = np.arange(0, N)
@@ -570,13 +558,11 @@ class CavityResonator:
for i in range(n_turn):
# Phasor decay during one turn
self.phasor_decay(self.ring.T0, ref_frame="rf")
- # Phasor evolution due to induced voltage by marco-particles during
- # one turn
+ # Phasor evolution due to induced voltage by marco-particles during one turn
sum_val = -2 * sum_tot * self.loss_factor
self.beam_phasor += sum_val
- # Replace phasor at t=0 (synchronous particle) of the first non empty
- # bunch.
+ # Replace phasor at t=0 (synchronous particle) of the first non empty bunch.
idx0 = self.valid_bunch_index[0]
self.phasor_decay(center[-1, idx0] + bin_length[idx0] / 2,
ref_frame="rf")
@@ -871,25 +857,23 @@ class CavityResonator:
"""
- # Generator power [W] - Eq. (4.1.2) [1] corrected with factor
- # (1+beta)**2 instead of (1+beta**2)
- self.Pg = (
- self.Vc**2 * (1 + self.beta)**2 /
- (2 * self.Rs * 4 * self.beta * np.cos(self.psi)**2) *
- ((np.cos(self.theta) + 2 * I0 * self.Rs /
- (self.Vc * (1 + self.beta)) * np.cos(self.psi)**2)**2 +
- (np.sin(self.theta) + 2 * I0 * self.Rs /
- (self.Vc *
- (1 + self.beta)) * np.cos(self.psi) * np.sin(self.psi))**2))
+ # Generator power [W] - Eq. (4.1.2) [1] corrected with factor (1+beta)**2 instead of (1+beta**2)
+ self.Pg = self.Vc**2 * (1 + self.beta)**2 / (
+ 2 * self.Rs * 4 * self.beta * np.cos(self.psi)**2) * (
+ (np.cos(self.theta) + 2 * I0 * self.Rs /
+ (self.Vc * (1 + self.beta)) * np.cos(self.psi)**2)**2 +
+ (np.sin(self.theta) + 2 * I0 * self.Rs /
+ (self.Vc *
+ (1 + self.beta)) * np.cos(self.psi) * np.sin(self.psi))**2)
# Generator voltage at resonance [V] - Eq. (3.2.2) [1]
- self.Vgr = (2 * self.beta**(1 / 2) / (1 + self.beta) *
- (2 * self.Rs * self.Pg)**(1 / 2))
+ self.Vgr = 2 * self.beta**(1 / 2) / (1 + self.beta) * (
+ 2 * self.Rs * self.Pg)**(1 / 2)
# Generator phase at resonance [rad] - from Eq. (4.1.1)
- self.theta_gr = (np.arctan(
+ self.theta_gr = np.arctan(
(self.Vc * np.sin(self.theta) +
self.Vbr(I0) * np.cos(self.psi) * np.sin(self.psi)) /
(self.Vc * np.cos(self.theta) +
- self.Vbr(I0) * np.cos(self.psi)**2)) - self.psi)
+ self.Vbr(I0) * np.cos(self.psi)**2)) - self.psi
# Generator voltage [V]
self.Vg = self.Vgr * np.cos(self.psi)
# Generator phase [rad]
@@ -916,14 +900,12 @@ class CavityResonator:
def make_legend_arrow(legend, orig_handle, xdescent, ydescent, width,
height, fontsize):
- p = mpatches.FancyArrow(
- 0,
- 0.5 * height,
- width,
- 0,
- length_includes_head=True,
- head_width=0.75 * height,
- )
+ p = mpatches.FancyArrow(0,
+ 0.5 * height,
+ width,
+ 0,
+ length_includes_head=True,
+ head_width=0.75 * height)
return p
fig = plt.figure()
@@ -937,40 +919,33 @@ class CavityResonator:
1,
alpha=0.5,
width=0.015,
- edgecolor="black",
+ edgecolor='black',
lw=2)
- arr2 = ax.arrow(
- self.psi + np.pi,
- 0,
- 0,
- self.Vb(I0) / self.Vc,
- alpha=0.5,
- width=0.015,
- edgecolor="red",
- lw=2,
- )
-
- arr3 = ax.arrow(
- self.theta_g,
- 0,
- 0,
- self.Vg / self.Vc,
- alpha=0.5,
- width=0.015,
- edgecolor="blue",
- lw=2,
- )
+ arr2 = ax.arrow(self.psi + np.pi,
+ 0,
+ 0,
+ self.Vb(I0) / self.Vc,
+ alpha=0.5,
+ width=0.015,
+ edgecolor='red',
+ lw=2)
+
+ arr3 = ax.arrow(self.theta_g,
+ 0,
+ 0,
+ self.Vg / self.Vc,
+ alpha=0.5,
+ width=0.015,
+ edgecolor='blue',
+ lw=2)
ax.set_rticks([]) # less radial ticks
- plt.legend(
- [arr1, arr2, arr3],
- ["Vc", "Vb", "Vg"],
- handler_map={
- mpatches.FancyArrow:
- HandlerPatch(patch_func=make_legend_arrow),
- },
- )
+ plt.legend([arr1, arr2, arr3], ['Vc', 'Vb', 'Vg'],
+ handler_map={
+ mpatches.FancyArrow:
+ HandlerPatch(patch_func=make_legend_arrow),
+ })
return fig
@@ -980,7 +955,7 @@ class CavityResonator:
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 is a wraper around lcbi_growth_rate to caluclate the CBI
@@ -1025,8 +1000,7 @@ class CavityResonator:
Vrf=self.Vc,
fr=self.fr,
RL=self.RL,
- QL=self.QL,
- )
+ QL=self.QL)
if modes is not None:
growth_rates = growth_rates[modes]
@@ -1034,9 +1008,9 @@ class CavityResonator:
if bool_return:
if growth_rate > 1 / self.ring.tau[2]:
- return True
- else:
return False
+ else:
+ return True
else:
return growth_rate, mu
@@ -1054,8 +1028,8 @@ class CavityResonator:
bool
"""
- return (2 * self.Vc * np.sin(self.theta) +
- self.Vbr(I0) * np.sin(2 * self.psi) > 0)
+ return 2 * self.Vc * np.sin(self.theta) + self.Vbr(I0) * np.sin(
+ 2 * self.psi) > 0
def plot_DC_Robinson_stability(self, detune_range=[-1e5, 1e5]):
"""
@@ -1159,9 +1133,8 @@ class CavityResonator:
# Goes from (-T1/2) to (T1/2 + DeltaT) in n_points steps
for i in range(n_points):
- phase = (self.m * self.ring.omega1 *
- (pos[i] + self.ring.T1 *
- (index + self.ring.h * self.nturn)))
+ phase = self.m * self.ring.omega1 * (
+ pos[i] + self.ring.T1 * (index + self.ring.h * self.nturn))
Vgene = np.real(self.generator_phasor_record[index] *
np.exp(1j * phase))
Vbeam = np.real(self.beam_phasor)
@@ -1170,10 +1143,9 @@ class CavityResonator:
self.phasor_decay(DeltaT, ref_frame="beam")
# Get back to t=0
- self.phasor_decay(
- -DeltaT * n_points + self.ring.T1 / 2 - index * self.ring.T1,
- ref_frame="beam",
- )
+ self.phasor_decay(-DeltaT * n_points + self.ring.T1 / 2 -
+ index * self.ring.T1,
+ ref_frame="beam")
return pos, voltage_rec