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

@@ -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

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)
+        super().append_to_model(Zydip)
\ No newline at end of file
+        
+    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