diff --git a/mbtrack2/tracking/__init__.py b/mbtrack2/tracking/__init__.py
index 02843e280f5cd945178b040b903c468b095765cc..fe95c6d605c281df813310b6f434a229d5bc268d 100644
--- a/mbtrack2/tracking/__init__.py
+++ b/mbtrack2/tracking/__init__.py
@@ -26,6 +26,14 @@ from mbtrack2.tracking.rf import (
RFCavity,
TunerLoop,
)
+from mbtrack2.tracking.rfdipole import (
+ DipoleCavityResonator,
+ DipoleDirectFeedback,
+ DipoleProportionalIntegralLoop,
+ DipoleProportionalLoop,
+ DipoleRFCavity,
+ DipoleTunerLoop,
+)
from mbtrack2.tracking.synchrotron import Synchrotron
from mbtrack2.tracking.wakepotential import (
LongRangeResistiveWall,
diff --git a/mbtrack2/tracking/parallel.py b/mbtrack2/tracking/parallel.py
index 58c45b116a85901e7cb173c58e0ac934805cffc7..bdd155002d72c857dcaee716c84d52b3319f2c51 100644
--- a/mbtrack2/tracking/parallel.py
+++ b/mbtrack2/tracking/parallel.py
@@ -143,7 +143,7 @@ class Mpi:
else:
return max(self.table[:, 0])
- def share_distributions(self, beam, dimensions="tau", n_bin=75):
+ def share_distributions(self, beam, dimensions="tau", dipole_plane=None,n_bin=75):
"""
Compute the bunch profiles and share it between the different bunches.
@@ -152,6 +152,9 @@ class Mpi:
beam : Beam object
dimension : str or list of str, optional
Dimensions in which the binning is done. The default is "tau".
+ dipole_plane : str or list of str, optional
+ Use for transverse high-Q resonator. The default is None.
+ "x" or "y"
n_bin : int or list of int, optional
Number of bins. The default is 75.
@@ -179,10 +182,21 @@ class Mpi:
if len(bunch) != 0:
bins, sorted_index, profile, center = bunch.binning(
dimension=dim, n_bin=n)
+ if dim == "tau" and dipole_plane is not None :
+ dipole = np.zeros((n - 1, ), dtype=np.float64)
+ for i in range(n - 1):
+ dipole[i] = bunch[dipole_plane][sorted_index == i].sum()
+ dipole = np.divide(dipole, profile, out=np.zeros_like(dipole), where=profile!=0)
+ #dipole = dipole / profile # here
+ #dipole[np.isnan(dipole)] = 0
+ #if self.bunch_num == 0:
+ # print(dipole)
else:
sorted_index = None
profile = np.zeros((n - 1, ), dtype=np.int64)
center = np.zeros((n - 1, ), dtype=np.float64)
+ if dim == "tau" and dipole_plane is not None :
+ dipole = np.zeros((n - 1, ), dtype=np.float64)
if beam.filling_pattern[self.bunch_num] is True:
beam.update_filling_pattern()
beam.update_distance_between_bunches()
@@ -201,6 +215,13 @@ class Mpi:
self.__setattr__(dim + "_sorted_index", sorted_index)
+ if dim == "tau" and dipole_plane is not None :
+ self.__setattr__("tau_dipole",
+ np.empty((self.size, n - 1), dtype=np.float64))
+ self.comm.Allgather(
+ [dipole, self.MPI.DOUBLE],
+ [self.__getattribute__("tau_dipole"), self.MPI.DOUBLE])
+
def share_means(self, beam):
"""
Compute the bunch means and share it between the different bunches.
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 3d37de9b07f80a2bb01472f405e282b334ce359c..33513c4cda904c9ccd83a91301de164e7aaa1a8e 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -334,6 +334,8 @@ class CavityResonator():
self.beam_phasor -= 2 * charge_per_mp * self.loss_factor * mp_per_bin
self.phasor_decay(bin_length, ref_frame="beam")
+ #if self.bunch_index == 0 and i < 20:
+ # print("#",self.nturn,energy_change[ind])
# phasor decay to be at t=0 of the current bunch (=-1*bins[-1])
self.phasor_decay(-1 * (center[-1] + bin_length/2),
@@ -345,6 +347,8 @@ class CavityResonator():
# save beam phasor value
self.beam_phasor_record[index] = self.beam_phasor
+ #if self.bunch_index == 0 and self.nturn==4:
+ # print("#",self.beam_phasor)
# phasor decay to be at t=0 of the next bunch
self.phasor_decay(self.ring.T1, ref_frame="beam")
@@ -1492,7 +1496,6 @@ class ProportionalIntegralLoop():
T = self.ring.T1 * self.every
alpha = np.cos(omega*T)-1
tmp = alpha * alpha - 2 * alpha
- print(tmp)
if tmp > 0:
self.IIRcoef = alpha + np.sqrt(tmp)
else:
diff --git a/mbtrack2/tracking/rfdipole.py b/mbtrack2/tracking/rfdipole.py
new file mode 100644
index 0000000000000000000000000000000000000000..adfcd424e7a8ebac85e026b9b718b5f0f4af2832
--- /dev/null
+++ b/mbtrack2/tracking/rfdipole.py
@@ -0,0 +1,1778 @@
+# -*- 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.tracking.element import Element
+
+
+class DipoleRFCavity(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/m]
+ theta : float
+ Phase of Cavity voltage
+ """
+ def __init__(self, ring, m, Vc, theta,plane="x",offset=0):
+ self.ring = ring
+ self.m = m
+ self.Vc = Vc
+ self.theta = theta
+ self.plane = plane
+ self.offset = offset
+
+ @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"] += (bunch["self.plane"] - self.offset) *self.Vc / self.ring.E0 * np.cos(
+ self.m * self.ring.omega1 * bunch["tau"] + self.theta)
+
+ def value(self, val):
+ return self.Vc / self.ring.E0 * np.cos(self.m * self.ring.omega1 *
+ val + self.theta)
+
+
+class DipoleCavityResonator():
+ """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/m].
+ 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,
+ R,
+ QL,
+ plane = "x",
+ Q=None,
+ detune = 0,
+ Ncav=1,
+ Vc=0,
+ theta=0,
+ local_beta=None,
+ n_bin=75):
+ self.ring = ring
+ self.feedback = []
+ self.m = m
+ self.Ncav = Ncav
+ if Ncav != 1:
+ self.R_per_cavity = R
+ else:
+ self.R = R
+ if Q is None:
+ self.Q = QL
+ else:
+ 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)
+ if local_beta is None:
+ self.weight = 1.0
+ else:
+ self.weight = local_beta / self.ring.optics.local_beta [self.plane]
+ if plane == "x" or plane == "xp":
+ self.plane = "x"
+ self.action_plane = "xp"
+ elif plane == "y" or plane == "yp":
+ self.plane = "y"
+ self.action_plane = "xp"
+ else:
+ raise ValueError("Plane must be: x, xp, y or yp")
+
+ 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
+
+ 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]
+ dipole = beam.mpi.tau_dipole[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
+ #if self.bunch_index == 0:
+ # print("#",center)
+ 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
+ dipole = np.zeros(self.n_bin-1, dtype=np.float64)
+ for i in range(self.n_bin - 1):
+ dipole[i] = bunch[self.plane][sorted_index == i].sum()
+ dipole = np.divide(dipole, profile, out=np.zeros_like(dipole), where=profile!=0)
+ #dipole = dipole / profile # Source of RuntimeWarning
+ #dipole[np.isnan(dipole)] = 0
+ 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)
+ continue
+
+ kick = bunch["tau"] * 0
+ # kick = bunch[self.action_plane] * 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)
+ """if beam.mpi_switch:
+ if (self.bunch_index < 2 and self.nturn < 1) :
+ print("####",self.nturn,index,self.bunch_index,(self.beam_phasor),dipole[0],center[0],flush=True)
+ else:
+ if (self.nturn < 1) :
+ print("####",self.nturn,index,self.bunch_index,(self.beam_phasor),dipole[0],center[0])"""
+
+ if index != self.bunch_index:
+ self.phasor_evol(profile,dipole, bin_length, charge_per_mp)
+ else:
+ # modify beam phasor
+ for i, dipole0 in enumerate(dipole):
+ mp_per_bin = profile[i]
+
+ if mp_per_bin == 0:
+ self.phasor_decay(bin_length)
+ 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
+ kick[ind] = np.real(self.beam_phasor) / self.ring.E0 * self.weight
+ #if self.nturn < 1:
+ # if (self.bunch_index < 2 and i < 2) or (self.bunch_index < 2 and i > self.n_bin-2) :
+ # print("#d#",self.nturn,index,i,kick[ind][0],(self.beam_phasor),flush=True)
+ self.beam_phasor -= 2j * charge_per_mp * self.loss_factor * mp_per_bin * dipole0
+ self.phasor_decay(bin_length)
+ #if (self.bunch_index < 2 and i < 2) or (self.bunch_index < 2 and i > self.n_bin-2) :
+ # print("##",self.nturn,i,(self.beam_phasor),dipole0,center[i])
+
+
+ # phasor decay to be at t=0 of the current bunch (=-1*bins[-1])
+ self.phasor_decay(-1 * (center[-1] + bin_length/2))
+
+ if index == self.bunch_index:
+ bunch[self.action_plane] += kick # apply kick
+ #if self.bunch_index == 0:
+ # print("#",self.bunch_index,index,self.beam_phasor)
+
+ # save beam phasor value
+ self.beam_phasor_record[index] = self.beam_phasor
+ #if self.bunch_index == 0 and index == 0 :
+ # print("##",self.bunch_index,index,self.beam_phasor)
+
+ # phasor decay to be at t=0 of the next bunch
+ self.phasor_decay(self.ring.T1)
+ #if self.nturn < 1 and self.bunch_index == 0 :
+ #if self.nturn <1 and index < 2:
+ # print("#",self.nturn,index,self.beam_phasor,flush=True)
+
+ # 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
+ 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]
+ dipole = beam.mpi.tau_dipole[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)
+ dipole = np.zeros(self.n_bin-1, dtype=np.float64)
+ for i in range(self.n_bin - 1):
+ dipole[i] = bunch[self.plane][sorted_index == i].sum()
+ dipole = np.divide(dipole, profile, out=np.zeros_like(dipole), where=profile!=0)
+ if j == 0:
+ self.profile_save = np.zeros((
+ len(beam),
+ len(profile),
+ ))
+ self.center_save = np.zeros((
+ len(beam),
+ len(center),
+ ))
+ self.dipole_save = np.zeros((
+ len(beam),
+ len(dipole),
+ ))
+ self.profile_save[j, :] = profile
+ self.center_save[j, :] = center
+ self.dipole_save[j, :] = dipole
+ else:
+ profile = self.profile_save[j, :]
+ center = self.center_save[j, :]
+ dipole = self.dipole_save[j, :]
+
+ bin_length = center[1] - center[0]
+ charge_per_mp = bunch.charge_per_mp
+
+ self.phasor_decay(center[0] - bin_length/2)
+ self.phasor_evol(profile,dipole,
+ bin_length,
+ charge_per_mp)
+ self.phasor_decay(-1 * (center[-1] + bin_length/2))
+ self.phasor_decay((self.distance[index] * self.ring.T1))
+
+ 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)
+
+ def phasor_evol(self, profile, dipole, bin_length, charge_per_mp):
+ """
+ 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].
+
+ References
+ ----------
+ [1] mbtrack2 manual.
+
+ """
+
+ delta = self.wr
+
+ n_bin = len(profile)
+
+ # Phasor decay during crossing time
+ deltaT = n_bin * bin_length
+ self.phasor_decay(deltaT)
+
+ # 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 * dipole * var)
+ sum_val = -2j * sum_tot * charge_per_mp * self.loss_factor
+ self.beam_phasor += sum_val
+
+ """
+ for i, dipole0 in enumerate(dipole):
+ mp_per_bin = profile[i]
+
+ if mp_per_bin == 0:
+ self.phasor_decay(bin_length)
+ continue
+
+ self.beam_phasor -= 2j * charge_per_mp * self.loss_factor * mp_per_bin * dipole0
+ self.phasor_decay(bin_length)
+ """
+
+ 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))
+ dipole = 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]
+ dipole = beam.mpi.tau_dipole[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)
+ dipole = np.zeros(self.n_bin-1, dtype=np.float64)
+ for i in range(self.n_bin - 1):
+ dipole[i] = bunch[self.plane][sorted_index == i].sum()
+ dipole = dipole / profile # Source of RuntimeWarning
+ dipole[np.isnan(dipole)] = 0
+ 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*dipole) * 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, (DipoleProportionalIntegralLoop, DipoleDirectFeedback)):
+ return FB.ig_phasor_record
+ return np.zeros(self.ring.h)
+
+ @property
+ def DFB_ig_phasor(self):
+ """Last current generator phasor of each bunch in [A]"""
+ for FB in self.feedback:
+ if isinstance(FB, DipoleDirectFeedback):
+ 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.R / (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 R_per_cavity(self):
+ """Shunt impedance of a single cavity in [Ohm], defined as
+ 0.5*Vc*Vc/Pc."""
+ return self._R_per_cavity
+
+ @R_per_cavity.setter
+ def R_per_cavity(self, value):
+ self._R_per_cavity = value
+
+ @property
+ def R(self):
+ """Shunt impedance [ohm]"""
+ return self.R_per_cavity * self.Ncav
+
+ @R.setter
+ def R(self, value):
+ self.R_per_cavity = value / self.Ncav
+
+ @property
+ def RT(self):
+ """Loaded shunt impedance [ohm]"""
+ return self.R / (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.R)
+
+ 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)
+
+ 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
+
+ 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.R / (1 + self.beta)
+
+ 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)
+
+ def Z(self, f):
+ """Cavity impedance in [Ohm] for a given frequency f in [Hz]"""
+ return self.RT * (self.fr / f) / (1 + 1j * self.QL * (self.fr / f - f / self.fr))
+
+ def set_optimal_detune(self, I0):
+ """
+ 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) / self.Vc * np.sin(self.theta))
+
+ 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.R * np.cos(self.theta) / self.Vc)
+
+ 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.R * 4 * self.beta * np.cos(self.psi)**2) * (
+ (np.cos(self.theta) + 2 * I0 * self.R /
+ (self.Vc * (1 + self.beta)) * np.cos(self.psi)**2)**2 +
+ (np.sin(self.theta) + 2 * I0 * self.R /
+ (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.R * 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
+
+ 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
+
+ def is_CBI_stable(self,I0,tune=None,mode=[0],max_freq=5,bool_return=False):
+ """
+ Check Coupled-Bunch-Instability stability
+ Effect of Direct RF feedback is not included.
+
+ This method caluclates the CBI growth rate from own impedance.
+
+ Parameters
+ ----------
+ I0 : float
+ Beam current in [A].
+ tune : float
+ fractional number of longitudinal tune
+ mode : float or list of float
+ Coupled Bunch Instability mode number
+ bool_return : bool
+ if True, return bool
+
+ Returns
+ -------
+ bool
+
+ """
+ planeIdx = 0
+ if self.plane == "y":
+ planeIdx = 1
+ if tune is None:
+ tune = self.ring.tune[planeIdx]
+ if tune > 1:
+ tune = tune - int(tune)
+ h = self.ring.h
+
+ cof = self.ring.f0*I0/2.0/self.ring.E0 * self.weight * self.ring.optics.local_beta[planeIdx]
+ if isinstance(mode,list):
+ CBImode = mode
+ else:
+ CBImode = [mode]
+
+ gr = np.zeros(len(CBImode))
+ gr_bool = np.zeros(len(CBImode),dtype=bool)
+ count = 0
+ for i in CBImode:
+ #fp = np.array([0,1,2,3,4])*self.ring.f1+i*self.ring.f0*(i+tune)
+ fp = self.ring.f0*(np.arange(max_freq)*h+(i+tune))
+ #fm = np.array([1,2,3,4,5])*self.ring.f1-i*self.ring.f0*(i+tune)
+ fm = self.ring.f0*(np.arange(1,max_freq+1)*h-(i+tune))
+ sumZ = np.sum(np.real(self.Z(fm))) - np.sum(np.real(self.Z(fp)))
+
+ gr[count] = cof*sumZ
+ gr_bool[count] = ( gr[count] > 1/self.ring.tau[planeIdx] )
+ count +=1
+
+ if bool_return:
+ return gr_bool
+ else:
+ return gr
+
+
+ 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)
+
+ 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)
+
+ 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)
+
+ 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, (DipoleProportionalIntegralLoop, DipoleDirectFeedback)):
+ FB.init_Ig2Vg_matrix()
+
+ 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
+
+
+class DipoleProportionalLoop():
+ """
+ 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
+
+ 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 DipoleTunerLoop():
+ """
+ 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 DipoleProportionalIntegralLoop():
+ """
+ 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.
+ 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).
+
+ 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)-->(-)-->(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)
+
+ # IIR filter
+ ## IIR_cutoff ; cutoff frequecny in Hz. If 0, cutoff frequency is infinity.
+ self.IIR_init(IIR_cutoff)
+ self.init_FFconst()
+
+ # Pre caclulation for Ig2Vg
+ 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]) * 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()
+
+ 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 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.
+ 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
+
+ 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))
+
+ 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
+
+ def IIR_init(self,cutoff):
+ """
+ f = 1/2/pi/T arccos((2-2r-r*r)/2(1-r))
+ T: Sample time
+ r: IIRcoef
+ f: cutoff freqeuncy
+ """
+ if cutoff == 0:
+ self.IIRcoef = 1.0
+ else:
+ omega = 2.0 * np.pi * cutoff
+ T = self.ring.T1 * self.every
+ alpha = np.cos(omega*T)-1
+ tmp = alpha * alpha - 2 * alpha
+ if tmp > 0:
+ self.IIRcoef = alpha + np.sqrt(tmp)
+ else:
+ self.IIRcoef = T * cutoff * 2 * np.pi
+ self.IIRout = self.cav_res.Vc
+
+ def IIR(self,input):
+ self.IIRout = (1-self.IIRcoef)*self.IIRout+self.IIRcoef*input
+ return self.IIRout
+
+ @property
+ def IIRcutoff(self):
+ """Return IIR cutoff frequency."""
+ T = self.ring.T1 * self.every
+ return 1.0/2.0/np.pi/T * np.arccos((2-2*self.IIRcoef-self.IIRcoef*self.IIRcoef)/2/(1-self.IIRcoef))
+
+class DipoleDirectFeedback(DipoleProportionalIntegralLoop):
+ """
+ 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(DipoleDirectFeedback, 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))
+
+ 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()
+
+ def track(self):
+ """
+ Tracking method for the Direct RF feedback.
+
+ Returns
+ -------
+ None.
+
+ """
+ super(DipoleDirectFeedback, 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
diff --git a/mbtrack2/utilities/beamloading.py b/mbtrack2/utilities/beamloading.py
index 1d208c2a7625349b31b2d254afbe738a0b57180e..bd0f7f9f7e68e71dde7ce165e5a630157aa9adef 100644
--- a/mbtrack2/utilities/beamloading.py
+++ b/mbtrack2/utilities/beamloading.py
@@ -295,10 +295,10 @@ class BeamLoadingEquilibrium():
x0 = x0 + [self.cavity_list[0].theta]
if CM:
- print("The initial center of mass offset is " +
+ print("# The initial center of mass offset is " +
str(self.center_of_mass() * 1e12) + " ps")
else:
- print("The initial energy balance is " +
+ print("# The initial energy balance is " +
str(self.energy_balance()) + " eV")
sol = root(lambda x: self.to_solve(x, CM),
@@ -319,12 +319,12 @@ class BeamLoadingEquilibrium():
self.PHI = sol.x[1::2]
if CM:
- print("The final center of mass offset is " +
+ print("# The final center of mass offset is " +
str(self.center_of_mass() * 1e12) + " ps")
else:
- print("The final energy balance is " + str(self.energy_balance()) +
+ print("# The final energy balance is " + str(self.energy_balance()) +
" eV")
- print("The algorithm has converged: " + str(sol.success))
+ print("# The algorithm has converged: " + str(sol.success))
if plot:
self.plot_rho(self.B1 / 4, self.B2 / 4)