Skip to content
Snippets Groups Projects
Commit 11b686ad authored by Alexis GAMELIN's avatar Alexis GAMELIN
Browse files

Update methods for ImpedanceModel and BeamLoadingEquilibrium 

The energy_loss and power_loss_spectrum methods for ImpedanceModel can now take any bunch spectrum in optional agrument.
Add bunch_spectrum and plot_bunch_spectrum to BeamLoadingEquilibrium.
Apply formatters and update notebook
parent 4f07219d
No related branches found
No related tags found
2 merge requests!13v0.7.0,!11Update methods for ImpedanceModel and BeamLoadingEquilibrium
Hide whitespace changes
Inline Side-by-side
examples/mbtrack2_cavity_resonator.ipynb
+ 139
− 147
View file @ 11b686ad
Source diff could not be displayed: it is too large. Options to address this: view the blob.
mbtrack2/impedance/impedance_model.py
+ 117
− 25
View file @ 11b686ad
...@@ -305,7 +305,7 @@ class ImpedanceModel(): ...@@ -305,7 +305,7 @@ class ImpedanceModel():
Z_type : str, optional Z_type : str, optional
Type of impedance to plot. Type of impedance to plot.
component : str, optional component : str, optional
Component to plot, can be "real" or "imag". Component to plot, can be "real" or "imag".
sigma : float, optional sigma : float, optional
RMS bunch length in [s] to use for the spectral density. If equal RMS bunch length in [s] to use for the spectral density. If equal
to None, the spectral density is not plotted. to None, the spectral density is not plotted.
...@@ -496,29 +496,51 @@ class ImpedanceModel(): ...@@ -496,29 +496,51 @@ class ImpedanceModel():
return summary return summary
def energy_loss(self, sigma, M, bunch_spacing, I, n_points=10e6): def energy_loss(self,
M,
bunch_spacing,
I,
sigma=None,
bunch_spectrum=None,
freq_spectrum=None,
n_points=10e6):
""" """
Compute the beam and bunch loss factor and energy losses for each type Compute the beam and bunch loss factor and energy losses for each type
of element in the model assuming Gaussian bunches and constant spacing of element in the model.
between bunches.
Assumtions:
- Constant spacing between bunches.
- All bunches have the same bunch distribution.
- Gaussian bunches if sigma is given.
Parameters Parameters
---------- ----------
sigma : float
RMS bunch length in [s].
M : int M : int
Number of bunches in the beam. Number of bunches in the beam.
bunch_spacing : float bunch_spacing : float
Time between two bunches in [s]. Time between two bunches in [s].
I : float I : float
Total beam current in [A]. Total beam current in [A].
sigma : float, optional
RMS bunch length in [s].
If None, freq_spectrum and bunch_spectrum must be given.
Default is None.
bunch_spectrum : array
Bunch spectrum to consider (i.e. FT of bunch profile).
Not used if sigma is not None.
Default is None.
freq_spectrum : array
Frequency points corresponding to bunch_spectrum in [Hz].
Not used if sigma is not None.
Default is None.
n_points : float, optional n_points : float, optional
Number of points used in the frequency spectrums. Number of points used in the frequency spectrums.
Default is 10e6.
Returns Returns
------- -------
summary : Dataframe summary : Dataframe
Contains the beam and bunch loss factor and energy loss for the Contains the beam and bunch loss factor and energy loss for the
full model and for each type of different component. full model and for each type of different component.
""" """
...@@ -532,9 +554,33 @@ class ImpedanceModel(): ...@@ -532,9 +554,33 @@ class ImpedanceModel():
fmin = -1 * fmax fmin = -1 * fmax
f = np.linspace(fmin, fmax, int(n_points)) f = np.linspace(fmin, fmax, int(n_points))
beam_spect = beam_spectrum(f, M, bunch_spacing, sigma=sigma) if sigma is None:
if bunch_spectrum is None:
bunch_spect = gaussian_bunch_spectrum(f, sigma) raise ValueError(
"Either sigma or bunch_spectrum and freq_spectrum" +
"must be specified.")
else:
if freq_spectrum is None:
raise ValueError(
"freq_spectrum must be given if bunch_spectrum" +
" is specified.")
if freq_spectrum[0] >= 0:
freq_spectrum = np.concatenate(
(-1 * freq_spectrum[:0:-1], freq_spectrum))
bunch_spectrum = np.concatenate(
(bunch_spectrum[:0:-1], bunch_spectrum))
bunch_spectrum_interp = interp1d(freq_spectrum,
bunch_spectrum,
bounds_error=False,
fill_value=0)
bunch_spect = bunch_spectrum_interp(f)
beam_spect = beam_spectrum(f,
M,
bunch_spacing,
bunch_spectrum=bunch_spect)
else:
beam_spect = beam_spectrum(f, M, bunch_spacing, sigma=sigma)
bunch_spect = gaussian_bunch_spectrum(f, sigma)
attr_list = self.sum_names attr_list = self.sum_names
...@@ -566,42 +612,60 @@ class ImpedanceModel(): ...@@ -566,42 +612,60 @@ class ImpedanceModel():
return summary return summary
def power_loss_spectrum(self, def power_loss_spectrum(self,
sigma,
M, M,
bunch_spacing, bunch_spacing,
I, I,
sigma=None,
bunch_spectrum=None,
freq_spectrum=None,
n_points=10e6, n_points=10e6,
max_overlap=False, max_overlap=False,
plot=False): plot=False):
""" """
Compute the power loss spectrum of the summed longitudinal impedance Compute the power loss spectrum of the summed longitudinal impedance
as in Eq. (4) of [1]. as in Eq. (4) of [1].
Assume Gaussian bunches and constant spacing between bunches. Assumtions:
- Constant spacing between bunches.
- All bunches have the same bunch distribution.
- Gaussian bunches if sigma is given.
Parameters Parameters
---------- ----------
sigma : float
RMS bunch length in [s].
M : int M : int
Number of bunches in the beam. Number of bunches in the beam.
bunch_spacing : float bunch_spacing : float
Time between two bunches in [s]. Time between two bunches in [s].
I : float I : float
Total beam current in [A]. Total beam current in [A].
sigma : float, optional
RMS bunch length in [s].
If None, freq_spectrum and bunch_spectrum must be given.
Default is None.
bunch_spectrum : array
Bunch spectrum to consider (i.e. FT of bunch profile).
Not used if sigma is not None.
Default is None.
freq_spectrum : array
Frequency points corresponding to bunch_spectrum in [Hz].
Not used if sigma is not None.
Default is None.
n_points : float, optional n_points : float, optional
Number of points used in the frequency spectrum. Number of points used in the frequency spectrum.
Default is 10e6.
max_overlap : bool, optional max_overlap : bool, optional
If True, the bunch spectrum (scaled to the number of bunches) is If True, the bunch spectrum (scaled to the number of bunches) is
used instead of the beam spectrum to compute the maximum value of used instead of the beam spectrum to compute the maximum value of
the power loss spectrum at each frequency. Should only be used to the power loss spectrum at each frequency. Should only be used to
maximise the power loss at a given frequency (e.g. for HOMs) and maximise the power loss at a given frequency (e.g. for HOMs) and
not for the full spectrum. not for the full spectrum.
Default is False.
plot : bool, optional plot : bool, optional
If True, plots: If True, plots:
- the overlap between the real part of the longitudinal impedance - the overlap between the real part of the longitudinal impedance
and the beam spectrum. and the beam spectrum.
- the power loss spectrum. - the power loss spectrum.
Default is False.
Returns Returns
------- -------
...@@ -612,8 +676,8 @@ class ImpedanceModel(): ...@@ -612,8 +676,8 @@ class ImpedanceModel():
References References
---------- ----------
[1] : L. Teofili, et al. "A Multi-Physics Approach to Simulate the RF [1] : L. Teofili, et al. "A Multi-Physics Approach to Simulate the RF
Heating 3D Power Map Induced by the Proton Beam in a Beam Intercepting Heating 3D Power Map Induced by the Proton Beam in a Beam Intercepting
Device", in IPAC'18, 2018, doi:10.18429/JACoW-IPAC2018-THPAK093 Device", in IPAC'18, 2018, doi:10.18429/JACoW-IPAC2018-THPAK093
""" """
...@@ -629,10 +693,38 @@ class ImpedanceModel(): ...@@ -629,10 +693,38 @@ class ImpedanceModel():
double_sided_impedance(impedance) double_sided_impedance(impedance)
frequency = np.linspace(fmin, fmax, int(n_points)) frequency = np.linspace(fmin, fmax, int(n_points))
if sigma is None:
if bunch_spectrum is None:
raise ValueError(
"Either sigma or bunch_spectrum and freq_spectrum" +
"must be specified.")
else:
if freq_spectrum is None:
raise ValueError(
"freq_spectrum must be given if bunch_spectrum" +
" is specified.")
if freq_spectrum[0] >= 0:
freq_spectrum = np.concatenate(
(-1 * freq_spectrum[:0:-1], freq_spectrum))
bunch_spectrum = np.concatenate(
(bunch_spectrum[:0:-1], bunch_spectrum))
bunch_spectrum_interp = interp1d(freq_spectrum,
bunch_spectrum,
bounds_error=False,
fill_value=0)
bunch_spectrum = bunch_spectrum_interp(frequency)
else:
bunch_spectrum = None
if max_overlap is False: if max_overlap is False:
spectrum = beam_spectrum(frequency, M, bunch_spacing, sigma) spectrum = beam_spectrum(frequency, M, bunch_spacing, sigma,
bunch_spectrum)
else: else:
spectrum = gaussian_bunch_spectrum(frequency, sigma) * M if bunch_spectrum is None:
spectrum = gaussian_bunch_spectrum(frequency, sigma) * M
else:
spectrum = bunch_spectrum * M
pmax = np.floor(fmax / self.ring.f0) pmax = np.floor(fmax / self.ring.f0)
pmin = np.floor(fmin / self.ring.f0) pmin = np.floor(fmin / self.ring.f0)
......
mbtrack2/utilities/beamloading.py
+ 39
− 23
View file @ 11b686ad
...@@ -6,6 +6,7 @@ Module where the BeamLoadingEquilibrium class is defined. ...@@ -6,6 +6,7 @@ Module where the BeamLoadingEquilibrium class is defined.
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from scipy.constants import c from scipy.constants import c
from scipy.fft import rfft, rfftfreq
from scipy.integrate import trapz from scipy.integrate import trapz
from scipy.optimize import root from scipy.optimize import root
...@@ -13,8 +14,8 @@ from mbtrack2.utilities.spectrum import gaussian_bunch ...@@ -13,8 +14,8 @@ from mbtrack2.utilities.spectrum import gaussian_bunch
class BeamLoadingEquilibrium(): class BeamLoadingEquilibrium():
"""Class used to compute beam equilibrium profile for a given storage ring """Class used to compute beam equilibrium profile for a given storage ring
and a list of RF cavities of any harmonic. The class assumes an uniform and a list of RF cavities of any harmonic. The class assumes an uniform
filling of the storage ring. Based on [1] which is an extension of [2]. filling of the storage ring. Based on [1] which is an extension of [2].
Parameters Parameters
...@@ -31,16 +32,16 @@ class BeamLoadingEquilibrium(): ...@@ -31,16 +32,16 @@ class BeamLoadingEquilibrium():
N : int or float, optinal N : int or float, optinal
Number of points used for longitudinal grid. Number of points used for longitudinal grid.
Default is 1000. Default is 1000.
References References
---------- ----------
[1] Gamelin, A., & Yamamoto, N. (2021). Equilibrium bunch density [1] Gamelin, A., & Yamamoto, N. (2021). Equilibrium bunch density
distribution with multiple active and passive RF cavities. distribution with multiple active and passive RF cavities.
IPAC'21 (MOPAB069). IPAC'21 (MOPAB069).
[2] Venturini, M. (2018). Passive higher-harmonic rf cavities with general [2] Venturini, M. (2018). Passive higher-harmonic rf cavities with general
settings and multibunch instabilities in electron storage rings. settings and multibunch instabilities in electron storage rings.
Physical Review Accelerators and Beams, 21(11), 114404. Physical Review Accelerators and Beams, 21(11), 114404.
""" """
def __init__(self, def __init__(self,
...@@ -73,6 +74,8 @@ class BeamLoadingEquilibrium(): ...@@ -73,6 +74,8 @@ class BeamLoadingEquilibrium():
self.z0 = np.linspace(self.B1, self.B2, self.N) self.z0 = np.linspace(self.B1, self.B2, self.N)
self.tau0 = self.z0 / c self.tau0 = self.z0 / c
self.rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0) / c self.rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0) / c
self.sampling = self.tau0[1] - self.tau0[0]
self.freq0 = rfftfreq(self.N, self.sampling)
# Define constants for scaled potential u(z) # Define constants for scaled potential u(z)
self.u0 = self.ring.U0 / (self.ring.ac * self.ring.sigma_delta**2 * self.u0 = self.ring.U0 / (self.ring.ac * self.ring.sigma_delta**2 *
...@@ -151,7 +154,7 @@ class BeamLoadingEquilibrium(): ...@@ -151,7 +154,7 @@ class BeamLoadingEquilibrium():
"""System of non-linear equation to solve to find the form factors F """System of non-linear equation to solve to find the form factors F
and PHI at equilibrum. and PHI at equilibrum.
The system is composed of Eq. (B6) and (B7) of [1] for each cavity. The system is composed of Eq. (B6) and (B7) of [1] for each cavity.
If auto_set_MC_theta == True, the system also find the main cavity If auto_set_MC_theta == True, the system also find the main cavity
phase to impose energy balance or cancel center of mass offset. phase to impose energy balance or cancel center of mass offset.
If CM is True, the system imposes zero center of mass offset, If CM is True, the system imposes zero center of mass offset,
if False, the system imposes energy balance. if False, the system imposes energy balance.
...@@ -269,7 +272,7 @@ class BeamLoadingEquilibrium(): ...@@ -269,7 +272,7 @@ class BeamLoadingEquilibrium():
"""Solve system of non-linear equation to find the form factors F """Solve system of non-linear equation to find the form factors F
and PHI at equilibrum. Can be used to compute the equilibrium bunch and PHI at equilibrum. Can be used to compute the equilibrium bunch
profile. profile.
Parameters Parameters
---------- ----------
x0 : initial guess x0 : initial guess
...@@ -280,7 +283,7 @@ class BeamLoadingEquilibrium(): ...@@ -280,7 +283,7 @@ class BeamLoadingEquilibrium():
CM : if True, the system imposes zero center of mass offset, CM : if True, the system imposes zero center of mass offset,
if False, the system imposes energy balance. if False, the system imposes energy balance.
verbiose : if True, print out informations about convergence. verbiose : if True, print out informations about convergence.
Returns Returns
------- -------
sol : OptimizeResult object representing the solution sol : OptimizeResult object representing the solution
...@@ -332,13 +335,13 @@ class BeamLoadingEquilibrium(): ...@@ -332,13 +335,13 @@ class BeamLoadingEquilibrium():
def PTBL_threshold(self, I0, m=None, MC_index=0, HC_index=1): def PTBL_threshold(self, I0, m=None, MC_index=0, HC_index=1):
""" """
Compute the periodic transient beam loading (PTBL) instability Compute the periodic transient beam loading (PTBL) instability
threshold based on [1]. threshold based on [1].
The beam_equilibrium method should be called before the PTBL_threshold The beam_equilibrium method should be called before the PTBL_threshold
one. one.
Eq. (22) and (23) have been modified for cosine voltage convention Eq. (22) and (23) have been modified for cosine voltage convention
for the main cavity. for the main cavity.
Parameters Parameters
...@@ -346,7 +349,7 @@ class BeamLoadingEquilibrium(): ...@@ -346,7 +349,7 @@ class BeamLoadingEquilibrium():
I0 : float I0 : float
Beam total current in [A]. Beam total current in [A].
m : int, optional m : int, optional
Number of bunches in the beam. Number of bunches in the beam.
The default is None. If None, uniform filling is assumed. The default is None. If None, uniform filling is assumed.
MC_index : int, optional MC_index : int, optional
Index of the main cavity in cavity_list. The default is 0. Index of the main cavity in cavity_list. The default is 0.
...@@ -356,18 +359,18 @@ class BeamLoadingEquilibrium(): ...@@ -356,18 +359,18 @@ class BeamLoadingEquilibrium():
Returns Returns
------- -------
eta : float eta : float
Amplification factor, if bigger than 1 then the beam is possibly Amplification factor, if bigger than 1 then the beam is possibly
in the PTBL regime, Eq. (22) of [1]. in the PTBL regime, Eq. (22) of [1].
RQth : float RQth : float
R/Q threshold, Eq. (24) of [1]. R/Q threshold, Eq. (24) of [1].
f : float f : float
f function, Eq. (23) of [1]. f function, Eq. (23) of [1].
Reference Reference
--------- ---------
[1] : He, T. (2022). Novel perturbation method for judging the [1] : He, T. (2022). Novel perturbation method for judging the
stability of the equilibrium solution in the presence of passive stability of the equilibrium solution in the presence of passive
harmonic cavities. Physical Review Accelerators and Beams, 25(9), harmonic cavities. Physical Review Accelerators and Beams, 25(9),
094402. 094402.
""" """
...@@ -402,14 +405,27 @@ class BeamLoadingEquilibrium(): ...@@ -402,14 +405,27 @@ class BeamLoadingEquilibrium():
def R_factor(self): def R_factor(self):
""" """
Touschek lifetime ratio as defined in [1]. Touschek lifetime ratio as defined in [1].
Reference Reference
--------- ---------
[1] : Byrd, J. M., and M. Georgsson. "Lifetime increase using passive [1] : Byrd, J. M., and M. Georgsson. "Lifetime increase using passive
harmonic cavities in synchrotron light sources." Physical Review harmonic cavities in synchrotron light sources." Physical Review
Special Topics-Accelerators and Beams 4.3 (2001): 030701. Special Topics-Accelerators and Beams 4.3 (2001): 030701.
""" """
rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0) / c rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0) / c
R = trapz(rho0**2, self.tau0) / trapz(self.rho0**2, self.tau0) R = trapz(rho0**2, self.tau0) / trapz(self.rho0**2, self.tau0)
return R return R
@property
def bunch_spectrum(self):
"""Return the bunch equilibrium spectrum."""
dft_rho = rfft(self.rho0) * self.sampling / self.N
return np.abs(dft_rho) / np.max(np.abs(dft_rho))
def plot_bunch_spectrum(self):
"""Plot the bunch equilibrium spectrum"""
plt.plot(self.freq0 * 1e-9, self.bunch_spectrum)
plt.xlabel(r"Frequency [GHz]")
plt.xlim([0, self.freq0.max() * 1e-9 / 3])
plt.title("Equilibrium bunch spectrum")
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

The place to be for your source code