# -*- coding: utf-8 -*-
"""
Define resistive wall elements based on the WakeField class.

@author: Alexis Gamelin
@date: 14/01/2020
"""

import numpy as np
from scipy.constants import mu_0, epsilon_0, c
from scipy.integrate import quad
from mbtrack2.collective_effects.wakefield import WakeField, Impedance, WakeFunction

def skin_depth(frequency, rho, mu_r = 1, epsilon_r = 1):
    """
    General formula for the skin depth.
    
    Parameters
    ----------
    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 : array of float
        Skin depth in [m].
    
    """
    
    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.
    
    Impedance from approximated formulas from Eq. (2.77) of Chao book [1].
    Wake function formulas from [2].
        
    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].
    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, 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(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) ) * self.length
            ind = np.isnan(wl)
            wl[ind] = 0
    
        return -1*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.
        
        Exact expression (Eq. (25) from [1]) is corrected by factor (c * 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 = -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) * 
                                      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 ( (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)
        
    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