# -*- coding: utf-8 -*-
"""
Module to compute longitudinal instabilities in a double RF system.
"""
import numpy as np
from math import factorial
from copy import deepcopy
from scipy.constants import c
from scipy.integrate import quad, trapezoid
from scipy.optimize import root, root_scalar
from mbtrack2 import BeamLoadingEquilibrium, gaussian_bunch
from pycolleff.longitudinal_equilibrium import LongitudinalEquilibrium, ImpedanceSource
# from robinsonmodes.pycolleff_to_mbtrack2 import synchrotron_to_pycolleff, cavityresonator_to_pycolleff, impedance_to_pycolleff
from albums.pycolleff_to_mbtrack2 import synchrotron_to_pycolleff, cavityresonator_to_pycolleff, impedance_to_pycolleff
[docs]
class RobinsonModes():
"""
Class to compute instabilities in a double RF system mainly based on the
theory and algorithm developped in [1].
The original algorithm has been adapted to work with other methods to solve
the self-consistent problem (i.e. Haissinski equation) [3,4,5,6] and to give
an estimate of the threshold of l=1/PTBL instability [7].
See the solve method for details.
Parameters
----------
ring : Synchrotron object
Ring parameters.
cavity_list : list of CavityResonator objects
A list which contains:
- a CavityResonator for the main cavity in 1st position.
- a CavityResonator for the harmonic cavity in 2nd position.
I0 : float
Beam current in [A].
tau_boundary : float, optinal
Integration boundary in [s].
If None, 0.1 * RF period is used.
The default is None.
Attributes
----------
xi : float
Methods
-------
solve()
Solve for a given settings of the double RF system.
R_factor(method)
Touschek lifetime ratio.
References
----------
[1] : Bosch, R. A., K. J. Kleman, and J. J. Bisognano. "Robinson
instabilities with a higher-harmonic cavity." Physical Review Special
Topics-Accelerators and Beams 4.7 (2001): 074401.
[2] : Bosch, R. A., and C. S. Hsue. "Suppression of longitudinal
coupled-bunch instabilities by a passive higher harmonic cavity."
Proceedings of International Conference on Particle Accelerators. IEEE, 1993.
[3] : Gamelin, A., Yamamoto, N. (2021). Equilibrium bunch density
distribution with multiple active and passive RF cavities.
IPAC'21 (MOPAB069).
[4] : Venturini, M. (2018). Passive higher-harmonic rf cavities with general
settings and multibunch instabilities in electron storage rings.
Physical Review Accelerators and Beams, 21(11), 114404.
[5] : Alves, Murilo B., and Fernando H. de Sá. "Equilibrium of longitudinal
bunch distributions in electron storage rings with arbitrary impedance
sources and generic filling patterns." Physical Review Accelerators and
Beams 26.9 (2023): 094402.
[6] : de Sá, F., & Alves, M. (2023). pycolleff and cppcolleff: modules for
impedance analysis and wake-field induced instabilities evaluation.
(Version 0.1.0) [Computer software]. https://doi.org/10.5281/zenodo.7974571
[7] : He, Tianlong. "Novel perturbation method for judging the stability of
the equilibrium solution in the presence of passive harmonic cavities."
Physical Review Accelerators and Beams 25.9 (2022): 094402.
"""
[docs]
def __init__(self,
ring,
cavity_list,
I0,
tau_boundary=None):
self.ring = ring
self.cavity_list = []
for item in cavity_list:
self.cavity_list.append(deepcopy(item))
self.I0 = I0
self.n_cavity = len(cavity_list)
if tau_boundary is None:
self.tau_boundary = self.ring.T1*0.1
else:
self.tau_boundary = tau_boundary
[docs]
def solve(self,
passive_harmonic_cavity=True,
mode_coupling=True,
auto_set_MC_theta=True,
optimal_tunning=True,
f_HOM=0,
Z_HOM=0,
method="Venturini",
**kwargs):
"""
Solve for a given settings of the double RF system.
Parameters
----------
passive_harmonic_cavity : bool, optional
If True, solve for a passive harmonic cavity.
If False, solve for an active harmonic cavity.
The default is True.
mode_coupling : bool, optional
If True, compute Robinson instabilities taking into account
dipole-quadrupole mode coupling (only dipole and quadrupole modes).
If False, compute Robinson instabilities without taking into
account mode coupling (modes from dipole to octupole).
The default is True.
auto_set_MC_theta : bool, optional
For passive harmonic cavity only:
-if True, set the main cavity phase to cancel the losses in the
harmonic cavity.
-if False, keep input main cavity phase setting.
The default is True.
optimal_tunning : bool, optional
For passive harmonic cavity:
-if True, set main cavity tuning to optimal.
-if False, keep input main cavity tuning setting.
For active harmonic cavity:
-if True, set main and harmonic cavity tuning to optimal.
-if False, keep input main and harmonic cavity tuning setting.
The default is True.
f_HOM : float or list of float, optional
Frequencies of the HOMs to evaluate in [Hz].
The default is 0.
Z_HOM : float or list of float, optional
Impedances of the HOMs at f_HOM in [ohm].
The default is 0.
method : {"Bosch","Venturini", "Alves"}, optional
Choose the method used to solve the self-consistent problem (i.e.
Haissinski equation):
- "Bosch" corresponds to an implementation of the original
algorithm in [1]. It computes only the amplitude form factor
and bunch length. The result can depend heavily on tau_boundary.
- "Venturini" corresponds to the method described in [3-4].
It can deal with passive and active HC.
- "Alves" corresponds to the method described in [5].
It can deal with arbitrary filling pattern and broad band impedance.
The pycolleff package is needed for this method [6].
Only for passive HC.
For both "Venturini" and "Alves", the bunch profile is computed
which gives both amplitude and phase form factors.
The default is "Venturini".
Kwargs for method "Bosch"
-------------------------
Passive harmonic cavity:
max_counter : int
Maximum number of iteration for the computation of form factors.
Default is 200.
Active harmonic cavity:
auto_set_for_xi : bool
If True, set harmonic cavity setting to achive xi value.
specified as in [1].
If False, keep input main and harmonic cavity setting.
Default is False.
xi : float
Value to achive if auto_set_for_xi is True.
Default is None.
Kwargs for method "Venturini"
-----------------------------
F_init : array-like of float of shape (2,)
Initial form factors (MC, HC) to consider for synchronous phase
calculation if auto_set_MC_theta is True.
Default is [1,1].
set_MC_phase_HCpassive : bool
If True, set the main cavity phase to cancel the losses in the
harmonic cavity using approxmiate relation.
Default is False.
Kwargs for method "Alves"
-------------------------
filling : array-like of float of shape (ring.h,)
Filling pattern of the beam.
Sum must be equal to one.
Default is None, corresponding to uniform filling.
impedance : Impedance object
Additional (broad band) longitudinal impedance to consider for the
Haissiski equation.
Default is None.
F_init : array-like of float of shape (2,)
Initial form factors (MC, HC) to consider for synchronous phase
calculation if auto_set_MC_theta is True.
Default is [1,1].
set_MC_phase_HCpassive : bool
If True, set the main cavity phase to cancel the losses in the
harmonic cavity using approximate relation.
Default is False.
use_GaussLMCI : bool
If True, the Gaussian LMCI method is used to compute Robinson
instabilities instead of Bosch's equations.
Default is False.
**kwargs passed to LongitudinalEquilibrium.calc_longitudinal_equilibrium:
- niter
- tol
- beta
- m
- print_flag
**kwargs passed to LongitudinalEquilibrium.calc_mode_coupling:
- max_azi
- max_rad
Returns
-------
bunch_length : float or None
RMS bunch length in [s].
zero_frequency : bool or None
True if zero frequency instability is predicted.
robinson : array of bool of shape (4,) or None
Bool of the arrays are True if the following instability is predicted:
If mode_coupling:
- Coupled dipole mode
- Coupled quadrupole mode
- Fast mode coupling by dipole mode
- Fast mode coupling by quadrupole mode
else:
- Dipole mode
- Quadrupole mode
- Sextupole mode
- Octupole mode
HOM : bool
True if HOM instability is predicted.
Omega : array of float or None
If mode_coupling:
- Coupled dipole angular frequency in unit of 2*pi*[Hz]
- Coupled quadrupole angular frequency in unit of 2*pi*[Hz]
else:
- Dipole angular frequency in unit of 2*pi*[Hz]
- Quadrupole angular frequency in unit of 2*pi*[Hz]
- Sextupole angular frequency in unit of 2*pi*[Hz]
- Octupole angular frequency in unit of 2*pi*[Hz]
PTBL : bool
True if PTBL instability is predicted.
None is returned if method is "Bosch".
converged : bool or array of bool
Return True values if calculation has converged. In details:
If form factor calculation or Haissinski equation resolution fails:
return False and other results are None
If zero_frequency is True:
return True and other instability results than zero_frequency
are None
If mode_coupling is True:
return array of bool of shape (2,) for coupled dipole and
quadrupole modes.
If mode_coupling is False:
return array of bool of shape (4,) for dipole to octupole modes.
"""
(bunch_length, a, b, c, converged) = self._converge_form_factors(passive_harmonic_cavity, auto_set_MC_theta, optimal_tunning, method, **kwargs)
if converged:
(zero_frequency, robinson, HOM, Omega, converged) = self._solve_core(bunch_length, a, b, c, mode_coupling, f_HOM, Z_HOM, method, **kwargs)
else:
return (None, None, None, None, None, None, converged)
if method == "Bosch":
omega_r = self._robinson_frequency()
landau = self._landau_threshold(omega_r, c, b)
coupled_mode1, _, _ = self._coupled_bunch_mode1(omega_r, bunch_length, landau)
return (bunch_length, zero_frequency, robinson, HOM, Omega, np.any(coupled_mode1), converged)
elif method == "Venturini":
(eta, _, _) = self.BLE.PTBL_threshold(self.I0)
if eta > 1:
PTBL = True
else:
PTBL = False
return (bunch_length, zero_frequency, robinson, HOM, Omega, PTBL, converged)
elif method == "Alves":
PTBL = self._PTBL_alves(**kwargs)
return (bunch_length, zero_frequency, robinson, HOM, Omega, PTBL, converged)
[docs]
def solve_equilibrium_only(self,
passive_harmonic_cavity=True,
auto_set_MC_theta=True,
optimal_tunning=True,
method="Venturini",
**kwargs):
(bunch_length, _, _, _, converged) = self._converge_form_factors(passive_harmonic_cavity, auto_set_MC_theta, optimal_tunning, method, **kwargs)
return (bunch_length, self.R_factor(method), self.xi, converged)
[docs]
def R_factor(self, method):
"""
Touschek lifetime ratio as defined in [1].
Parameters
----------
method : {"Bosch","Venturini", "Alves"}
Method used in self.solve.
Returns
-------
R : float
Touschek lifetime ratio for methods "Veturini" or "Alves".
1 if method is "Bosch" or solve has not converged.
Reference
---------
[1] : Byrd, J. M., and M. Georgsson. "Lifetime increase using passive
harmonic cavities in synchrotron light sources." Physical Review
Special Topics-Accelerators and Beams 4.3 (2001): 030701.
"""
if method == "Bosch":
R = 1
elif method == "Venturini":
try:
R = self.BLE.R_factor
except AttributeError:
R = 1
elif method == "Alves":
try:
rho0 = gaussian_bunch(self.LE.zgrid/c, self.ring.sigma_0)/c
R = trapezoid(rho0**2, self.LE.zgrid)/trapezoid(self.LE.distributions[0]**2, self.LE.zgrid)
except AttributeError:
R = 1
else:
raise AttributeError
return R
@property
def xi(self):
"""
Ratio of the harmonic cavity 'force' over the main cavity one.
Near Eq. (7) in [1].
"""
MC = self.cavity_list[0]
HC = self.cavity_list[1]
return -1 * HC.m * HC.Vc * np.sin(HC.theta) / (MC.Vc * np.sin(MC.theta))
[docs]
def _potential_decomposition(self):
# Eq. (4-6) in [1]
a_sum = 0
b_sum = 0
c_sum = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
a_sum += cavity.m*cavity.Vc*np.sin(cavity.theta)
b_sum += cavity.m**2*cavity.Vc*np.cos(cavity.theta)
c_sum += cavity.m**3*cavity.Vc*np.sin(cavity.theta)
a = self.ring.ac * self.ring.omega1 / (2 * self.ring.E0 * self.ring.T0) * a_sum
b = self.ring.ac * self.ring.omega1**2 / (6 * self.ring.E0 * self.ring.T0) * b_sum
c = - self.ring.ac * self.ring.omega1**3 / (24 * self.ring.E0 * self.ring.T0) * c_sum
return (a, b, c)
[docs]
def _bunch_length(self, a, b, c):
# Eq. (8) in [1]
U0 = self.ring.ac**2 * self.ring.sigma_delta**2/2
U = lambda tau: a*tau**2 + b*tau**3 + c*tau**4
numerator = lambda tau: tau**2*np.exp(-U(tau)/(2*U0))
denominator = lambda tau: np.exp(-U(tau)/(2.*U0))
num = quad(numerator, -self.tau_boundary, self.tau_boundary)[0]
den = quad(denominator, -self.tau_boundary, self.tau_boundary)[0]
bunch_length = np.sqrt(num/den)
return bunch_length
[docs]
def _robinson_frequency(self):
# Eq. (10) in [1]
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
sum_val += cavity.m*self.F[i]*cavity.Vc*np.sin(cavity.theta)
omega2 = self.ring.ac * self.ring.omega1 / (self.ring.E0 * self.ring.T0) * sum_val
return np.sqrt(omega2)
[docs]
def _dipole_coupled_bunch_growth_rate(self, f_HOM, Z_HOM, bunch_length, omega_r):
# Eq. (22) in [1]
from_factor = self._form_factors(f_HOM, 1, bunch_length)
gr = ((self.I0 * self.ring.eta() * 2 * np.pi * f_HOM * Z_HOM * from_factor**2) /
(2 * self.ring.E0 * self.ring.T0 * omega_r))
return gr
[docs]
def _dipole_coupled_bunch(self, f_HOM, Z_HOM, bunch_length, omega_r, landau):
HOM = False
if isinstance(f_HOM, (int, float)):
f_HOM = [f_HOM]
if isinstance(Z_HOM, (int, float)):
Z_HOM = [Z_HOM]
if len(f_HOM) != len(Z_HOM):
raise ValueError("f_HOM and Z_HOM must have the same length.")
else:
N = len(f_HOM)
for i in range(N):
gr_HOM = self._dipole_coupled_bunch_growth_rate(f_HOM[i], Z_HOM[i], bunch_length, omega_r)
# Check if growth rate is stronger than synchrotron damping
gr = gr_HOM - 1/self.ring.tau[2]
if gr > 0:
if gr_HOM > landau[0]:
HOM = True
break
return HOM
[docs]
def _zero_frequency_no_coupling(self, mu, bunch_length):
# Eq. (14) in [1]
coef1 = self.ring.ac * self.ring.omega1 * self.I0 / (self.ring.E0 * self.ring.T0)
coef2 = (self.ring.omega1 * bunch_length)**(2*mu-2)/(2**(mu-1)*factorial(mu))
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
sum_val += cavity.m**(2*mu - 1) * cavity.RL * self.F[i]**2 * np.sin(-2*cavity.psi)
return coef1 * coef2 * sum_val
[docs]
def _zero_frequency_instability(self, omega_r, bunch_length, mode_coupling):
# Eq. (14) in [1]
if mode_coupling is False:
if omega_r**2 < self._zero_frequency_no_coupling(1, bunch_length):
return True
elif omega_r**2 < self._zero_frequency_no_coupling(2, bunch_length):
return True
elif omega_r**2 < self._zero_frequency_no_coupling(3, bunch_length):
return True
elif omega_r**2 < self._zero_frequency_no_coupling(4, bunch_length):
return True
else:
return False
else:
# These conditions are the same for no coupling and with coupling (comparison with A_tilde and B_tilde results in same formula) except for addition of coupled zero-frequency condition
if omega_r**2 - self._A_tilde(0) < 0:
return True
elif 4*omega_r**2 - self._B_tilde(0, bunch_length) < 0:
return True
# Eq. (B.12) in [1]
elif (omega_r**2 - self._A_tilde(0))*(4*omega_r**2 - self._B_tilde(0, bunch_length)) + self._d_tilde(0, bunch_length) < 0:
return True
else:
return False
[docs]
def _landau_threshold(self, omega_r, c, b):
# Eq. (19) in [1]
dipole_landau = 0.78*self.ring.ac**2*self.ring.sigma_delta**2/omega_r*np.abs(3*c/omega_r**2 - (3*b/omega_r**2)**2)
landau_threshold = np.zeros((4,1))
landau_threshold[0] = dipole_landau
landau_threshold[1] = 2.24/0.78*dipole_landau
landau_threshold[2] = 4.12/0.78*dipole_landau
landau_threshold[3] = 6.36/0.78*dipole_landau
return landau_threshold
[docs]
def _robinson_damping_rate(self, Omega, mu, bunch_length):
# Eq. (16) in [1]
coef1 = 8 * self.ring.ac * self.I0 / (self.ring.E0 * self.ring.T0)
coef2 = mu * (self.ring.omega1 * bunch_length)**(2*mu - 2) / (2**mu * factorial(mu-1))
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
prod = np.tan(-cavity.psi) * np.cos(self._phi_pm(cavity, Omega, "+"))**2 * np.cos(self._phi_pm(cavity, Omega, "-"))**2
sum_val += cavity.m**(2*mu - 2) * self.F[i]**2 * cavity.RL * cavity.QL * prod
return coef1 * coef2 * sum_val
[docs]
def _coupled_bunch_mode1_damping_rate(self, Omega, mu, bunch_length):
# Eq. (15) in [1]
coef1 = self.ring.ac * self.ring.omega1 * self.I0 / (Omega * self.ring.E0 * self.ring.T0)
coef2 = mu*(self.ring.omega1 * bunch_length)**(2*mu - 2) / (2**mu * factorial(mu-1))
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
diff = np.cos(self._phi_pm(cavity, Omega, "-"))**2 - np.cos(self._phi_pm(cavity, Omega, "+"))**2
sum_val += cavity.m**(2*mu - 1) * self.F[i]**2 * cavity.RL * diff
return coef1 * coef2 * sum_val
[docs]
def _robinson_no_coupling(self, omega_r, bunch_length, landau_threshold):
# Follow III.A.(vii) in [1].
# Calculate AC robinson growth rates
modes = np.array([1,2,3,4])
# Output array
ac_robinson = np.zeros(len(modes),dtype=bool)
Omega_modes = np.zeros(len(modes),dtype=float)
converged = np.zeros(len(modes),dtype=bool)
# Calculate the Robinson frequency
for mu in modes:
# Initial guess based on zero current
Omega0 = mu*omega_r
# Eq. (13) in [1]
def f(Omega):
coef1 = self.ring.ac * self.ring.omega1 * self.I0 / (self.ring.E0 * self.ring.T0)
coef2 = mu * (self.ring.omega1 * bunch_length)**(2*mu - 2) / (2**mu * factorial(mu-1))
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.sin(2*self._phi_pm(cavity, Omega, "-")) + np.sin(2*self._phi_pm(cavity, Omega, "+"))
sum_val += cavity.m**(2*mu - 1) * cavity.RL * self.F[i]**2 * brackets
inner = (mu * omega_r)**2 - coef1 * coef2 * sum_val
return Omega - np.sqrt(inner)
sol = root_scalar(f, x0=Omega0, xtol=1e-6)
if sol.converged:
# Save the final result
Omega = sol.root
Omega_modes[mu-1] = Omega
converged[mu-1] = True
# When the Robinson frequency is found calculate the damping rate - Eq. (16) in [1]
alpha_r = self._robinson_damping_rate(Omega, mu, bunch_length)
# Add radiation damping to the damping rate
alpha_r_incl_rad_damp = alpha_r + mu/self.ring.tau[2]
# Complex frequency shift
delta_Omega = Omega - 1j*alpha_r - mu*omega_r
# Check if the instability occur
if alpha_r_incl_rad_damp < 0:
# Check if Landau damping is overcome
if np.abs(delta_Omega) > landau_threshold[mu-1]:
ac_robinson[mu-1] = True
else:
converged[mu-1] = False
return ac_robinson, Omega_modes, converged
[docs]
def _coupled_bunch_mode1(self, omega_r, bunch_length, landau_threshold):
# Follow III.A.(viii) in [1].
# Calculate coupled bunch mode +-1 growth rates
modes = ["+1","-1"]
mu = 1
# Output array
coupled_mode1 = np.zeros(len(modes),dtype=bool)
mode1_modes = np.zeros(len(modes),dtype=float)
converged = np.zeros(len(modes),dtype=bool)
# Calculate the coupled bunch mode +-1 frequency
for j, mode in enumerate(modes):
if mode == "+1":
added = self.ring.omega0
elif mode == "-1":
added = - self.ring.omega0
# Initial guess based on zero current
Omega0 = mu*omega_r
# Eq.13 with Omega -> Omega +- omega0
def f(Omega):
coef1 = self.ring.ac * self.ring.omega1 * self.I0 / (self.ring.E0 * self.ring.T0)
coef2 = mu * (self.ring.omega1 * bunch_length)**(2*mu - 2) / (2**mu * factorial(mu-1))
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.sin(2*self._phi_pm(cavity, Omega + added, "-")) + np.sin(2*self._phi_pm(cavity, Omega + added, "+"))
sum_val += cavity.m**(2*mu - 1) * cavity.RL * self.F[i]**2 * brackets
inner = (mu * omega_r)**2 - coef1 * coef2 * sum_val
return Omega - np.sqrt(inner)
sol = root_scalar(f, x0=Omega0, xtol=1e-6)
if sol.converged:
# Save the final result
Omega = sol.root
mode1_modes[j] = Omega
converged[j] = True
# When the coupled bunch mode +-1 frequency is found calculate the damping rate - Eq. 15
alpha_r = self._coupled_bunch_mode1_damping_rate(Omega + added, mu, bunch_length)
# Add radiation damping to the damping rate
alpha_r_incl_rad_damp = alpha_r + mu/self.ring.tau[2]
# Complex frequency shift
delta_Omega = Omega - 1j*alpha_r - mu*omega_r
# Check if the instability occur
if alpha_r_incl_rad_damp < 0:
# Check if Landau damping is overcome
if np.abs(delta_Omega) > landau_threshold[mu-1]:
coupled_mode1[j] = True
else:
converged[j] = False
return coupled_mode1, mode1_modes, converged
[docs]
def _robinson_coupling(self, omega_r, bunch_length, landau_threshold, abs_val=True):
# Follow III.A.(ix) in [1].
methods = ["hybr", "broyden1"] # These two methods seem enough to reach converge in most cases.
coupled_dipole = False
coupled_quadrupole = False
fast_mode_coupling_bydip = False
fast_mode_coupling_byquad = False
Omega_dip = np.nan
Omega_quad = np.nan
# Coupled dipole
def f_dipole(x):
y0 = x[0] - np.sqrt(self._Omega_B13(x[0], x[1], bunch_length, omega_r, "-", abs_val))
y1 = x[1] - self._ar_B11(x[0], bunch_length, omega_r)
return np.array([y0, y1])
# Initial guess based on zero current
x0 = np.array([omega_r, 1/self.ring.tau[2]])
for method in methods:
try:
sol = root(f_dipole, x0, tol=1e-6, method=method)
except ValueError:
sol = {"success":False}
if sol["success"]:
break
if sol["success"]:
# Save the final result
Omega_dip = sol.x[0]
alpha_r_incl_rad_damp_dip = sol.x[1]
delta_Omega_dip = Omega_dip - 1j*alpha_r_incl_rad_damp_dip - omega_r
converged_dip = True
# Check fast mode-coupling instability Eq. inside square root of B13
root_val = self._Omega_B13(Omega_dip, alpha_r_incl_rad_damp_dip, bunch_length, omega_r, "-", return_root=True)
if root_val < 0:
fast_mode_coupling_bydip = True
# Check if coupled-dipole instability exists
if alpha_r_incl_rad_damp_dip < 0:
# Check if Landau damping is overcome
if np.abs(delta_Omega_dip) > landau_threshold[0]:
coupled_dipole = True
else:
converged_dip = False
# Coupled quadrupole
def f_quadrupole(x):
y0 = x[0] - np.sqrt(self._Omega_B13(x[0], x[1], bunch_length, omega_r, "+", abs_val))
y1 = x[1] - self._ar_B11(x[0], bunch_length, omega_r)
return np.array([y0, y1])
# Initial guess based on zero current
x0 = np.array([2*omega_r, 2/self.ring.tau[2]])
for method in methods:
try:
sol = root(f_quadrupole, x0, tol=1e-6, method=method)
except ValueError:
sol = {"success":False}
if sol["success"]:
break
if sol["success"]:
# Save the final result
Omega_quad = sol.x[0]
alpha_r_incl_rad_damp_quad = sol.x[1]
delta_Omega_quad = Omega_quad - 1j*alpha_r_incl_rad_damp_quad - 2*omega_r
converged_quad = True
# Check fast mode-coupling instability Eq. inside square root of B13
root_val = self._Omega_B13(Omega_quad, alpha_r_incl_rad_damp_quad, bunch_length, omega_r, "+", return_root=True)
if root_val < 0:
fast_mode_coupling_byquad = True
# Check if coupled-quadrupole instability exists
if alpha_r_incl_rad_damp_quad < 0:
# Check if Landau damping is overcome
if np.abs(delta_Omega_quad) > landau_threshold[1]:
coupled_quadrupole = True
else:
converged_quad = False
robinson = np.array([coupled_dipole, coupled_quadrupole, fast_mode_coupling_bydip, fast_mode_coupling_byquad])
Omega = np.array([Omega_dip, Omega_quad])
converged = np.array([converged_dip, converged_quad])
return (robinson, Omega, converged)
[docs]
def _conv_ff_passive_bosch(self,
auto_set_MC_theta=True,
optimal_tunning=True,
**kwargs
):
# Follow III.A. from (i) to (v) in [1] with added information from [2].
max_counter = kwargs.get("max_counter", 200)
self.F = np.zeros(self.n_cavity)
self.F[0] = 1
self.F[1:] = 0.1
MC = self.cavity_list[0]
diff = 1
counter = 0
max_counter = max_counter
converged = True
old_F = np.zeros_like(self.F)
while(diff > 0.01):
counter += 1
if counter > max_counter:
converged = False
break
if auto_set_MC_theta:
# Calculate synchronous phase
delta = 0
for i in range(self.n_cavity):
if i == 0:
continue
cavity = self.cavity_list[i]
delta += cavity.Vb(self.I0) * self.F[i] * np.cos(cavity.psi)
MC.theta = np.arccos((self.ring.U0 + delta) / MC.Vc)
if optimal_tunning:
# Optimal detuning for main cavity
MC.set_optimal_detune(self.I0, self.F[0])
# Get theta & Vc (from (iii) near Eq. (21) in [1])
for i in range(self.n_cavity):
if i == 0:
continue
cavity = self.cavity_list[i]
cavity.theta = np.arctan(np.sin(-2*cavity.psi)/(-2*np.cos(cavity.psi)**2))
cavity.Vc = -2*self.I0*self.F[i]*cavity.RL*np.cos(cavity.psi)**2 / np.cos(cavity.theta)
# Taylor expansion coeff
a, b, c = self._potential_decomposition()
# Bunch length
bunch_length = self._bunch_length(a, b, c)
# Form factors
diff = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
old_F[i] = self.F[i]
self.F[i] = self._form_factors(self.ring.f1, cavity.m, bunch_length)
diff += np.abs((self.F[i] - old_F[i])/self.F[i])
if diff > 0.01:
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
# From 4. in [2].
self.F[i] = (self.F[i]*0.1 + old_F[i]*0.9)
return (bunch_length, a, b, c, converged)
[docs]
def _conv_ff_active_bosch(self, optimal_tunning=True, **kwargs):
# Follow V.A. from (i) to (iv) in [1].
# kwargs: auto_set_for_xi, xi
auto_set_for_xi = kwargs.get("auto_set_for_xi", False)
self.F = np.zeros(self.n_cavity)
self.F[0] = 1
self.F[1:] = 0.1
if auto_set_for_xi:
xi = kwargs.get("xi", None)
if xi is None:
raise ValueError("xi must be given if auto_set_for_xi=True")
MC = self.cavity_list[0]
HC = self.cavity_list[1]
MC.theta = np.arccos(self.ring.U0/(MC.Vc*(1 - 1/HC.m**2)))
HC.theta = np.arctan(HC.m * xi * np.tan(MC.theta)) - np.pi
HC.Vc = - xi * MC.Vc * np.sin(MC.theta) / (HC.m * np.sin(HC.theta))
(a, b, c) = self._potential_decomposition()
bunch_length = self._bunch_length(a, b, c)
for i, cavity in enumerate(self.cavity_list):
self.F[i] = self._form_factors(self.ring.f1, cavity.m, bunch_length)
if optimal_tunning:
for i, cavity in enumerate(self.cavity_list):
cavity.set_optimal_detune(self.I0, self.F[i])
return (bunch_length, a, b, c, True)
[docs]
def _conv_ff_venturini(self,
passive_harmonic_cavity=True,
auto_set_MC_theta=True,
optimal_tunning=True,
**kwargs):
self.F = kwargs.get("F_init", np.ones(self.n_cavity))
set_MC_phase_HCpassive = kwargs.get("set_MC_phase_HCpassive", False)
MC = self.cavity_list[0]
HC = self.cavity_list[1]
if set_MC_phase_HCpassive:
# Calculate synchronous phase
delta = 0
for i in range(self.n_cavity):
if i == 0:
continue
cavity = self.cavity_list[i]
delta += cavity.Vb(self.I0) * self.F[i] * np.cos(cavity.psi)
MC.theta = np.arccos((self.ring.U0 + delta) / MC.Vc)
if optimal_tunning:
# Optimal detuning for main cavity
MC.set_optimal_detune(self.I0)
if optimal_tunning and not passive_harmonic_cavity:
# Optimal detuning for active harmonic cavity
HC.set_optimal_detune(self.I0)
MC.set_generator(self.I0)
if passive_harmonic_cavity:
HC.Vg = 0
HC.theta_g = 0
else:
HC.set_generator(self.I0)
BLE = BeamLoadingEquilibrium(self.ring, self.cavity_list, self.I0,
auto_set_MC_theta=auto_set_MC_theta,
B1=-1*self.tau_boundary*c,
B2=self.tau_boundary*c)
try:
sol = BLE.beam_equilibrium()
bunch_length = BLE.std_rho()/c
converged = sol.success
self.F = BLE.F
self.PHI = BLE.PHI
self.BLE = BLE
# Get theta & Vc from phasor addition to take into account F&PHI and generator voltage
for i, cavity in enumerate(self.cavity_list):
Vc_phasor = -1*cavity.Vb(self.I0) * self.F[i] * np.exp(1j*(cavity.psi - self.PHI[i])) + cavity.Vg * np.exp(1j*cavity.theta_g)
cavity.Vc = np.abs(Vc_phasor)
cavity.theta = np.angle(Vc_phasor)
# Taylor expansion coeff
a, b, cc = self._potential_decomposition()
except:
converged = False
bunch_length = None
a = None
b = None
cc = None
return (bunch_length, a, b, cc, converged)
[docs]
def _conv_ff_alves(self,
passive_harmonic_cavity=True,
auto_set_MC_theta=True,
optimal_tunning=True,
**kwargs):
if passive_harmonic_cavity is False:
raise NotImplementedError("Alves method with active harmonic cavities is not implemented.")
# kwargs
niter = kwargs.get('niter', 200)
tol = kwargs.get('tol', 1e-8)
beta = kwargs.get('beta', 0.1)
m = kwargs.get('m', 3)
print_flag = kwargs.get('print_flag', False)
filling = kwargs.get("filling", None)
impedance = kwargs.get("impedance", None)
self.F = kwargs.get("F_init", np.ones(self.n_cavity))
self.PHI = np.zeros(self.n_cavity)
set_MC_phase_HCpassive = kwargs.get("set_MC_phase_HCpassive", False)
MC = self.cavity_list[0]
HC = self.cavity_list[1]
if set_MC_phase_HCpassive:
# Calculate synchronous phase
delta = 0
for i in range(self.n_cavity):
if i == 0:
continue
cavity = self.cavity_list[i]
delta += cavity.Vb(self.I0) * self.F[i] * np.cos(cavity.psi)
MC.theta = np.arccos((self.ring.U0 + delta) / MC.Vc)
if optimal_tunning:
# Optimal detuning for main cavity
MC.set_optimal_detune(self.I0)
MC.set_generator(self.I0)
HC.Vg = 0
HC.theta_g = 0
if filling is None:
filling = np.ones(self.ring.h) / self.ring.h
identical_bunches = True
else:
identical_bunches = False
nb = sum(filling != 0)
ring_pce = synchrotron_to_pycolleff(self.ring, self.I0, MC.Vc, nb)
mc_pce = cavityresonator_to_pycolleff(MC, Impedance=False)
hc_pce = cavityresonator_to_pycolleff(HC, Impedance=True)
sources = [mc_pce, hc_pce]
if impedance is not None:
impedance_pce = impedance_to_pycolleff(impedance)
sources += [impedance_pce]
longeq = LongitudinalEquilibrium(ring=ring_pce,
impedance_sources=sources,
fillpattern=filling)
longeq.identical_bunches = identical_bunches
longeq.feedback_on = True
longeq.feedback_method = 0 # 0 = Phasor, 1 = LS
# Numerical parameters from M. Alves example
longeq.zgrid = np.linspace(-1, 1, 2001) * ring_pce.rf_lamb / 2
longeq.max_mode = 1000*ring_pce.harm_num
longeq.min_mode0_ratio = 1e-10
try:
_, converged = longeq.calc_longitudinal_equilibrium(niter,
tol,
beta,
m,
print_flag)
z0, sigmaz = longeq.calc_moments(longeq.zgrid, longeq.distributions)
bunch_length = sigmaz[0]/c
if auto_set_MC_theta:
ring_pce = ring_pce = synchrotron_to_pycolleff(self.ring, self.I0, MC.Vc, nb)
longeq = LongitudinalEquilibrium(ring=ring_pce,
impedance_sources=sources,
fillpattern=filling)
longeq.identical_bunches = identical_bunches
longeq.main_ref_phase_offset = z0[0]/c*self.ring.omega1
longeq.feedback_on = True
longeq.feedback_method = 0 # 0 = Phasor, 1 = LS
longeq.zgrid = np.linspace(-1, 1, 2001) * ring_pce.rf_lamb / 2
longeq.max_mode = 1000*ring_pce.harm_num
longeq.min_mode0_ratio = 1e-10
_, converged = longeq.calc_longitudinal_equilibrium(niter,
tol,
beta,
m,
print_flag)
except:
converged = False
if converged:
z0, sigmaz = longeq.calc_moments(longeq.zgrid, longeq.distributions)
bunch_length = sigmaz[0]/c
self.LE = longeq
self.LE.ring.bunlen = sigmaz[0]
# FF
zgrid = longeq.zgrid - z0[0]
rho = longeq.distributions[0]
for i, cavity in enumerate(self.cavity_list):
exp = np.exp(1j * self.ring.k1 * cavity.m * zgrid)
FF = trapezoid(exp * rho, x=zgrid)
self.F[i] = np.real(FF)
self.PHI[i] = np.imag(FF)
# Get Vg and theta_g
MC.Vg = longeq.main_gen_amp_mon
MC.theta_g = -1*longeq.main_gen_phase_mon
# Get theta & Vc from phasor addition to take into account F&PHI and generator voltage
# (instead of (iii) near Eq. 21)
for i, cavity in enumerate(self.cavity_list):
Vc_phasor = -1*cavity.Vb(self.I0) * self.F[i] * np.exp(1j*(cavity.psi - self.PHI[i])) + cavity.Vg * np.exp(1j*cavity.theta_g)
cavity.Vc = np.abs(Vc_phasor)
cavity.theta = np.angle(Vc_phasor)
# Taylor expansion coeff
a, b, cc = self._potential_decomposition()
else:
bunch_length = None
a = None
b = None
cc = None
return (bunch_length, a, b, cc, converged)
[docs]
def _solve_core(self, bunch_length, a, b, c, mode_coupling, f_HOM, Z_HOM, method, **kwargs):
use_GaussLMCI = kwargs.get('use_GaussLMCI', False)
zero_frequency = None
robinson = None
HOM = None
Omega = None
omega_r = self._robinson_frequency()
landau = self._landau_threshold(omega_r, c, b)
HOM = self._dipole_coupled_bunch(f_HOM, Z_HOM, bunch_length, omega_r, landau)
zero_frequency = self._zero_frequency_instability(omega_r, bunch_length, mode_coupling)
if zero_frequency:
return (zero_frequency, robinson, HOM, Omega, True)
if method == "Alves" and use_GaussLMCI:
robinson, Omega, converged = self._FMCI_alves(**kwargs)
else:
if mode_coupling:
robinson, Omega, converged = self._robinson_coupling(omega_r, bunch_length, landau)
else:
robinson, Omega, converged = self._robinson_no_coupling(omega_r, bunch_length, landau)
return (zero_frequency, robinson, HOM, Omega, converged)
[docs]
def _omega_n(self, cavity):
# cavity resonant frequency - in text near II.A in [1]
return 2 * cavity.QL * cavity.m * self.ring.omega1 / (2 * cavity.QL + np.tan(- cavity.psi))
[docs]
def _phi_pm(self, cavity, Omega, sign):
# Near Eq. (13) in [1]
omega_n = self._omega_n(cavity)
if sign == "+":
return np.arctan(2*cavity.QL*(cavity.m*self.ring.omega1+Omega-omega_n)/omega_n)
elif sign == "-":
return np.arctan(2*cavity.QL*(cavity.m*self.ring.omega1-Omega-omega_n)/omega_n)
[docs]
def _A_tilde(self, Omega):
# Eq. (B8) in [1]
coef = self.ring.ac * self.ring.omega1 * self.I0 / (2* self.ring.E0 * self.ring.T0)
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.sin(2*self._phi_pm(cavity, Omega, "-")) + np.sin(2*self._phi_pm(cavity, Omega, "+"))
sum_val += cavity.m * cavity.RL * self.F[i]**2 * brackets
return coef * sum_val
[docs]
def _B_tilde(self, Omega, bunch_length):
# Eq. (B8) in [1]
coef1 = self.ring.ac * self.ring.omega1 * self.I0 / (2* self.ring.E0 * self.ring.T0)
coef2 = (self.ring.omega1 * bunch_length)**2
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.sin(2*self._phi_pm(cavity, Omega, "-")) + np.sin(2*self._phi_pm(cavity, Omega, "+"))
sum_val += cavity.m**3 * cavity.RL * self.F[i]**2 * brackets
return coef1 * coef2 * sum_val
[docs]
def _D_tilde(self, Omega, bunch_length):
# Eq. (B8) in [1]
coef2 = self.ring.ac * self.ring.omega1 * self.I0 / (2* self.ring.E0 * self.ring.T0)
coef1 = (self.ring.omega1 * bunch_length)
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.sin(2*self._phi_pm(cavity, Omega, "-")) - np.sin(2*self._phi_pm(cavity, Omega, "+"))
sum_val += cavity.m**2 * cavity.RL * self.F[i]**2 * brackets
return coef1 * coef2 * sum_val
[docs]
def _a_tilde(self, Omega):
# Eq. (B8) in [1]
coef = self.ring.ac * self.ring.omega1 * self.I0 / (self.ring.E0 * self.ring.T0)
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.cos(self._phi_pm(cavity, Omega, "-"))**2 - np.cos(self._phi_pm(cavity, Omega, "+"))**2
sum_val += cavity.m * cavity.RL * self.F[i]**2 * brackets
return coef * sum_val + 2 * Omega / self.ring.tau[2]
[docs]
def _b_tilde(self, Omega, bunch_length):
# Eq. (B8) in [1]
coef1 = self.ring.ac * self.ring.omega1 * self.I0 / (self.ring.E0 * self.ring.T0)
coef2 = (self.ring.omega1 * bunch_length)**2
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.cos(self._phi_pm(cavity, Omega, "-"))**2 - np.cos(self._phi_pm(cavity, Omega, "+"))**2
sum_val += cavity.m**3 * cavity.RL * self.F[i]**2 * brackets
return coef1 * coef2 * sum_val + 4 * Omega / self.ring.tau[2]
[docs]
def _d_tilde(self, Omega, bunch_length):
# Eq. (B8) in [1]
coef2 = self.ring.ac * self.ring.omega1 * self.I0 / (self.ring.E0 * self.ring.T0)
coef1 = (self.ring.omega1 * bunch_length)
sum_val = 0
for i in range(self.n_cavity):
cavity = self.cavity_list[i]
brackets = np.cos(self._phi_pm(cavity, Omega, "-"))**2 + np.cos(self._phi_pm(cavity, Omega, "+"))**2
sum_val += cavity.m**2 * cavity.RL * self.F[i]**2 * brackets
return coef1 * coef2 * sum_val
[docs]
def _ar_B11(self, Omega, bunch_length, omega_r):
# Eq. (B11) in [1]
a = self._a_tilde(Omega)
b = self._b_tilde(Omega, bunch_length)
d = self._d_tilde(Omega, bunch_length)
A = self._A_tilde(Omega)
B = self._B_tilde(Omega, bunch_length)
D = self._D_tilde(Omega, bunch_length)
num = a * (Omega**2 - (2 * omega_r)**2 + B) + b*(Omega**2 - omega_r**2 + A) - 2*D*d
den = 2 * Omega * (2 * Omega**2 - 5*omega_r**2 + A + B)
return num/den
[docs]
def _Omega_B13(self, Omega, alpha_r, bunch_length, omega_r, sign, abs_val=True, return_root=False):
# Eq. (B13) in [1]
a = self._a_tilde(Omega)
b = self._b_tilde(Omega, bunch_length)
d = self._d_tilde(Omega, bunch_length)
A = self._A_tilde(Omega)
B = self._B_tilde(Omega, bunch_length)
D = self._D_tilde(Omega, bunch_length)
part1 = (5*omega_r**2 - A - B)/2
root_val = (3*omega_r**2 + A - B)**2 / 4 + D**2 - d**2 + (a - 2*Omega*alpha_r)*(b - 2*Omega*alpha_r)
if return_root:
return root_val
if abs_val:
root_val = np.abs(root_val)
if sign == "+":
return part1 + np.sqrt(root_val)
elif sign == "-":
return part1 - np.sqrt(root_val)
[docs]
def _PTBL_alves(self, **kwargs):
# kwargs
max_azi = kwargs.get("max_azi", 10)
max_rad = kwargs.get("max_rad", 10)
# Synchrotron tune from bunch length and maching condition
sync_freq = self.ring.sigma_delta * self.ring.eta() * c
sync_freq /= self.LE.ring.bunlen * 2 * np.pi
self.LE.ring.sync_tune = sync_freq / self.ring.f0
# All impedance source as ImpedanceDFT
for idx, _ in enumerate(self.LE.impedance_sources):
self.LE.impedance_sources[idx].calc_method = ImpedanceSource.Methods.ImpedanceDFT
# Calculate Vlasov's equation eigen-frequencies:
eigenfreq, *_ = self.LE.calc_mode_coupling(
w=[-10 * self.ring.omega1, +10 * self.ring.omega1],
cbmode=1,
max_azi=max_azi,
max_rad=max_rad,
use_fokker=False,
delete_m0=True,
delete_m0k0=True,
reduced=True,
)
# Find most unstable mode:
idx = np.argmax(eigenfreq.imag)
grate = eigenfreq.imag[idx]
# Subtract radiation damping rate
grate -= 1 / self.ring.tau[2]
PTBL = grate > 0
return PTBL
[docs]
def _FMCI_alves(self, **kwargs):
# kwargs
max_azi = kwargs.get("max_azi", 2)
max_rad = kwargs.get("max_rad", 0)
# Synchrotron tune from bunch length and maching condition
sync_freq = self.ring.sigma_delta * self.ring.eta() * c
sync_freq /= self.LE.ring.bunlen * 2 * np.pi
self.LE.ring.sync_tune = sync_freq / self.ring.f0
# All impedance source as ImpedanceDFT
for idx, _ in enumerate(self.LE.impedance_sources):
self.LE.impedance_sources[idx].calc_method = ImpedanceSource.Methods.ImpedanceDFT
# Calculate Vlasov's equation eigen-frequencies:
eigenfreq, *_ = self.LE.calc_mode_coupling(
w=[-300 * self.ring.omega1, +300 * self.ring.omega1],
cbmode=0,
max_azi=max_azi,
max_rad=max_rad,
use_fokker=False,
delete_m0=True,
delete_m0k0=True,
reduced=True,
)
# Find most unstable mode:
idx = np.argmax(eigenfreq.imag)
grate = eigenfreq.imag[idx]
# Subtract radiation damping rate
grate -= 1 / self.ring.tau[2]
FMCI = grate > 0
robinson = np.array([False, False, FMCI, FMCI])
Omega = np.array([eigenfreq.real[0], eigenfreq.real[1]])
converged = np.array([True, True])
return (robinson, Omega, converged)