# -*- coding: utf-8 -*-
"""
This module handles radio-frequency (RF) cavitiy elements.
"""
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.legend_handler import HandlerPatch
from mbtrack2.instability import lcbi_growth_rate
from mbtrack2.tracking.element import Element
[docs]class RFCavity(Element):
"""
Perfect RF cavity class for main and harmonic RF cavities.
Use cosine definition.
Parameters
----------
ring : Synchrotron object
m : int
Harmonic number of the cavity
Vc : float
Amplitude of cavity voltage [V]
theta : float
Phase of Cavity voltage
"""
def __init__(self, ring, m, Vc, theta):
self.ring = ring
self.m = m
self.Vc = Vc
self.theta = theta
[docs] @Element.parallel
def track(self, bunch):
"""
Tracking method for the element.
No bunch to bunch interaction, so written for Bunch objects and
@Element.parallel is used to handle Beam objects.
Parameters
----------
bunch : Bunch or Beam object
"""
bunch["delta"] += self.Vc / self.ring.E0 * np.cos(
self.m * self.ring.omega1 * bunch["tau"] + self.theta)
[docs] def value(self, val):
return self.Vc / self.ring.E0 * np.cos(self.m * self.ring.omega1 *
val + self.theta)
[docs]class CavityResonator():
"""Cavity resonator class for active or passive RF cavity with beam
loading or HOM, based on [1,2].
Use cosine definition.
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
m : int or float
Harmonic number of the cavity.
Rs : float
Shunt impedance of the cavities in [Ohm], defined as 0.5*Vc*Vc/Pc.
If Ncav = 1, used for the total shunt impedance.
If Ncav > 1, used for the shunt impedance per cavity.
Q : float
Quality factor of the cavity.
QL : float
Loaded quality factor of the cavity.
detune : float
Detuing of the cavity in [Hz], defined as (fr - m*ring.f1).
Ncav : int, optional
Number of cavities.
Vc : float, optinal
Total cavity voltage objective value in [V].
theta : float, optional
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
beam.mpi.share_distributions method.
The default is 75.
Attributes
----------
beam_phasor : complex
Beam phasor in [V].
beam_phasor_record : array of complex
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
Cavity total phase in [rad].
loss_factor : float
Cavity loss factor in [V/C].
Rs_per_cavity : float
Shunt impedance of a single cavity in [Ohm], defined as 0.5*Vc*Vc/Pc.
beta : float
Coupling coefficient of the cavity.
fr : float
Resonance frequency of the cavity in [Hz].
wr : float
Angular resonance frequency in [Hz.rad].
psi : float
Tuning angle in [rad].
filling_time : float
Cavity filling time in [s].
Pc : float
Power dissipated in the cavity walls in [W].
Pg : float
Generator power in [W].
Vgr : float
Generator voltage at resonance in [V].
theta_gr : float
Generator phase at resonance in [rad].
Vg : float
Generator voltage in [V].
theta_g : float
Generator phase in [rad].
tracking : bool
True if the tracking has been initialized.
bunch_index : int
Number of the tracked bunch in the current core.
distance : array
Distance between bunches.
valid_bunch_index : array
Methods
-------
Vbr(I0)
Return beam voltage at resonance in [V].
Vb(I0)
Return beam voltage in [V].
Pb(I0)
Return power transmitted to the beam in [W].
Pr(I0)
Return power reflected back to the generator in [W].
Z(f)
Cavity impedance in [Ohm] for a given frequency f in [Hz].
set_optimal_coupling(I0)
Set coupling to optimal value.
set_optimal_detune(I0)
Set detuning to optimal conditions.
set_generator(I0)
Set generator parameters.
plot_phasor(I0)
Plot phasor diagram.
is_DC_Robinson_stable(I0)
Check DC Robinson stability.
plot_DC_Robinson_stability()
Plot DC Robinson stability limit.
init_tracking(beam)
Initialization of the tracking.
track(beam)
Tracking method.
phasor_decay(time)
Compute the beam phasor decay during a given time span.
phasor_evol(profile, bin_length, charge_per_mp)
Compute the beam phasor evolution during the crossing of a bunch.
VRF(z, I0)
Return the total RF voltage.
dVRF(z, I0)
Return derivative of total RF voltage.
ddVRF(z, I0)
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.
sample_voltage()
Sample the voltage seen by a zero charge particle during an RF period.
References
----------
[1] Wilson, P. B. (1994). Fundamental-mode rf design in e+ e− storage ring
factories. In Frontiers of Particle Beams: Factories with e+ e-Rings
(pp. 293-311). Springer, Berlin, Heidelberg.
[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.
"""
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:
self.Rs_per_cavity = Rs
else:
self.Rs = Rs
self.Q = Q
self.QL = QL
self.detune = detune
self.Vc = Vc
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
self.Vgr = 0
self.theta_gr = 0
self.Pg = 0
self.n_bin = int(n_bin)
[docs] def init_tracking(self, beam):
"""
Initialization of the tracking.
Parameters
----------
beam : Beam object
"""
if beam.mpi_switch:
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
self.tracking = True
self.nturn = 0
[docs] def track(self, beam):
"""
Track a Beam object through the CavityResonator object.
Can be used with or without mpi.
If used with mpi, beam.mpi.share_distributions must be called before.
The beam phasor is given at t=0 (synchronous particle) of the first
non empty bunch.
Parameters
----------
beam : Beam object
"""
if self.tracking is False:
self.init_tracking(beam)
for index, bunch in enumerate(beam):
if beam.filling_pattern[index]:
if beam.mpi_switch:
# get rank of bunch n° index
rank = beam.mpi.bunch_to_rank(index)
# 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]
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:
(bins, sorted_index, profile,
center) = bunch.binning(n_bin=self.n_bin)
bin_length = center[1] - center[0]
charge_per_mp = bunch.charge_per_mp
self.bunch_index = index
else:
# Update filling pattern
beam.update_filling_pattern()
beam.update_distance_between_bunches()
# save beam phasor value
self.beam_phasor_record[index] = self.beam_phasor
# phasor decay to be at t=0 of the next bunch
self.phasor_decay(self.ring.T1, ref_frame="beam")
continue
energy_change = bunch["tau"] * 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:
self.phasor_evol(profile,
bin_length,
charge_per_mp,
ref_frame="beam")
else:
# modify beam phasor
for i, center0 in enumerate(center):
mp_per_bin = 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))
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
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])
self.phasor_decay(-1 * (center[-1] + bin_length/2),
ref_frame="beam")
if index == self.bunch_index:
# apply kick
bunch["delta"] += energy_change
# save beam phasor value
self.beam_phasor_record[index] = self.beam_phasor
# 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
[docs] def init_phasor_track(self, beam):
"""
Initialize the beam phasor for a given beam distribution using a
tracking like method.
Follow the same steps as the track method but in the "rf" reference
frame and without any modifications on the beam.
Parameters
----------
beam : Beam object
"""
if self.tracking is False:
self.init_tracking(beam)
n_turn = int(self.filling_time / self.ring.T0 * 10)
for i in range(n_turn):
for j, bunch in enumerate(beam.not_empty):
index = self.valid_bunch_index[j]
if beam.mpi_switch:
# get shared bunch profile for current bunch
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]
else:
if i == 0:
# get bunch profile for current bunch
(bins, sorted_index, profile,
center) = bunch.binning(n_bin=self.n_bin)
if j == 0:
self.profile_save = np.zeros((
len(beam),
len(profile),
))
self.center_save = np.zeros((
len(beam),
len(center),
))
self.profile_save[j, :] = profile
self.center_save[j, :] = center
else:
profile = self.profile_save[j, :]
center = self.center_save[j, :]
bin_length = center[1] - center[0]
charge_per_mp = bunch.charge_per_mp
self.phasor_decay(center[0] - bin_length/2, ref_frame="rf")
self.phasor_evol(profile,
bin_length,
charge_per_mp,
ref_frame="rf")
self.phasor_decay(-1 * (center[-1] + bin_length/2),
ref_frame="rf")
self.phasor_decay((self.distance[index] * self.ring.T1),
ref_frame="rf")
[docs] def phasor_decay(self, time, ref_frame="beam"):
"""
Compute the beam phasor decay during a given time span, assuming that
no particles are crossing the cavity during the time span.
Parameters
----------
time : float
Time span in [s], can be positive or negative.
ref_frame : string, optional
Reference frame to be used, can be "beam" or "rf".
"""
if ref_frame == "beam":
delta = self.wr
elif ref_frame == "rf":
delta = (self.wr - self.m * self.ring.omega1)
self.beam_phasor = self.beam_phasor * np.exp(
(-1 / self.filling_time + 1j*delta) * time)
[docs] def phasor_evol(self,
profile,
bin_length,
charge_per_mp,
ref_frame="beam"):
"""
Compute the beam phasor evolution during the crossing of a bunch using
an analytic formula [1].
Assume that the phasor decay happens before the beam loading.
Parameters
----------
profile : array
Longitudinal profile of the bunch in [number of macro-particle].
bin_length : float
Length of a bin in [s].
charge_per_mp : float
Charge per macro-particle in [C].
ref_frame : string, optional
Reference frame to be used, can be "beam" or "rf".
References
----------
[1] mbtrack2 manual.
"""
if ref_frame == "beam":
delta = self.wr
elif ref_frame == "rf":
delta = (self.wr - self.m * self.ring.omega1)
n_bin = len(profile)
# Phasor decay during crossing time
deltaT = n_bin * bin_length
self.phasor_decay(deltaT, ref_frame)
# Phasor evolution due to induced voltage by marco-particles
k = np.arange(0, n_bin)
var = np.exp(
(-1 / self.filling_time + 1j*delta) * (n_bin-k) * bin_length)
sum_tot = np.sum(profile * var)
sum_val = -2 * sum_tot * charge_per_mp * self.loss_factor
self.beam_phasor += sum_val
[docs] def init_phasor(self, beam):
"""
Initialize the beam phasor for a given beam distribution using an
analytic formula [1].
No modifications on the Beam object.
Parameters
----------
beam : Beam object
References
----------
[1] mbtrack2 manual.
"""
# Initialization
if self.tracking is False:
self.init_tracking(beam)
N = self.n_bin - 1
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)
profile = np.zeros((N, self.ring.h))
center = np.zeros((N, self.ring.h))
# Gather beam distribution data
for j, bunch in enumerate(beam.not_empty):
index = self.valid_bunch_index[j]
if beam.mpi_switch:
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]
else:
(bins, sorted_index, profile[:, index],
center[:, index]) = bunch.binning(n_bin=self.n_bin)
bin_length[index] = center[1, index] - center[0, index]
charge_per_mp[index] = bunch.charge_per_mp
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)
# Compute matrix coefficients
k = np.arange(0, N)
Tkj = np.zeros((N, self.ring.h))
for j in range(self.ring.h):
sum_t = np.array(
[T[n] + N * bin_length[n] for n in range(j + 1, self.ring.h)])
Tkj[:, j] = (N-k) * bin_length[j] + T[j] + np.sum(sum_t)
var = np.exp((-1 / self.filling_time + 1j*delta) * Tkj)
sum_tot = np.sum((profile*charge_per_mp) * var)
# Use the formula n_turn times
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
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.
idx0 = self.valid_bunch_index[0]
self.phasor_decay(center[-1, idx0] + bin_length[idx0] / 2,
ref_frame="rf")
@property
def generator_phasor(self):
"""Generator phasor in [V]"""
return self.Vg * np.exp(1j * self.theta_g)
@property
def cavity_phasor(self):
"""Cavity total phasor in [V]"""
return self.generator_phasor + self.beam_phasor
@property
def cavity_phasor_record(self):
"""Last cavity phasor value of each bunch in [V]"""
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 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):
"""Cavity total voltage in [V]"""
return np.abs(self.cavity_phasor)
@property
def cavity_phase(self):
"""Cavity total phase in [rad]"""
return np.angle(self.cavity_phasor)
@property
def beam_voltage(self):
"""Beam loading voltage in [V]"""
return np.abs(self.beam_phasor)
@property
def beam_phase(self):
"""Beam loading phase in [rad]"""
return np.angle(self.beam_phasor)
@property
def loss_factor(self):
"""Cavity loss factor in [V/C]"""
return self.wr * self.Rs / (2 * self.Q)
@property
def m(self):
"""Harmonic number of the cavity"""
return self._m
@m.setter
def m(self, value):
self._m = value
@property
def Ncav(self):
"""Number of cavities"""
return self._Ncav
@Ncav.setter
def Ncav(self, value):
self._Ncav = value
@property
def Rs_per_cavity(self):
"""Shunt impedance of a single cavity in [Ohm], defined as
0.5*Vc*Vc/Pc."""
return self._Rs_per_cavity
@Rs_per_cavity.setter
def Rs_per_cavity(self, value):
self._Rs_per_cavity = value
@property
def Rs(self):
"""Shunt impedance [ohm]"""
return self.Rs_per_cavity * self.Ncav
@Rs.setter
def Rs(self, value):
self.Rs_per_cavity = value / self.Ncav
@property
def RL(self):
"""Loaded shunt impedance [ohm]"""
return self.Rs / (1 + self.beta)
@property
def Q(self):
"""Quality factor"""
return self._Q
@Q.setter
def Q(self, value):
self._Q = value
@property
def QL(self):
"""Loaded quality factor"""
return self._QL
@QL.setter
def QL(self, value):
self._QL = value
self._beta = self.Q / self.QL - 1
self.update_feedback()
@property
def beta(self):
"""Coupling coefficient"""
return self._beta
@beta.setter
def beta(self, value):
self.QL = self.Q / (1+value)
@property
def detune(self):
"""Cavity detuning [Hz] - defined as (fr - m*f1)"""
return self._detune
@detune.setter
def detune(self, value):
self._detune = value
self._fr = self.detune + self.m * self.ring.f1
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):
"""Resonance frequency of the cavity in [Hz]"""
return self._fr
@fr.setter
def fr(self, value):
self.detune = value - self.m * self.ring.f1
@property
def wr(self):
"""Angular resonance frequency in [Hz.rad]"""
return self._wr
@wr.setter
def wr(self, value):
self.detune = (value - self.m * self.ring.f1) * 2 * np.pi
@property
def psi(self):
"""Tuning angle in [rad]"""
return self._psi
@psi.setter
def psi(self, value):
delta = (self.ring.f1 * self.m * np.tan(value) /
self.QL)**2 + 4 * (self.ring.f1 * self.m)**2
fr = (self.ring.f1 * self.m * np.tan(value) / self.QL +
np.sqrt(delta)) / 2
self.detune = fr - self.m * self.ring.f1
@property
def filling_time(self):
"""Cavity filling time in [s]"""
return 2 * self.QL / self.wr
@property
def Pc(self):
"""Power dissipated in the cavity walls in [W]"""
return self.Vc**2 / (2 * self.Rs)
[docs] def Pb(self, I0):
"""
Return power transmitted to the beam in [W] - near Eq. (4.2.3) in [1].
Parameters
----------
I0 : float
Beam current in [A].
Returns
-------
float
Power transmitted to the beam in [W].
"""
return I0 * self.Vc * np.cos(self.theta)
[docs] def Pr(self, I0):
"""
Power reflected back to the generator in [W].
Parameters
----------
I0 : float
Beam current in [A].
Returns
-------
float
Power reflected back to the generator in [W].
"""
return self.Pg - self.Pb(I0) - self.Pc
[docs] def Vbr(self, I0):
"""
Return beam voltage at resonance in [V].
Parameters
----------
I0 : float
Beam current in [A].
Returns
-------
float
Beam voltage at resonance in [V].
"""
return 2 * I0 * self.Rs / (1 + self.beta)
[docs] def Vb(self, I0):
"""
Return beam voltage in [V].
Parameters
----------
I0 : float
Beam current in [A].
Returns
-------
float
Beam voltage in [V].
"""
return self.Vbr(I0) * np.cos(self.psi)
[docs] def Z(self, f):
"""Cavity impedance in [Ohm] for a given frequency f in [Hz]"""
return self.RL / (1 + 1j * self.QL * (self.fr / f - f / self.fr))
[docs] def set_optimal_detune(self, I0, F=1):
"""
Set detuning to optimal conditions - second Eq. (4.2.1) in [1].
Parameters
----------
I0 : float
Beam current in [A].
"""
self.psi = np.arctan(-self.Vbr(I0) * F / self.Vc * np.sin(self.theta))
[docs] def set_optimal_coupling(self, I0):
"""
Set coupling to optimal value - Eq. (4.2.3) in [1].
Parameters
----------
I0 : float
Beam current in [A].
"""
self.beta = 1 + (2 * I0 * self.Rs * np.cos(self.theta) / self.Vc)
[docs] def set_generator(self, I0):
"""
Set generator parameters (Pg, Vgr, theta_gr, Vg and theta_g) for a
given current and set of parameters.
Parameters
----------
I0 : float
Beam current in [A].
"""
# 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)
# Generator phase at resonance [rad] - from Eq. (4.1.1)
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
# Generator voltage [V]
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
[docs] def plot_phasor(self, I0):
"""
Plot phasor diagram showing the vector addition of generator and beam
loading voltage.
Parameters
----------
I0 : float
Beam current in [A].
Returns
-------
Figure.
"""
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)
return p
fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
ax.set_rmax(
max([1.2,
self.Vb(I0) / self.Vc * 1.2, self.Vg / self.Vc * 1.2]))
arr1 = ax.arrow(self.theta,
0,
0,
1,
alpha=0.5,
width=0.015,
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)
ax.set_rticks([]) # less radial ticks
plt.legend([arr1, arr2, arr3], ['Vc', 'Vb', 'Vg'],
handler_map={
mpatches.FancyArrow:
HandlerPatch(patch_func=make_legend_arrow),
})
return fig
[docs] 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
[docs] def is_DC_Robinson_stable(self, I0):
"""
Check DC Robinson stability - Eq. (6.1.1) [1]
Parameters
----------
I0 : float
Beam current in [A].
Returns
-------
bool
"""
return 2 * self.Vc * np.sin(self.theta) + self.Vbr(I0) * np.sin(
2 * self.psi) > 0
[docs] def plot_DC_Robinson_stability(self, detune_range=[-1e5, 1e5]):
"""
Plot DC Robinson stability limit.
Parameters
----------
detune_range : list or array, optional
Range of tuning to plot in [Hz].
Returns
-------
Figure.
"""
old_detune = self.psi
x = np.linspace(detune_range[0], detune_range[1], 1000)
y = []
for i in range(0, x.size):
self.detune = x[i]
y.append(-self.Vc * (1 + self.beta) /
(self.Rs * np.sin(2 * self.psi)) *
np.sin(self.theta)) # droite de stabilité
fig = plt.figure()
ax = plt.gca()
ax.plot(x, y)
ax.set_xlabel("Detune [Hz]")
ax.set_ylabel("Threshold current [A]")
ax.set_title("DC Robinson stability limit")
self.psi = old_detune
return fig
[docs] def VRF(self, z, I0, F=1, PHI=0):
"""Total RF voltage taking into account form factor amplitude F and form factor phase PHI"""
return self.Vg * np.cos(self.ring.k1 * self.m * z +
self.theta_g) - self.Vb(I0) * F * np.cos(
self.ring.k1 * self.m * z + self.psi - PHI)
[docs] def dVRF(self, z, I0, F=1, PHI=0):
"""Return derivative of total RF voltage taking into account form factor amplitude F and form factor phase PHI"""
return -1 * self.Vg * self.ring.k1 * self.m * np.sin(
self.ring.k1 * self.m * z +
self.theta_g) + self.Vb(I0) * F * self.ring.k1 * self.m * np.sin(
self.ring.k1 * self.m * z + self.psi - PHI)
[docs] def ddVRF(self, z, I0, F=1, PHI=0):
"""Return the second derivative of total RF 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)
[docs] 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)
[docs] 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()
[docs] def sample_voltage(self, n_points=1e4, index=0):
"""
Sample the voltage seen by a zero charge particle during an RF period.
Parameters
----------
n_points : int or float, optional
Number of sample points. The default is 1e4.
index : int, optional
RF bucket number to sample. Be carful if index > 0 as no new beam
loading is not taken into account here.
The default is 0.
Returns
-------
pos : array of float
Array of position from -T1/2 to T1/2.
voltage_rec : array of float
Recoring of the voltage.
"""
# Init
n_points = int(n_points)
index = 0
voltage_rec = np.zeros(n_points)
pos = np.linspace(-self.ring.T1 / 2, self.ring.T1 / 2, n_points)
DeltaT = self.ring.T1 / (n_points-1)
# From t=0 of first non empty bunch to -T1/2
self.phasor_decay(-self.ring.T1 / 2 + index * self.ring.T1,
ref_frame="beam")
# 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))
Vgene = np.real(self.generator_phasor_record[index] *
np.exp(1j * phase))
Vbeam = np.real(self.beam_phasor)
Vtot = Vgene + Vbeam
voltage_rec[i] = Vtot
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")
return pos, voltage_rec
[docs]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.
Must be supperior or equal to 1.
"""
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)
if delay < 1:
raise ValueError("delay must be >= 1.")
self.volt_delay = np.ones(self.delay) * self.cav_res.Vc
self.phase_delay = np.ones(self.delay) * self.cav_res.theta
[docs] 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
[docs]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
[docs] 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
[docs]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 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 list
Proportional gain (Pgain) and integral gain (Igain) of the feedback
system in the form of a list [Pgain, Igain].
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 the cavity voltage
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
of the tracking.
In case of small Pgain (QL ~ 1e4), cavity voltage drop may cause
beam loss due to static Robinson.
Default is True.
Methods
-------
track()
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()
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
-----
Assumption : delay >> every ~ sample_num
Adjusting ig(Vg) parameter to keep the cavity voltage constant (to target values).
The "track" method is
1) Monitoring the cavity voltage phasor
Mean cavity voltage 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 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-->(rot)-->IIR-->(-)-->(V->I,fac)-->PI-->Ig --> Vg
|
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.
The practical clock period is 13ns.
==> Igain_PF = Igain_mbtrack * Trf / 13ns = Igain_mbtrack * 0.153
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,
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
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 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
self.sample_list = range(0, self.ring.h, self.every)
self.IIR_init(IIR_cutoff)
self.init_FFconst()
# Pre caclulation for Ig2Vg
self.init_Ig2Vg_matrix()
[docs] 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]) * np.exp(
-1j * self.cav_res.theta)
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:
self.Ig2Vg()
[docs] 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]
[docs] def init_FFconst(self):
"""Initialize feedforward constant."""
if self.FF:
self.FFconst = np.mean(self.ig_phasor)
else:
self.FFconst = 0
[docs] def Ig2Vg_matrix(self):
"""
Return Vg from Ig using matrix formalism.
Warning: 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
[docs] 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.
"""
self.cav_res.generator_phasor_record = self.Ig2Vg_matrix()
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))
[docs] def Vg2Ig(self, Vg):
"""
Return Ig from Vg (assuming constant Vg).
Eq.25 of ref [2] assuming the dVg/dt = 0.
"""
return Vg * (1 - 1j * np.tan(self.cav_res.psi)) / self.cav_res.RL
[docs] 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
[docs] 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))
[docs]class DirectFeedback(ProportionalIntegralLoop):
"""
Direct RF FB (DFB) via generator current (ig) using PI controller [1,2].
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.
The loop must be added to the CavityResonator object using:
cav_res.feedback.append(loop)
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.
[2] : N. Yamamoto et al. (2023). Stability survey of a double RF system
with RF feedback loops for bunch lengthening in a low-emittance synchrotron
ring. In Proc. IPAC'23. doi:10.18429/JACoW-IPAC2023-WEPL161
"""
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 : DFB_delay or DFB_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))
[docs] 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
self.init_FFconst()
[docs] 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()
[docs] 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
[docs] 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