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

Merge branch 'RW_LR' into Prepare_version

parents c5a79d24 f3296289
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
collective_effects/resistive_wall.py
+ 72
− 47
View file @ f0f97c3a
......@@ -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):
......
collective_effects/resonator.py
+ 43
− 49
View file @ f0f97c3a
......@@ -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))
......
collective_effects/wakefield.py
+ 4
− 2
View file @ f0f97c3a
......@@ -690,14 +690,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):
......
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