diff --git a/mbtrack2/impedance/impedance_model.py b/mbtrack2/impedance/impedance_model.py
index bc3da3c0f8f7b3bb2794e580668c1efdc1783239..cf586387aea1e5dd7d78f23ee23ec373df9eab9a 100644
--- a/mbtrack2/impedance/impedance_model.py
+++ b/mbtrack2/impedance/impedance_model.py
@@ -102,8 +102,6 @@ class ImpedanceModel():
None.
"""
- self.wakefields.append(wakefield)
- self.positions.append(positions)
if name is None:
name = wakefield.name
if name is None:
@@ -112,6 +110,8 @@ class ImpedanceModel():
self.names.append(name)
else:
raise ValueError("This name is already taken.")
+ self.wakefields.append(wakefield)
+ self.positions.append(positions)
def add_global(self, wakefield, name=None):
"""
diff --git a/mbtrack2/impedance/wakefield.py b/mbtrack2/impedance/wakefield.py
index 5de4a8fffb9ecf70a7feb9f14a1a1347ba8d1e5a..31712518626ca947d6faedfee0219d01d1a90876 100644
--- a/mbtrack2/impedance/wakefield.py
+++ b/mbtrack2/impedance/wakefield.py
@@ -277,6 +277,8 @@ class WakeFunction(ComplexData):
Compute a wake function from wake potential data.
plot()
Plot the wake function data.
+ loss_factor(sigma)
+ Compute the loss factor or the kick factor for a Gaussian bunch.
"""
def __init__(self,
@@ -454,6 +456,34 @@ class WakeFunction(ComplexData):
label = r"$W_{" + self.component_type + r"}$" + unit
ax.set_ylabel(label)
return ax
+
+ def loss_factor(self, sigma):
+ """
+ Compute the loss factor or the kick factor for a Gaussian bunch.
+ Formulas from Eq. (5.6), Eq. (5.12) and Eq. (5.17) of [1].
+
+ Parameters
+ ----------
+ sigma : float
+ RMS bunch length in [s]
+
+ Returns
+ -------
+ kloss: float
+ Loss factor in [V/C] or kick factor in [V/C/m] depanding on the
+ impedance type.
+
+ References
+ ----------
+ [1] : Zotter, Bruno W., and Semyon A. Kheifets. Impedances and wakes
+ in high-energy particle accelerators. World Scientific, 1998.
+
+ """
+ time = self.data.index
+ S = 1/(2*np.sqrt(np.pi)*sigma)*np.exp(-1*time**2/(4*sigma**2))
+ kloss = trapz(x = time, y = self.data["real"]*S)
+
+ return kloss
class Impedance(ComplexData):
"""
diff --git a/mbtrack2/tracking/parallel.py b/mbtrack2/tracking/parallel.py
index 30f835b413adfd02aef6a9f4aeabf860f0fa0cc5..b71cf4d66d6f5661bee12c17a177f02c8cf689cb 100644
--- a/mbtrack2/tracking/parallel.py
+++ b/mbtrack2/tracking/parallel.py
@@ -48,13 +48,16 @@ class Mpi:
Compute the bunch profiles and share it between the different bunches.
share_means(beam)
Compute the bunch means and share it between the different bunches.
+ share_stds(beam)
+ Compute the bunch standard deviations and share it between the
+ different bunches.
References
----------
[1] L. Dalcin, P. Kler, R. Paz, and A. Cosimo, Parallel Distributed
Computing using Python, Advances in Water Resources, 34(9):1124-1139, 2011.
- """
+ """
def __init__(self, filling_pattern):
from mpi4py import MPI
self.MPI = MPI
@@ -214,4 +217,29 @@ class Mpi:
else:
mean = np.zeros((6,), dtype=np.float64)
self.comm.Allgather([mean, self.MPI.DOUBLE], [self.mean_all, self.MPI.DOUBLE])
+
+ def share_stds(self, beam):
+ """
+ Compute the bunch standard deviations and share it between the
+ different bunches.
+
+ Parameters
+ ----------
+ beam : Beam object
+
+ """
+ if(beam.mpi_switch == False):
+ print("Error, mpi is not initialised.")
+
+ bunch = beam[self.bunch_num]
+
+ charge_all = self.comm.allgather(bunch.charge)
+ self.charge_all = charge_all
+
+ self.std_all = np.empty((self.size, 6), dtype=np.float64)
+ if len(bunch) != 0:
+ std = bunch.std
+ else:
+ std = np.zeros((6,), dtype=np.float64)
+ self.comm.Allgather([std, self.MPI.DOUBLE], [self.std_all, self.MPI.DOUBLE])
\ No newline at end of file
diff --git a/mbtrack2/tracking/particles.py b/mbtrack2/tracking/particles.py
index 9dfd6910fb1baa41950d9dea10eaf7b4b3063740..e663ee0ac563c871951452098a1d4565a1226237 100644
--- a/mbtrack2/tracking/particles.py
+++ b/mbtrack2/tracking/particles.py
@@ -119,12 +119,12 @@ class Bunch:
current = 0
self._mp_number = int(mp_number)
- self.dtype = np.dtype([('x',np.float),
- ('xp',np.float),
- ('y',np.float),
- ('yp',np.float),
- ('tau',np.float),
- ('delta',np.float)])
+ self.dtype = np.dtype([('x', float),
+ ('xp', float),
+ ('y', float),
+ ('yp', float),
+ ('tau', float),
+ ('delta', float)])
self.particles = np.zeros(self.mp_number, self.dtype)
self.track_alive = track_alive
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index aa4eb68562a3165ba2a3e2f81d19dbcde970670d..69a1c20f28f74503f9a11ab24ccc68a187418086 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -200,8 +200,8 @@ class CavityResonator():
self.detune = detune
self.Vc = Vc
self.theta = theta
- self.beam_phasor = np.zeros(1, dtype=np.complex)
- self.beam_phasor_record = np.zeros((self.ring.h), dtype=np.complex)
+ self.beam_phasor = np.zeros(1, dtype=complex)
+ self.beam_phasor_record = np.zeros((self.ring.h), dtype=complex)
self.tracking = False
self.Vg = 0
self.theta_g = 0
diff --git a/mbtrack2/tracking/wakepotential.py b/mbtrack2/tracking/wakepotential.py
index 14ac15167630cda3e042d5b623fe3e9b80a3532a..19bee813701ab7a7d27299a7f271e2ad23c99f3b 100644
--- a/mbtrack2/tracking/wakepotential.py
+++ b/mbtrack2/tracking/wakepotential.py
@@ -70,7 +70,7 @@ class WakePotential(Element):
"""
- def __init__(self, ring, wakefield, n_bin=65):
+ def __init__(self, ring, wakefield, n_bin=80):
self.wakefield = wakefield
self.types = self.wakefield.wake_components
self.n_types = len(self.wakefield.wake_components)
@@ -793,7 +793,7 @@ class LongRangeResistiveWall(Element):
for wake_type in self.types:
kick = self.get_kick(rank, wake_type)
if wake_type == "Wlong":
- bunch["delta"] += kick / self.ring.E0
+ bunch["delta"] -= kick / self.ring.E0
elif wake_type == "Wxdip":
bunch["xp"] += kick / self.ring.E0
elif wake_type == "Wydip":
diff --git a/mbtrack2/utilities/misc.py b/mbtrack2/utilities/misc.py
index 93332b83e708819e13fa424c3da624997bc90542..2a2bd27b36860ef39f1a6d6bd3c485741ecf13fc 100644
--- a/mbtrack2/utilities/misc.py
+++ b/mbtrack2/utilities/misc.py
@@ -9,13 +9,14 @@ from pathlib import Path
from scipy.interpolate import interp1d
from mbtrack2.impedance.wakefield import Impedance
from mbtrack2.utilities.spectrum import spectral_density
-
-def effective_impedance(ring, imp, m, mu, sigma, M, tuneS, xi=None,
+
+
+def effective_impedance(ring, imp, m, mu, sigma, M, tuneS, xi=None,
mode="Hermite"):
"""
Compute the effective (longitudinal or transverse) impedance.
Formulas from Eq. (1) and (2) p238 of [1].
-
+
Parameters
----------
ring : Synchrotron object
@@ -42,61 +43,112 @@ def effective_impedance(ring, imp, m, mu, sigma, M, tuneS, xi=None,
Zeff : float
effective impedance in [ohm] or in [ohm/m] depanding on the impedance
type.
-
+
References
----------
[1] : Handbook of accelerator physics and engineering, 3rd printing.
"""
-
+
if not isinstance(imp, Impedance):
raise TypeError("{} should be an Impedance object.".format(imp))
-
+
fmin = imp.data.index.min()
fmax = imp.data.index.max()
if fmin > 0:
double_sided_impedance(imp)
-
- if mode == "Hermite":
+
+ if mode in ["Hermite", "Legendre", "Sinusoidal", "Sacherer", "Chebyshev"]:
def h(f):
return spectral_density(frequency=f, sigma=sigma, m=m,
- mode="Hermite")
+ mode=mode)
else:
raise NotImplementedError("Not implemanted yet.")
-
- pmax = fmax/(ring.f0 * M) - 1
- pmin = fmin/(ring.f0 * M) + 1
-
- p = np.arange(pmin,pmax+1)
-
+
if imp.component_type == "long":
+ pmax = fmax/(ring.f0 * M) - 1
+ pmin = fmin/(ring.f0 * M) + 1
+ p = np.arange(pmin, pmax+1)
+
fp = ring.f0*(p*M + mu + m*tuneS)
- fp = fp[np.nonzero(fp)] # Avoid division by 0
- num = np.sum( imp(fp) * h(fp) / (fp*2*np.pi) )
- den = np.sum( h(fp) )
+ fp = fp[np.nonzero(fp)] # Avoid division by 0
+ num = np.sum(imp(fp) * h(fp) / (fp*2*np.pi))
+ den = np.sum(h(fp))
Zeff = num/den
-
+
elif imp.component_type == "xdip" or imp.component_type == "ydip":
if imp.component_type == "xdip":
- tuneXY = ring.tune[0]
- if xi is None :
+ tuneXY = ring.tune[0]-np.floor(ring.tune[0])
+ if xi is None:
xi = ring.chro[0]
elif imp.component_type == "ydip":
- tuneXY = ring.tune[1]
+ tuneXY = ring.tune[1]-np.floor(ring.tune[1])
if xi is None:
xi = ring.chro[1]
+ pmax = fmax/(ring.f0 * M) - 1
+ pmin = fmin/(ring.f0 * M) + 1
+ p = np.arange(pmin, pmax+1)
fp = ring.f0*(p*M + mu + tuneXY + m*tuneS)
f_xi = xi/ring.eta*ring.f0
- num = np.sum( imp(fp) * h(fp - f_xi) )
- den = np.sum( h(fp - f_xi) )
+ num = np.sum(imp(fp) * h(fp - f_xi))
+ den = np.sum(h(fp - f_xi))
Zeff = num/den
else:
raise TypeError("Effective impedance is only defined for long, xdip"
" and ydip impedance type.")
-
+
return Zeff
-def yokoya_elliptic(x_radius , y_radius):
+def head_tail_form_factor(ring, imp, m, sigma, tuneS, xi=None, mode="Hermite", mu=0):
+ M = 1
+ if not isinstance(imp, Impedance):
+ raise TypeError("{} should be an Impedance object.".format(imp))
+
+ fmin = imp.data.index.min()
+ fmax = imp.data.index.max()
+ if fmin > 0:
+ double_sided_impedance(imp)
+
+ if mode in ["Hermite", "Legendre", "Sinusoidal", "Sacherer", "Chebyshev"]:
+ def h(f):
+ return spectral_density(frequency=f, sigma=sigma, m=m,
+ mode=mode)
+ else:
+ raise NotImplementedError("Not implemanted yet.")
+
+ pmax = np.floor(fmax/(ring.f0 * M))
+ pmin = np.ceil(fmin/(ring.f0 * M))
+
+ p = np.arange(pmin, pmax+1)
+
+ if imp.component_type == "long":
+ fp = ring.f0*(p*M + mu + m*tuneS)
+ fp = fp[np.nonzero(fp)] # Avoid division by 0
+ den = np.sum(h(fp))
+
+ elif imp.component_type == "xdip" or imp.component_type == "ydip":
+ if imp.component_type == "xdip":
+ tuneXY = ring.tune[0]-np.floor(ring.tune[0])
+ if xi is None:
+ xi = ring.chro[0]
+ elif imp.component_type == "ydip":
+ tuneXY = ring.tune[1]-np.floor(ring.tune[0])
+ if xi is None:
+ xi = ring.chro[1]
+ fp = ring.f0*(p*M + mu + tuneXY + m*tuneS)
+ f_xi = xi/ring.eta*ring.f0
+ den = np.sum(h(fp - f_xi))
+ else:
+ raise TypeError("Effective impedance is only defined for long, xdip"
+ " and ydip impedance type.")
+
+ return den
+
+
+def tune_shift_from_effective_impedance(Zeff):
+ pass
+
+def yokoya_elliptic(x_radius, y_radius):
"""
Compute Yokoya factors for an elliptic beam pipe.
Function adapted from N. Mounet IW2D.
@@ -127,7 +179,7 @@ def yokoya_elliptic(x_radius , y_radius):
else:
small_semiaxis = x_radius
large_semiaxis = y_radius
-
+
path_to_file = Path(__file__).parent
file = path_to_file / "data" / "Yokoya_elliptic_from_Elias_USPAS.csv"
@@ -140,7 +192,7 @@ def yokoya_elliptic(x_radius , y_radius):
# interpolate Yokoya file at the correct ratio
yoklong = 1
-
+
if y_radius < x_radius:
yokydip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
@@ -150,10 +202,11 @@ def yokoya_elliptic(x_radius , y_radius):
yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
yokydip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
yokxquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
- yokyquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
+ yokyquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)
+
def beam_loss_factor(impedance, frequency, spectrum, ring):
"""
Compute "beam" loss factor using the beam spectrum, uses a sum instead of
@@ -172,7 +225,7 @@ def beam_loss_factor(impedance, frequency, spectrum, ring):
-------
kloss_beam : float
Beam loss factor in [V/C].
-
+
References
----------
[1] : Handbook of accelerator physics and engineering, 3rd printing.
@@ -180,24 +233,25 @@ def beam_loss_factor(impedance, frequency, spectrum, ring):
"""
pmax = np.floor(impedance.data.index.max()/ring.f0)
pmin = np.floor(impedance.data.index.min()/ring.f0)
-
+
if pmin >= 0:
double_sided_impedance(impedance)
pmin = -1*pmax
-
- p = np.arange(pmin+1,pmax)
+
+ p = np.arange(pmin+1, pmax)
pf0 = p*ring.f0
ReZ = np.real(impedance(pf0))
spectral_density = np.abs(spectrum)**2
- # interpolation of the spectrum is needed to avoid problems liked to
+ # interpolation of the spectrum is needed to avoid problems liked to
# division by 0
# computing the spectrum directly to the frequency points gives
# wrong results
spect = interp1d(frequency, spectral_density)
kloss_beam = ring.f0 * np.sum(ReZ*spect(pf0))
-
+
return kloss_beam
+
def double_sided_impedance(impedance):
"""
Add negative frequency points to single sided impedance spectrum following
@@ -209,31 +263,30 @@ def double_sided_impedance(impedance):
Single sided impedance.
"""
fmin = impedance.data.index.min()
-
+
if fmin >= 0:
negative_index = impedance.data.index*-1
negative_data = impedance.data.set_index(negative_index)
-
+
imp_type = impedance.component_type
-
+
if imp_type == "long":
negative_data["imag"] = -1*negative_data["imag"]
-
+
elif (imp_type == "xdip") or (imp_type == "ydip"):
negative_data["real"] = -1*negative_data["real"]
-
+
elif (imp_type == "xquad") or (imp_type == "yquad"):
negative_data["real"] = -1*negative_data["real"]
-
+
else:
raise ValueError("Wrong impedance type")
-
+
try:
negative_data = negative_data.drop(0)
except KeyError:
pass
-
+
all_data = impedance.data.append(negative_data)
all_data = all_data.sort_index()
impedance.data = all_data
-
diff --git a/mbtrack2/utilities/read_impedance.py b/mbtrack2/utilities/read_impedance.py
index 195f0118c376b845cad4cfd363e8b0eab1de2e93..a0550fee0c07f8e4c9e5b07fe27812afde017b2d 100644
--- a/mbtrack2/utilities/read_impedance.py
+++ b/mbtrack2/utilities/read_impedance.py
@@ -66,7 +66,7 @@ def read_CST(file, component_type='long', divide_by=None, imp=True):
return result
-def read_IW2D(file, file_type='Zlong'):
+def read_IW2D(file, file_type='Zlong', output=False):
"""
Read IW2D file format into an Impedance object or a WakeField object.
@@ -76,6 +76,9 @@ def read_IW2D(file, file_type='Zlong'):
Path to the file to read.
file_type : str, optional
Type of the Impedance or WakeField object.
+ output : bool, optional
+ If True, print out the interpolated values.
+ Default is False.
Returns
-------
@@ -99,8 +102,9 @@ def read_IW2D(file, file_type='Zlong'):
if np.any(df.isna()):
index = df.isna().values
df = df.interpolate()
- print("Nan values have been interpolated to:")
- print(df[index])
+ if output:
+ print("Nan values have been interpolated to:")
+ print(df[index])
result = WakeFunction(variable=df.index,
function=df["Wake"],
component_type=file_type[1:])
@@ -109,7 +113,7 @@ def read_IW2D(file, file_type='Zlong'):
return result
-def read_IW2D_folder(folder, suffix, select="WZ"):
+def read_IW2D_folder(folder, suffix, select="WZ", output=False):
"""
Read IW2D results into a WakeField object.
@@ -123,6 +127,9 @@ def read_IW2D_folder(folder, suffix, select="WZ"):
select : str, optional
Select which object to load. "W" for WakeFunction, "Z" for Impedance
and "WZ" or "ZW" for both.
+ output : bool, optional
+ If True, print out the interpolated values.
+ Default is False.
Returns
-------
@@ -148,7 +155,7 @@ def read_IW2D_folder(folder, suffix, select="WZ"):
for key, item in types.items():
for component in components:
name = data_folder / (key + component + suffix)
- res = read_IW2D(file=name, file_type=key + component)
+ res = read_IW2D(file=name, file_type=key + component, output=output)
list_for_wakefield.append(res)
wake = WakeField(list_for_wakefield)
@@ -156,7 +163,7 @@ def read_IW2D_folder(folder, suffix, select="WZ"):
return wake
-def read_ABCI(file, azimuthal=False):
+def read_ABCI(file, azimuthal=False, output=False):
"""
Read ABCI output files [1].
@@ -170,6 +177,9 @@ def read_ABCI(file, azimuthal=False):
If False, it is loaded from the "TRANSVERSE" data. In that case, a -1
factor is applied on the wake to agree with mbtrack2 sign convention.
The default is False.
+ output : bool, optional
+ If True, print out the loaded components header.
+ Default is False.
Returns
-------
@@ -221,8 +231,7 @@ def read_ABCI(file, azimuthal=False):
f' TITLE: {source} WAKE POTENTIAL \n': 'Wxdip',
f' TITLE: REAL PART OF {source} IMPEDANCE \n': 'Zxdip_re',
f' TITLE: IMAGINARY PART OF {source} IMPEDANCE \n': 'Zxdip_im'}
- impedance_type_specification_line = {
- ' DATE: TIME: MROT: NO. OF POINTS:\n'}
+
wake_list = []
start = True
with open(file) as f:
@@ -238,7 +247,8 @@ def read_ABCI(file, azimuthal=False):
# read the header
header = [next(f) for _ in range(4)]
header.insert(0, line)
- print(header)
+ if output:
+ print(header)
# read the body until next TITLE field or end of file
body = []
while True:
@@ -263,8 +273,10 @@ def read_ABCI(file, azimuthal=False):
comp = _read_temp(tmp.name, abci_dict[header[0]])
wake_list.append(comp)
elif abci_dict[header[0]][1:] == "long" and impedance_type == 1:
- comp = _read_temp(tmp.name, abci_dict[header[0]]+'dip')
- wake_list.append(comp)
+ pass
+ # Wlongxdip & Wlongydip not implemented yet
+ # comp = _read_temp(tmp.name, abci_dict[header[0]]+'dip')
+ # wake_list.append(comp)
else:
comp_x = _read_temp(tmp.name, "Wxdip")
comp_y = _read_temp(tmp.name, "Wydip")
@@ -285,8 +297,10 @@ def read_ABCI(file, azimuthal=False):
comp = _read_temp(tmp1.name, "Zlong", tmp2.name)
wake_list.append(comp)
elif abci_dict[header[0]][1:-3] == "long" and impedance_type == 1:
- comp = _read_temp(tmp1.name, "Zlongdip", tmp2.name)
- wake_list.append(comp)
+ pass
+ # Zlongxdip & Zlongydip not implemented yet
+ # comp = _read_temp(tmp1.name, "Zlongdip", tmp2.name)
+ # wake_list.append(comp)
else:
comp_x = _read_temp(tmp1.name, "Zxdip", tmp2.name)
comp_y = _read_temp(tmp1.name, "Zydip", tmp2.name)
diff --git a/mbtrack2/utilities/spectrum.py b/mbtrack2/utilities/spectrum.py
index 91fba8d805e4ad9dce16ca12fd6d4770c9cededa..df78d64110cf84b337ef9b29bb6d043e835512d5 100644
--- a/mbtrack2/utilities/spectrum.py
+++ b/mbtrack2/utilities/spectrum.py
@@ -4,12 +4,14 @@ Module where bunch and beam spectrums and profile are defined.
"""
import numpy as np
+from scipy.special import jv, spherical_jn
-def spectral_density(frequency, sigma, m = 1, mode="Hermite"):
+
+def spectral_density(frequency, sigma, m=1, k=0, mode="Hermite"):
"""
Compute the spectral density of different modes for various values of the
head-tail mode number, based on Table 1 p238 of [1].
-
+
Parameters
----------
frequency : list or numpy array
@@ -18,26 +20,43 @@ def spectral_density(frequency, sigma, m = 1, mode="Hermite"):
RMS bunch length in [s]
m : int, optional
head-tail (or azimutal/synchrotron) mode number
+ k : int, optional
+ radial mode number (such that |q|=m+2k, where |q| is the head-tail mode number)
mode: str, optional
type of the mode taken into account for the computation:
- -"Hermite" modes for Gaussian bunches
+ -"Hermite" modes for Gaussian bunches (typical for electrons)
+ -"Chebyshev" for airbag bunches
+ -"Legendre" for parabolic bunches (typical for protons)
+ -"Sacherer" or "Sinusoidal" simplifying approximation of Legendre modes from [3]
Returns
-------
numpy array
-
+
References
----------
[1] : Handbook of accelerator physics and engineering, 3rd printing.
+ [2] : Ng, K. Y. (2005). Physics of intensity dependent beam instabilities. WORLD SCIENTIFIC. https://doi.org/10.1142/5835
+ [3] : Sacherer, F. J. (1972). Methods for computing bunched beam instabilities. CERN Tech. rep. CERN/SI-BR/72-5 https://cds.cern.ch/record/2291670?ln=en
"""
-
+
if mode == "Hermite":
- return (2*np.pi*frequency*sigma)**(2*m)*np.exp(
- -1*(2*np.pi*frequency*sigma)**2)
+ return 1/(np.math.factorial(m)*2**m)*(2*np.pi*frequency*sigma)**(2*m)*np.exp(
+ -(2*np.pi*frequency*sigma)**2)
+ elif mode == "Chebyshev":
+ tau_l = 4*sigma
+ return (jv(m, 2*np.pi*frequency*tau_l))**2
+ elif mode == "Legendre":
+ tau_l = 4*sigma
+ return (spherical_jn(m, np.abs(2*np.pi*frequency*tau_l)))**2
+ elif mode == "Sacherer" or mode == "Sinusoidal":
+ y = 4*2*np.pi*frequency*sigma/np.pi
+ return (2*(m+1)/np.pi*1/np.abs(y**2-(m+1)**2)*np.sqrt(1+(-1)**m*np.cos(np.pi*y)))**2
else:
raise NotImplementedError("Not implemanted yet.")
-
-def gaussian_bunch_spectrum(frequency, sigma):
+
+
+def gaussian_bunch_spectrum(frequency, sigma):
"""
Compute a Gaussian bunch spectrum [1].
@@ -52,7 +71,7 @@ def gaussian_bunch_spectrum(frequency, sigma):
-------
bunch_spectrum : array
Bunch spectrum sampled at points given in frequency.
-
+
References
----------
[1] : Gamelin, A. (2018). Collective effects in a transient microbunching
@@ -60,7 +79,8 @@ def gaussian_bunch_spectrum(frequency, sigma):
"""
return np.exp(-1/2*(2*np.pi*frequency)**2*sigma**2)
-def gaussian_bunch(time, sigma):
+
+def gaussian_bunch(time, sigma):
"""
Compute a Gaussian bunch profile.
@@ -77,10 +97,10 @@ def gaussian_bunch(time, sigma):
Bunch profile in [s**-1] sampled at points given in time.
"""
return np.exp(-1/2*(time**2/sigma**2))/(sigma*np.sqrt(2*np.pi))
-
-
-def beam_spectrum(frequency, M, bunch_spacing, sigma=None,
- bunch_spectrum=None):
+
+
+def beam_spectrum(frequency, M, bunch_spacing, sigma=None,
+ bunch_spectrum=None):
"""
Compute the beam spectrum assuming constant spacing between bunches [1].
@@ -107,13 +127,13 @@ def beam_spectrum(frequency, M, bunch_spacing, sigma=None,
[1] Rumolo, G - Beam Instabilities - CAS - CERN Accelerator School:
Advanced Accelerator Physics Course - 2014, Eq. 9
"""
-
+
if bunch_spectrum is None:
bunch_spectrum = gaussian_bunch_spectrum(frequency, sigma)
-
- beam_spectrum = (bunch_spectrum * np.exp(1j*np.pi*frequency *
- bunch_spacing*(M-1)) *
- np.sin(M*np.pi*frequency*bunch_spacing) /
+
+ beam_spectrum = (bunch_spectrum * np.exp(1j*np.pi*frequency *
+ bunch_spacing*(M-1)) *
+ np.sin(M*np.pi*frequency*bunch_spacing) /
np.sin(np.pi*frequency*bunch_spacing))
-
- return beam_spectrum
\ No newline at end of file
+
+ return beam_spectrum