diff --git a/mbtrack2/impedance/impedance_model.py b/mbtrack2/impedance/impedance_model.py
index 713fe464a7f336b543e3ba282748275c356433fa..c791a09e58e236e0de14bc84f9271a3ba8c4a428 100644
--- a/mbtrack2/impedance/impedance_model.py
+++ b/mbtrack2/impedance/impedance_model.py
@@ -305,7 +305,7 @@ class ImpedanceModel():
Z_type : str, optional
Type of impedance to plot.
component : str, optional
- Component to plot, can be "real" or "imag".
+ Component to plot, can be "real" or "imag".
sigma : float, optional
RMS bunch length in [s] to use for the spectral density. If equal
to None, the spectral density is not plotted.
@@ -496,29 +496,51 @@ class ImpedanceModel():
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
- of element in the model assuming Gaussian bunches and constant spacing
- between bunches.
+ Compute the beam and bunch loss factor and energy losses for each type
+ of element in the model.
+
+ Assumtions:
+ - Constant spacing between bunches.
+ - All bunches have the same bunch distribution.
+ - Gaussian bunches if sigma is given.
Parameters
----------
- sigma : float
- RMS bunch length in [s].
M : int
Number of bunches in the beam.
bunch_spacing : float
Time between two bunches in [s].
I : float
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
Number of points used in the frequency spectrums.
+ Default is 10e6.
Returns
-------
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.
"""
@@ -532,9 +554,24 @@ class ImpedanceModel():
fmin = -1 * fmax
f = np.linspace(fmin, fmax, int(n_points))
- beam_spect = beam_spectrum(f, M, bunch_spacing, sigma=sigma)
-
- bunch_spect = gaussian_bunch_spectrum(f, sigma)
+ 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_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
@@ -566,42 +603,60 @@ class ImpedanceModel():
return summary
def power_loss_spectrum(self,
- sigma,
M,
bunch_spacing,
I,
+ sigma=None,
+ bunch_spectrum=None,
+ freq_spectrum=None,
n_points=10e6,
max_overlap=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].
- 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
----------
- sigma : float
- RMS bunch length in [s].
M : int
Number of bunches in the beam.
bunch_spacing : float
Time between two bunches in [s].
I : float
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
Number of points used in the frequency spectrum.
+ Default is 10e6.
max_overlap : bool, optional
- If True, the bunch spectrum (scaled to the number of bunches) is
- used instead of the beam spectrum to compute the maximum value of
- 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
+ If True, the bunch spectrum (scaled to the number of bunches) is
+ used instead of the beam spectrum to compute the maximum value of
+ 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
not for the full spectrum.
+ Default is False.
plot : bool, optional
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.
- the power loss spectrum.
+ Default is False.
Returns
-------
@@ -612,8 +667,8 @@ class ImpedanceModel():
References
----------
- [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
+ [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
Device", in IPAC'18, 2018, doi:10.18429/JACoW-IPAC2018-THPAK093
"""
@@ -629,10 +684,31 @@ class ImpedanceModel():
double_sided_impedance(impedance)
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:
- spectrum = beam_spectrum(frequency, M, bunch_spacing, sigma)
+ spectrum = beam_spectrum(frequency, M, bunch_spacing, sigma, bunch_spectrum)
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)
pmin = np.floor(fmin / self.ring.f0)
diff --git a/mbtrack2/utilities/beamloading.py b/mbtrack2/utilities/beamloading.py
index 595fafd316940f8cdfa8f7d0cb53c4aa01f8187a..a74dec27f12e289ae85441ac0f527c519f15aa91 100644
--- a/mbtrack2/utilities/beamloading.py
+++ b/mbtrack2/utilities/beamloading.py
@@ -8,13 +8,14 @@ import numpy as np
from scipy.constants import c
from scipy.integrate import trapz
from scipy.optimize import root
+from scipy.fft import rfftfreq, rfft
from mbtrack2.utilities.spectrum import gaussian_bunch
class BeamLoadingEquilibrium():
- """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
+ """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
filling of the storage ring. Based on [1] which is an extension of [2].
Parameters
@@ -31,16 +32,16 @@ class BeamLoadingEquilibrium():
N : int or float, optinal
Number of points used for longitudinal grid.
Default is 1000.
-
+
References
----------
- [1] Gamelin, A., & Yamamoto, N. (2021). Equilibrium bunch density
- distribution with multiple active and passive RF cavities.
+ [1] Gamelin, A., & Yamamoto, N. (2021). Equilibrium bunch density
+ distribution with multiple active and passive RF cavities.
IPAC'21 (MOPAB069).
[2] 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.
-
+
"""
def __init__(self,
@@ -73,6 +74,8 @@ class BeamLoadingEquilibrium():
self.z0 = np.linspace(self.B1, self.B2, self.N)
self.tau0 = self.z0 / 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)
self.u0 = self.ring.U0 / (self.ring.ac * self.ring.sigma_delta**2 *
@@ -151,7 +154,7 @@ class BeamLoadingEquilibrium():
"""System of non-linear equation to solve to find the form factors F
and PHI at equilibrum.
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.
If CM is True, the system imposes zero center of mass offset,
if False, the system imposes energy balance.
@@ -269,7 +272,7 @@ class BeamLoadingEquilibrium():
"""Solve system of non-linear equation to find the form factors F
and PHI at equilibrum. Can be used to compute the equilibrium bunch
profile.
-
+
Parameters
----------
x0 : initial guess
@@ -280,7 +283,7 @@ class BeamLoadingEquilibrium():
CM : if True, the system imposes zero center of mass offset,
if False, the system imposes energy balance.
verbiose : if True, print out informations about convergence.
-
+
Returns
-------
sol : OptimizeResult object representing the solution
@@ -332,13 +335,13 @@ class BeamLoadingEquilibrium():
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].
-
- The beam_equilibrium method should be called before the PTBL_threshold
+
+ The beam_equilibrium method should be called before the PTBL_threshold
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.
Parameters
@@ -346,7 +349,7 @@ class BeamLoadingEquilibrium():
I0 : float
Beam total current in [A].
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.
MC_index : int, optional
Index of the main cavity in cavity_list. The default is 0.
@@ -356,18 +359,18 @@ class BeamLoadingEquilibrium():
Returns
-------
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].
RQth : float
R/Q threshold, Eq. (24) of [1].
f : float
f function, Eq. (23) of [1].
-
+
Reference
---------
- [1] : He, T. (2022). 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),
+ [1] : He, T. (2022). 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),
094402.
"""
@@ -402,14 +405,27 @@ class BeamLoadingEquilibrium():
def R_factor(self):
"""
Touschek lifetime ratio as defined in [1].
-
+
Reference
---------
- [1] : Byrd, J. M., and M. Georgsson. "Lifetime increase using passive
- harmonic cavities in synchrotron light sources." Physical Review
+ [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.
"""
rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0) / c
R = trapz(rho0**2, self.tau0) / trapz(self.rho0**2, self.tau0)
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")