diff --git a/collective_effects/resistive_wall.py b/collective_effects/resistive_wall.py
index 1bb385419d6324cd274945b6194f4cb3b0d513f2..c0ae255ff117f11cd76002aaf6b9a0cff6d73c3d 100644
--- a/collective_effects/resistive_wall.py
+++ b/collective_effects/resistive_wall.py
@@ -62,6 +62,10 @@ class CircularResistiveWall(WakeField):
exact : bool, optional
If False, approxmiated formulas are used for the wake function
computations.
+ atol : float, optional
+ Absolute tolerance used to enforce fundamental theorem of beam loading
+ for the exact expression of the longitudinal wake function.
+ Default is 1e-20.
References
----------
@@ -74,12 +78,15 @@ class CircularResistiveWall(WakeField):
"""
- def __init__(self, time, frequency, length, rho, radius, exact=False):
+ def __init__(self, time, frequency, length, rho, radius, exact=False,
+ atol=1e-20):
super().__init__()
self.length = length
self.rho = rho
self.radius = radius
+ self.Z0 = mu_0*c
+ self.t0 = (2*self.rho*self.radius**2 / self.Z0)**(1/3) / c
omega = 2*np.pi*frequency
Z1 = length*(1 + np.sign(frequency)*1j)*rho/(
@@ -87,7 +94,7 @@ class CircularResistiveWall(WakeField):
Z2 = c/omega*length*(1 + np.sign(frequency)*1j)*rho/(
np.pi*radius**3*skin_depth(frequency,rho))
- Wl = self.LongitudinalWakeFunction(time, exact)
+ Wl = self.LongitudinalWakeFunction(time, exact, atol)
Wt = self.TransverseWakeFunction(time, exact)
Zlong = Impedance(variable = frequency, function = Z1, impedance_type='long')
@@ -104,12 +111,15 @@ class CircularResistiveWall(WakeField):
super().append_to_model(Wxdip)
super().append_to_model(Wydip)
- def LongitudinalWakeFunction(self, time, exact=False):
+ def LongitudinalWakeFunction(self, time, exact=False, atol=1e-20):
"""
Compute the longitudinal wake function of a circular resistive wall
using Eq. (22), or approxmiated expression Eq. (24), of [1]. The
approxmiated expression is valid if the time is large compared to the
characteristic time t0.
+
+ If some time value is smaller than atol, then the fundamental theorem
+ of beam loading is applied: Wl(0) = Wl(0+)/2.
Parameters
----------
@@ -117,6 +127,10 @@ class CircularResistiveWall(WakeField):
Time points where the wake function is evaluated in [s].
exact : bool, optional
If True, the exact expression is used. The default is False.
+ atol : float, optional
+ Absolute tolerance used to enforce fundamental theorem of beam loading
+ for the exact expression of the longitudinal wake function.
+ Default is 1e-20.
Returns
-------
@@ -130,30 +144,17 @@ class CircularResistiveWall(WakeField):
and Methods in Physics Research Section A: Accelerators, Spectrometers,
Detectors and Associated Equipment 806 (2016): 221-230.
"""
-
- Z0 = mu_0*c
-
+ wl = np.zeros_like(time)
+ idx1 = time < 0
+ wl[idx1] = 0
if exact==True:
- self.t0 = (2*self.rho*self.radius**2 / Z0)**(1/3) / c
- factor = 4*Z0*c/(np.pi * self.radius**2) * self.length * -1
- wl = np.zeros_like(time)
-
- for i, t in enumerate(time):
- val, err = quad(lambda z:self.function(t, z), 0, np.inf)
- if t < 0:
- wl[i] = 0
- else:
- wl[i] = factor * ( np.exp(-t/self.t0) / 3 *
- np.cos( np.sqrt(3) * t / self.t0 )
- - np.sqrt(2) / np.pi * val )
+ idx2 = time > 20 * self.t0
+ idx3 = np.logical_not(np.logical_or(idx1,idx2))
+ wl[idx3] = self.__LongWakeExact(time[idx3], atol)
else:
- wl = (1/(4*np.pi * self.radius)
- * np.sqrt(Z0 * self.rho / (c * np.pi) )
- / time**(3/2) ) * self.length
- ind = np.isnan(wl)
- wl[ind] = 0
-
- return -1*wl
+ idx2 = np.logical_not(idx1)
+ wl[idx2] = self.__LongWakeApprox(time[idx2])
+ return wl
def TransverseWakeFunction(self, time, exact=False):
"""
@@ -183,35 +184,59 @@ class CircularResistiveWall(WakeField):
and Methods in Physics Research Section A: Accelerators, Spectrometers,
Detectors and Associated Equipment 806 (2016): 221-230.
"""
-
- Z0 = mu_0*c
-
+ wt = np.zeros_like(time)
+ idx1 = time < 0
+ wt[idx1] = 0
if exact==True:
- self.t0 = (2*self.rho*self.radius**2 / Z0)**(1/3) / c
- factor = -1 * (8 * Z0 * c**2 * self.t0) / (np.pi * self.radius**4) * self.length
- wt = np.zeros_like(time)
-
- for i, t in enumerate(time):
- val, err = quad(lambda z:self.function2(t, z), 0, np.inf)
- if t < 0:
- wt[i] = 0
- else:
- wt[i] = factor * ( 1 / 12 * (-1 * np.exp(-t/self.t0) *
+ idx2 = time > 20 * self.t0
+ idx3 = np.logical_not(np.logical_or(idx1,idx2))
+ wt[idx3] = self.__TransWakeExact(time[idx3])
+ else:
+ idx2 = np.logical_not(idx1)
+ wt[idx2] = self.__TransWakeApprox(time[idx2])
+ return wt
+
+ def __LongWakeExact(self, time, atol):
+ wl = np.zeros_like(time)
+ factor = 4*self.Z0*c/(np.pi * self.radius**2) * self.length
+ for i, t in enumerate(time):
+ val, err = quad(lambda z:self.__function(t, z), 0, np.inf)
+ wl[i] = factor * ( np.exp(-t/self.t0) / 3 *
+ np.cos( np.sqrt(3) * t / self.t0 )
+ - np.sqrt(2) / np.pi * val )
+ if np.isclose(0, t, atol=atol):
+ wl[i] = wl[i]/2
+ return wl
+
+ def __TransWakeExact(self, time):
+ wt = np.zeros_like(time)
+ factor = ((8 * self.Z0 * c**2 * self.t0) / (np.pi * self.radius**4) *
+ self.length)
+ for i, t in enumerate(time):
+ val, err = quad(lambda z:self.__function2(t, z), 0, np.inf)
+ wt[i] = factor * ( 1 / 12 * (-1 * np.exp(-t/self.t0) *
np.cos( np.sqrt(3) * t / self.t0 ) +
np.sqrt(3) * np.exp(-t/self.t0) *
np.sin( np.sqrt(3) * t / self.t0 ) ) -
np.sqrt(2) / np.pi * val )
- else:
- wt = (1 / (np.pi * self.radius**3) * np.sqrt(Z0 * c * self.rho / np.pi)
- / time**(1/2) * self.length * -1)
- ind = np.isnan(wt)
- wt[ind] = 0
- return -1*wt
-
- def function(self, t, x):
+ return wt
+
+ def __LongWakeApprox(self, t):
+ wl = - 1 * ( 1 / (4*np.pi * self.radius) *
+ np.sqrt(self.Z0 * self.rho / (c * np.pi) ) /
+ t ** (3/2) ) * self.length
+ return wl
+
+ def __TransWakeApprox(self, t):
+ wt = (1 / (np.pi * self.radius**3) *
+ np.sqrt(self.Z0 * c * self.rho / np.pi)
+ / t ** (1/2) * self.length)
+ return wt
+
+ def __function(self, t, x):
return ( (x**2 * np.exp(-1* (x**2) * t / self.t0) ) / (x**6 + 8) )
- def function2(self, t, x):
+ def __function2(self, t, x):
return ( (-1 * np.exp(-1* (x**2) * t / self.t0) ) / (x**6 + 8) )
class Coating(WakeField):
diff --git a/collective_effects/resonator.py b/collective_effects/resonator.py
index f11948ee3abe5b923d4d5f66bbcf2e483b64df09..41168b35793cc2b6e72b68ee4bae108ce3b17773 100644
--- a/collective_effects/resonator.py
+++ b/collective_effects/resonator.py
@@ -12,27 +12,29 @@ from mbtrack2.collective_effects.wakefield import (WakeField, Impedance,
WakeFunction)
class Resonator(WakeField):
- def __init__(self, Rs, fr, Q, plane, n_wake=1e6, n_imp=1e6, imp_freq_lim=100e9):
+ def __init__(self, time, frequency, Rs, fr, Q, plane, atol=1e-20):
"""
Resonator model WakeField element which computes the impedance and the
wake function in both longitudinal and transverse case.
Parameters
----------
+ time : array of float
+ Time points where the wake function will be evaluated in [s].
+ frequency : array of float
+ Frequency points where the impedance will be evaluated in [Hz].
Rs : float
Shunt impedance in [ohm].
fr : float
Resonance frequency in [Hz].
Q : float
Quality factor.
- plane : str
+ plane : str or list
Plane on which the resonator is used: "long", "x" or "y".
- n_wake : int or float, optional
- Number of points used in the wake function.
- n_imp : int or float, optional
- Number of points used in the impedance.
- imp_freq_lim : float, optional
- Maximum frequency used in the impedance.
+ atol : float, optional
+ Absolute tolerance used to enforce fundamental theorem of beam
+ loading for the exact expression of the longitudinal wake function.
+ Default is 1e-20.
References
----------
@@ -45,62 +47,54 @@ class Resonator(WakeField):
self.fr = fr
self.wr = 2 * np.pi * self.fr
self.Q = Q
- self.n_wake = int(n_wake)
- self.n_imp = int(n_imp)
- self.imp_freq_lim = imp_freq_lim
- self.plane = plane
-
- self.timestop = round(np.log(1000)/self.wr*2*self.Q, 12)
-
+ if isinstance(plane, str):
+ self.plane = [plane]
+ elif isinstance(plane, list):
+ self.plane = plane
+
if self.Q >= 0.5:
self.Q_p = np.sqrt(self.Q**2 - 0.25)
else:
self.Q_p = np.sqrt(0.25 - self.Q**2)
-
self.wr_p = (self.wr*self.Q_p)/self.Q
- if self.plane == "long":
-
- freq = np.linspace(start=1, stop=self.imp_freq_lim, num=self.n_imp)
- imp = Impedance(variable=freq,
- function=self.long_impedance(freq),
- impedance_type="long")
- super().append_to_model(imp)
-
- time = np.linspace(start=0, stop=self.timestop, num=self.n_wake)
- wake = WakeFunction(variable=time,
- function=self.long_wake_function(time),
- wake_type="long")
- super().append_to_model(wake)
-
- elif self.plane == "x" or self.plane == "y" :
-
- freq = np.linspace(start=1, stop=self.imp_freq_lim, num=self.n_imp)
- imp = Impedance(variable=freq,
- function=self.transverse_impedance(freq),
- impedance_type=self.plane + "dip")
- super().append_to_model(imp)
-
- time = np.linspace(start=0, stop=self.timestop, num=self.n_wake)
- wake = WakeFunction(variable=time,
- function=self.transverse_wake_function(time),
- wake_type=self.plane + "dip")
- super().append_to_model(wake)
- else:
- raise ValueError("Plane must be: long, x or y")
+ for dim in self.plane:
+ if dim == "long":
+ Zlong = Impedance(variable=frequency,
+ function=self.long_impedance(frequency),
+ impedance_type="long")
+ super().append_to_model(Zlong)
+ Wlong = WakeFunction(variable=time,
+ function=self.long_wake_function(time, atol),
+ wake_type="long")
+ super().append_to_model(Wlong)
+
+ elif dim == "x" or dim == "y":
+ Zdip = Impedance(variable=frequency,
+ function=self.transverse_impedance(frequency),
+ impedance_type=dim + "dip")
+ super().append_to_model(Zdip)
+ Wdip = WakeFunction(variable=time,
+ function=self.transverse_wake_function(time),
+ wake_type=dim + "dip")
+ super().append_to_model(Wdip)
+ else:
+ raise ValueError("Plane must be: long, x or y")
- def long_wake_function(self, t):
+ def long_wake_function(self, t, atol):
if self.Q >= 0.5:
- return ( (self.wr * self.Rs / self.Q) *
+ wl = ( (self.wr * self.Rs / self.Q) *
np.exp(-1* self.wr * t / (2 * self.Q) ) *
(np.cos(self.wr_p * t) -
np.sin(self.wr_p * t) / (2 * self.Q_p) ) )
-
elif self.Q < 0.5:
- return ( (self.wr * self.Rs / self.Q) *
+ wl = ( (self.wr * self.Rs / self.Q) *
np.exp(-1* self.wr * t / (2 * self.Q) ) *
(np.cosh(self.wr_p * t) -
np.sinh(self.wr_p * t) / (2 * self.Q_p) ) )
+ if np.any(np.abs(t) < atol):
+ wl[np.abs(t) < atol] = wl[np.abs(t) < atol]/2
+ return wl
def long_impedance(self, f):
return self.Rs / (1 + 1j * self.Q * (f/self.fr - self.fr/f))
diff --git a/collective_effects/wakefield.py b/collective_effects/wakefield.py
index 36b37f2e8d1060aa13ad6a8b9ec1363128cd202f..3af5ea02311e56c09036ec9d601d0dedc4fc11b8 100644
--- a/collective_effects/wakefield.py
+++ b/collective_effects/wakefield.py
@@ -691,14 +691,16 @@ class WakeField:
"""
Return an array of the impedance component names for the element.
"""
- return np.array([comp for comp in dir(self) if re.match(r'[Z]', comp)])
+ valid = ["Zlong", "Zxdip", "Zydip", "Zxquad", "Zyquad"]
+ return np.array([comp for comp in dir(self) if comp in valid])
@property
def wake_components(self):
"""
Return an array of the wake function component names for the element.
"""
- return np.array([comp for comp in dir(self) if re.match(r'[W]', comp)])
+ valid = ["Wlong", "Wxdip", "Wydip", "Wxquad", "Wyquad"]
+ return np.array([comp for comp in dir(self) if comp in valid])
@staticmethod
def add_wakefields(wake1, beta1, wake2, beta2):