From 29b70450631d9ae27c80d5552cbb49ed82dfb7c4 Mon Sep 17 00:00:00 2001
From: Gamelin Alexis <gamelin@synchrotron-soleil.fr>
Date: Thu, 12 Nov 2020 11:38:15 +0100
Subject: [PATCH] Add WakeFunction to CircularResisitiveWall and a Coating
element
The exact wake function for the transverse wake function is wrong (from paper) ! Need to find a better expression.
---
collective_effects/__init__.py | 2 +-
collective_effects/resistive_wall.py | 355 +++++++++++++++++++++++++--
2 files changed, 331 insertions(+), 26 deletions(-)
diff --git a/collective_effects/__init__.py b/collective_effects/__init__.py
index bdea55f..a2b820a 100644
--- a/collective_effects/__init__.py
+++ b/collective_effects/__init__.py
@@ -6,7 +6,7 @@ Created on Tue Jan 14 12:25:30 2020
"""
from mbtrack2.collective_effects.instabilities import mbi_threshold, cbi_threshold
-from mbtrack2.collective_effects.resistive_wall import skin_depth, CircularResistiveWall
+from mbtrack2.collective_effects.resistive_wall import (skin_depth, CircularResistiveWall, Coating)
from mbtrack2.collective_effects.resonator import Resonator, PureInductive, PureResistive
from mbtrack2.collective_effects.tapers import StupakovRectangularTaper, StupakovCircularTaper
from mbtrack2.collective_effects.wakefield import ComplexData, Impedance, WakeFunction, WakeField
diff --git a/collective_effects/resistive_wall.py b/collective_effects/resistive_wall.py
index 2df32f4..922cc94 100644
--- a/collective_effects/resistive_wall.py
+++ b/collective_effects/resistive_wall.py
@@ -8,58 +8,363 @@ Define resistive wall elements based on the WakeField class.
import numpy as np
from scipy.constants import mu_0, epsilon_0, c
-from mbtrack2.collective_effects.wakefield import WakeField, Impedance
+from scipy.integrate import quad
+from mbtrack2.collective_effects.wakefield import WakeField, Impedance, WakeFunction
-def skin_depth(omega, rho, mu_r = 1, epsilon_r = 1):
- """ General formula for the skin depth
+def skin_depth(frequency, rho, mu_r = 1, epsilon_r = 1):
+ """
+ General formula for the skin depth.
Parameters
----------
- omega: angular frequency in [Hz]
- rho: resistivity in [ohm.m]
- mu_r : relative magnetic permeability, optional
- epsilon_r : relative electric permittivity, optional
+ frequency : array of float
+ Frequency points in [Hz].
+ rho : float
+ Resistivity in [ohm.m].
+ mu_r : float, optional
+ Relative magnetic permeability.
+ epsilon_r : float, optional
+ Relative electric permittivity.
Returns
-------
- delta : skin depth in [m]
+ delta : array of float
+ Skin depth in [m].
"""
- delta = np.sqrt(2*rho/(np.abs(omega)*mu_r*mu_0))*np.sqrt(np.sqrt(1 +
- (rho*np.abs(omega)*epsilon_r*epsilon_0)**2 )
- + rho*np.abs(omega)*epsilon_r*epsilon_0)
+ delta = (np.sqrt(2*rho/(np.abs(2*np.pi*frequency)*mu_r*mu_0)) *
+ np.sqrt(np.sqrt(1 + (rho*np.abs(2*np.pi*frequency) *
+ epsilon_r*epsilon_0)**2 )
+ + rho*np.abs(2*np.pi*frequency)*epsilon_r*epsilon_0))
return delta
class CircularResistiveWall(WakeField):
"""
- Resistive wall WakeField element for a circular beam pipe, approximated
- formulas from Eq. (2.77) of Chao book [1].
+ Resistive wall WakeField element for a circular beam pipe.
+
+ Impedance from approximated formulas from Eq. (2.77) of Chao book [1].
+ Wake function formulas from [2].
+
+ !!! The exact formula for the transverse wake function is wrong !!!
Parameters
----------
- frequency: frequency points where the impedance will be evaluated in [Hz]
- length: beam pipe length in [m]
- rho: resistivity in [ohm.m]
- radius: beam pipe radius in [m]
- mu_r : relative magnetic permeability, optional
- epsilon_r : relative electric permittivity, optional
+ 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].
+ length : float
+ Beam pipe length in [m].
+ rho : float
+ Resistivity in [ohm.m].
+ radius : float
+ Beam pipe radius in [m].
+ exact : bool, optional
+ If False, approxmiated formulas are used for the wake function
+ computations.
+
+ References
+ ----------
+ [1] : Chao, A. W. (1993). Physics of collective beam instabilities in high
+ energy accelerators. Wiley.
+ [2] : Skripka, Galina, et al. "Simultaneous computation of intrabunch and
+ interbunch collective beam motions in storage rings." Nuclear Instruments
+ and Methods in Physics Research Section A: Accelerators, Spectrometers,
+ Detectors and Associated Equipment 806 (2016): 221-230.
+
"""
- def __init__(self, frequency, length, rho, radius, mu_r = 1, epsilon_r = 1):
+ def __init__(self, time, frequency, length, rho, radius, exact=False):
super().__init__()
+ self.length = length
+ self.rho = rho
+ self.radius = radius
+
omega = 2*np.pi*frequency
- Z1 = length*(1 + np.sign(omega)*1j)*rho/(
- 2*np.pi*radius*skin_depth(omega,rho))
- Z2 = c/omega*length*(1 + np.sign(omega)*1j)*rho/(
- np.pi*radius**3*skin_depth(omega,rho))
+ Z1 = length*(1 + np.sign(frequency)*1j)*rho/(
+ 2*np.pi*radius*skin_depth(frequency,rho))
+ Z2 = c/omega*length*(1 + np.sign(frequency)*1j)*rho/(
+ np.pi*radius**3*skin_depth(frequency,rho))
+
+ Wl = self.LongitudinalWakeFunction(time, exact)
+ Wt = self.TransverseWakeFunction(time, exact)
Zlong = Impedance(variable = frequency, function = Z1, impedance_type='long')
Zxdip = Impedance(variable = frequency, function = Z2, impedance_type='xdip')
Zydip = Impedance(variable = frequency, function = Z2, impedance_type='ydip')
+ Wlong = WakeFunction(variable = time, function = Wl, wake_type="long")
+ Wxdip = WakeFunction(variable = time, function = Wt, wake_type="xdip")
+ Wydip = WakeFunction(variable = time, function = Wt, wake_type="ydip")
+
+ super().append_to_model(Zlong)
+ super().append_to_model(Zxdip)
+ super().append_to_model(Zydip)
+ super().append_to_model(Wlong)
+ super().append_to_model(Wxdip)
+ super().append_to_model(Wydip)
+
+ def LongitudinalWakeFunction(self, time, exact=False):
+ """
+ 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.
+
+ Parameters
+ ----------
+ time : array of float
+ Time points where the wake function is evaluated in [s].
+ exact : bool, optional
+ If True, the exact expression is used. The default is False.
+
+ Returns
+ -------
+ wl : array of float
+ Longitudinal wake function in [V/C].
+
+ References
+ ----------
+ [1] : Skripka, Galina, et al. "Simultaneous computation of intrabunch and
+ interbunch collective beam motions in storage rings." Nuclear Instruments
+ and Methods in Physics Research Section A: Accelerators, Spectrometers,
+ Detectors and Associated Equipment 806 (2016): 221-230.
+ """
+
+ Z0 = mu_0*c
+
+ 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 )
+ else:
+ wl = (-1/(4*np.pi * self.radius)
+ * np.sqrt(Z0 * self.rho / (c * np.pi) )
+ / time**(3/2) )
+ ind = np.isnan(wl)
+ wl[ind] = 0
+
+ return wl
+
+ def TransverseWakeFunction(self, time, exact=False):
+ """
+ Compute the transverse wake function of a circular resistive wall
+ using Eq. (25), or approxmiated expression Eq. (26), of [1]. The
+ approxmiated expression is valid if the time is large compared to the
+ characteristic time t0.
+
+ Parameters
+ ----------
+ time : array of float
+ Time points where the wake function is evaluated in [s].
+ exact : bool, optional
+ If True, the exact expression is used. The default is False.
+
+ Returns
+ -------
+ wt : array of float
+ Transverse wake function in [V/C].
+
+ References
+ ----------
+ [1] : Skripka, Galina, et al. "Simultaneous computation of intrabunch and
+ interbunch collective beam motions in storage rings." Nuclear Instruments
+ and Methods in Physics Research Section A: Accelerators, Spectrometers,
+ Detectors and Associated Equipment 806 (2016): 221-230.
+ """
+
+ Z0 = mu_0*c
+
+ if exact==True:
+ self.t0 = (2*self.rho*self.radius**2 / Z0)**(1/3) / c
+ factor = 8*Z0*c/(np.pi * self.radius**4) * self.length * -1
+ 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) *
+ 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 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):
+ return ( (-1 * np.exp(-1* (x**2) * t / self.t0) ) / (x**6 + 8) )
+
+class Coating(WakeField):
+
+ def __init__(self, frequency, length, rho1, rho2, radius, thickness, approx=False):
+ """
+ WakeField element for a coated circular beam pipe.
+
+ The longitudinal and tranverse impedances are computed using formulas
+ from [1].
+
+ Parameters
+ ----------
+ f : array of float
+ Frequency points where the impedance is evaluated in [Hz].
+ length : float
+ Length of the beam pipe to consider in [m].
+ rho1 : float
+ Resistivity of the coating in [ohm.m].
+ rho2 : float
+ Resistivity of the bulk material in [ohm.m].
+ radius : float
+ Radius of the beam pipe to consier in [m].
+ thickness : float
+ Thickness of the coating in [m].
+ approx : bool, optional
+ If True, used approxmiated formula. The default is False.
+
+ References
+ ----------
+ [1] : Migliorati, M., E. Belli, and M. Zobov. "Impact of the resistive
+ wall impedance on beam dynamics in the Future Circular e+ e− Collider."
+ Physical Review Accelerators and Beams 21.4 (2018): 041001.
+
+ """
+ super().__init__()
+
+ self.length = length
+ self.rho1 = rho1
+ self.rho2 = rho2
+ self.radius = radius
+ self.thickness = thickness
+
+ Zl = self.LongitudinalImpedance(frequency, approx)
+ Zt = self.TransverseImpedance(frequency, approx)
+
+ Zlong = Impedance(variable = frequency, function = Zl, impedance_type='long')
+ Zxdip = Impedance(variable = frequency, function = Zt, impedance_type='xdip')
+ Zydip = Impedance(variable = frequency, function = Zt, impedance_type='ydip')
super().append_to_model(Zlong)
super().append_to_model(Zxdip)
- super().append_to_model(Zydip)
\ No newline at end of file
+ super().append_to_model(Zydip)
+
+ def LongitudinalImpedance(self, f, approx):
+ """
+ Compute the longitudinal impedance of a coating using Eq. (5), or
+ approxmiated expression Eq. (8), of [1]. The approxmiated expression
+ is valid if the skin depth of the coating is large compared to the
+ coating thickness.
+
+ Parameters
+ ----------
+ f : array of float
+ Frequency points where the impedance is evaluated in [Hz].
+ approx : bool
+ If True, used approxmiated formula.
+
+ Returns
+ -------
+ Zl : array
+ Longitudinal impedance values in [ohm].
+
+ References
+ ----------
+ [1] : Migliorati, M., E. Belli, and M. Zobov. "Impact of the resistive
+ wall impedance on beam dynamics in the Future Circular e+ e− Collider."
+ Physical Review Accelerators and Beams 21.4 (2018): 041001.
+
+ """
+
+ Z0 = mu_0*c
+ factor = Z0*f/(2*c*self.radius)*self.length
+ skin1 = skin_depth(f, self.rho1)
+ skin2 = skin_depth(f, self.rho2)
+
+ if approx == False:
+ alpha = skin1/skin2
+ tanh = np.tanh( (1 - 1j*np.sign(f)) * self.thickness / skin1 )
+ bracket = ( (np.sign(f) - 1j) * skin1 *
+ (alpha * tanh + 1) / (alpha + tanh) )
+ else:
+ valid_approx = self.thickness / skin1
+ if valid_approx < 0.01:
+ print("Approximation is not valid. Returning impedance anyway.")
+ bracket = ( (np.sign(f) - 1j) * skin2 - 2 * 1j * self.thickness *
+ (1 - self.rho2/self.rho1) )
+
+ Zl = factor * bracket
+
+ return Zl
+
+ def TransverseImpedance(self, f, approx):
+ """
+ Compute the transverse impedance of a coating using Eq. (6), or
+ approxmiated expression Eq. (9), of [1]. The approxmiated expression
+ is valid if the skin depth of the coating is large compared to the
+ coating thickness.
+
+ Parameters
+ ----------
+ f : array of float
+ Frequency points where the impedance is evaluated in [Hz].
+ approx : bool
+ If True, used approxmiated formula.
+
+ Returns
+ -------
+ Zt : array
+ Transverse impedance values in [ohm].
+
+ References
+ ----------
+ [1] : Migliorati, M., E. Belli, and M. Zobov. "Impact of the resistive
+ wall impedance on beam dynamics in the Future Circular e+ e− Collider."
+ Physical Review Accelerators and Beams 21.4 (2018): 041001.
+
+ """
+
+ Z0 = mu_0*c
+ factor = Z0/(2*np.pi*self.radius**3)*self.length
+ skin1 = skin_depth(f, self.rho1)
+ skin2 = skin_depth(f, self.rho2)
+
+ if approx == False:
+ alpha = skin1/skin2
+ tanh = np.tanh( (1 - 1j*np.sign(f)) * self.thickness / skin1 )
+ bracket = ( (1 - 1j*np.sign(f)) * skin1 *
+ (alpha * tanh + 1) / (alpha + tanh) )
+ else:
+ valid_approx = self.thickness / skin1
+ if valid_approx < 0.01:
+ print("Approximation is not valid. Returning impedance anyway.")
+ bracket = ( (1 - 1j*np.sign(f)) * skin2 - 2 * 1j * self.thickness
+ * np.sign(f) * (1 - self.rho2/self.rho1) )
+
+ Zt = factor * bracket
+
+ return Zt
+
+
+
+
+
+
+
\ No newline at end of file
--
GitLab