diff --git a/collective_effects/__init__.py b/collective_effects/__init__.py
index 67c3d4f6011d9fc22f490119b06933d996d50655..d14360ffa3ea6025eeba8dd2a0b1cead9fb592f3 100644
--- a/collective_effects/__init__.py
+++ b/collective_effects/__init__.py
@@ -7,9 +7,10 @@ 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.resonator import Resonator
 from mbtrack2.collective_effects.tapers import StupakovRectangularTaper, StupakovCircularTaper
-from mbtrack2.collective_effects.wakefield import ComplexData, Impedance, WakeFunction, ImpedanceModel, Wakefield
-from mbtrack2.collective_effects.wakefield import read_CST, read_IW2D, spectral_density
-from mbtrack2.collective_effects.wakefield import gaussian_bunch_spectrum, beam_spectrum, effective_impedance
-from mbtrack2.collective_effects.wakefield import yokoya_elliptic, beam_loss_factor
-from mbtrack2.collective_effects.wakefield import double_sided_impedance
\ No newline at end of file
+from mbtrack2.collective_effects.wakefield import ComplexData, Impedance, WakeFunction, ImpedanceModel, WakeField
+from mbtrack2.collective_effects.utilities import read_CST, read_IW2D, spectral_density
+from mbtrack2.collective_effects.utilities import gaussian_bunch_spectrum, beam_spectrum, effective_impedance
+from mbtrack2.collective_effects.utilities import yokoya_elliptic, beam_loss_factor
+from mbtrack2.collective_effects.utilities import double_sided_impedance
\ No newline at end of file
diff --git a/collective_effects/resistive_wall.py b/collective_effects/resistive_wall.py
index 240e30e268dfa956dd1a2dbc005e88cfe6897dca..2df32f47833f7caba38c24dc733523e1097e44be 100644
--- a/collective_effects/resistive_wall.py
+++ b/collective_effects/resistive_wall.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 """
-Define resistive wall elements based on the Wakefield class.
+Define resistive wall elements based on the WakeField class.
 
 @author: Alexis Gamelin
 @date: 14/01/2020
@@ -8,7 +8,7 @@ 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 mbtrack2.collective_effects.wakefield import WakeField, Impedance
 
 def skin_depth(omega, rho, mu_r = 1, epsilon_r = 1):
     """ General formula for the skin depth
@@ -32,9 +32,9 @@ def skin_depth(omega, rho, mu_r = 1, epsilon_r = 1):
     return delta
     
 
-class CircularResistiveWall(Wakefield):
+class CircularResistiveWall(WakeField):
     """
-    Resistive wall Wakefield element for a circular beam pipe, approximated
+    Resistive wall WakeField element for a circular beam pipe, approximated
     formulas from Eq. (2.77) of Chao book [1].
     
     Parameters
diff --git a/collective_effects/resonator.py b/collective_effects/resonator.py
index 410a287a18d3a526b05b477429e9fbf385f73049..c6a083d86542067dd1e8916ab7bcfb3ee6592d79 100644
--- a/collective_effects/resonator.py
+++ b/collective_effects/resonator.py
@@ -1,20 +1,20 @@
 # -*- coding: utf-8 -*-
 """
 This module defines the impedances and wake functions from the resonator model 
-based on the Wakefield class.
+based on the WakeField class.
 
 @author: Alexis Gamelin, Watanyu Foosang
 @date: 25/09/2020
 """
 
 import numpy as np
-from mbtrack2.collective_effects.wakefield import (Wakefield, Impedance, 
+from mbtrack2.collective_effects.wakefield import (WakeField, Impedance, 
                                                    WakeFunction)
 
-class Resonator(Wakefield):
+class Resonator(WakeField):
     def __init__(self, Rs, fr, Q, plane, n_wake=1e4, n_imp=1e4, imp_freq_lim=50e9):
         """
-        Resonator model Wakefield element which computes the impedance and the 
+        Resonator model WakeField element which computes the impedance and the 
         wake function in both longitudinal and transverse case.
 
         Parameters
diff --git a/collective_effects/tapers.py b/collective_effects/tapers.py
index 34dde1855a4ab3b3b2fcb0f220da80010286ab3b..8d944e2b84c772ceb99823fc57783c41daabf209 100644
--- a/collective_effects/tapers.py
+++ b/collective_effects/tapers.py
@@ -8,11 +8,11 @@ Created on Tue Mar 31 13:15:36 2020
 from scipy.constants import mu_0, epsilon_0, c, pi
 import numpy as np
 from scipy.integrate import trapz
-from mbtrack2.collective_effects.wakefield import Wakefield, Impedance
+from mbtrack2.collective_effects.wakefield import WakeField, Impedance
 
-class StupakovRectangularTaper(Wakefield):
+class StupakovRectangularTaper(WakeField):
     """
-    Rectangular vertical taper Wakefield element, using the low frequency 
+    Rectangular vertical taper WakeField element, using the low frequency 
     approxmiation. Assume constant taper angle. Formulas from [1].
     
     ! Valid for low w
@@ -136,9 +136,9 @@ class StupakovRectangularTaper(Wakefield):
         
         return -1j*pi*self.Z0/4*integral
     
-class StupakovCircularTaper(Wakefield):
+class StupakovCircularTaper(WakeField):
     """
-    Circular taper Wakefield element, using the low frequency 
+    Circular taper WakeField element, using the low frequency 
     approxmiation. Assume constant taper angle. Formulas from [1].
     
     ! Valid for low w
diff --git a/collective_effects/utilities.py b/collective_effects/utilities.py
new file mode 100644
index 0000000000000000000000000000000000000000..8f792b2ce17194998a0aea4c524dbcf6135b7d91
--- /dev/null
+++ b/collective_effects/utilities.py
@@ -0,0 +1,393 @@
+# -*- coding: utf-8 -*-
+"""
+This module defines utilities functions, helping to deals with the 
+collective_effects module.
+
+@author: Alexis Gamelin
+@date: 28/09/2020
+"""
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+from mbtrack2.collective_effects.wakefield import Impedance
+
+def read_CST(file, impedance_type='long', divide_by=None):
+    """
+    Read CST file format into an Impedance object.
+    
+    Parameters
+    ----------
+    file : str
+        Path to the file to read.
+    impedance_type : str, optional
+        Type of the Impedance object.
+    divide_by : float, optional
+        Divide the impedance by a value. Mainly used to normalize transverse 
+        impedance by displacement.
+        
+    Returns
+    -------
+    result : Impedance object
+        Data from file.
+    """
+    df = pd.read_csv(file, comment="#", header = None, sep = "\t", 
+                    names = ["Frequency","Real","Imaginary"])
+    df["Frequency"] = df["Frequency"]*1e9 
+    if divide_by is not None:
+        df["Real"] = df["Real"]/divide_by
+        df["Imaginary"] = df["Imaginary"]/divide_by
+    if impedance_type == "long":
+        df["Real"] = np.abs(df["Real"])
+    df.set_index("Frequency", inplace = True)
+    result = Impedance(variable = df.index,
+                       function = df["Real"] + 1j*df["Imaginary"],
+                       impedance_type=impedance_type)
+    return result
+
+def read_IW2D(file, impedance_type='long'):
+    """
+    Read IW2D file format into an Impedance object.
+    
+    Parameters
+    ----------
+    file : str
+        Path to the file to read.
+    impedance_type : str, optional
+        Type of the Impedance object.
+        
+    Returns
+    -------
+    result : Impedance object
+        Data from file.
+    """
+    df = pd.read_csv(file, sep = ' ', header = None, 
+                     names = ["Frequency","Real","Imaginary"], skiprows=1)
+    df.set_index("Frequency", inplace = True)
+    result = Impedance(variable = df.index,
+                       function = df["Real"] + 1j*df["Imaginary"],
+                       impedance_type=impedance_type)
+    return result
+
+def spectral_density(frequency, sigma, m = 1, mode="Hermite"):
+    """
+    Compute the spectral density of different modes for various values of the
+    head-tail mode number, based on Table 1 p238 of [1].
+    
+    Parameters
+    ----------
+    frequency : list or numpy array
+        sample points of the spectral density in [Hz]
+    sigma : float
+        RMS bunch length in [s]
+    m : int, optional
+        head-tail (or azimutal/synchrotron) mode number
+    mode: str, optional
+        type of the mode taken into account for the computation:
+        -"Hermite" modes for Gaussian bunches
+
+    Returns
+    -------
+    numpy array
+    
+    References
+    ----------
+    [1] : Handbook of accelerator physics and engineering, 3rd printing.
+    """
+    
+    if mode == "Hermite":
+        return (2*np.pi*frequency*sigma)**(2*m)*np.exp(
+                -1*(2*np.pi*frequency*sigma)**2)
+    else:
+        raise NotImplementedError("Not implemanted yet.")
+        
+def gaussian_bunch_spectrum(frequency, sigma): 
+    """
+    Compute a Gaussian bunch spectrum [1].
+
+    Parameters
+    ----------
+    frequency : array
+        sample points of the beam spectrum in [Hz].
+    sigma : float
+        RMS bunch length in [s].
+
+    Returns
+    -------
+    bunch_spectrum : array
+        Bunch spectrum sampled at points given in frequency.
+        
+    References
+    ----------
+    [1] : Gamelin, A. (2018). Collective effects in a transient microbunching 
+    regime and ion cloud mitigation in ThomX. p86, Eq. 4.19
+    """
+    return np.exp(-1/2*(2*np.pi*frequency)**2*sigma**2)
+        
+        
+def beam_spectrum(frequency, M, bunch_spacing, sigma=None, 
+                  bunch_spectrum=None): 
+    """
+    Compute the beam spectrum assuming constant spacing between bunches [1].
+
+    Parameters
+    ----------
+    frequency : list or numpy array
+        sample points of the beam spectrum in [Hz].
+    M : int
+        Number of bunches.
+    bunch_spacing : float
+        Time between two bunches in [s].
+    sigma : float, optional
+        If bunch_spectrum is None then a Gaussian bunch with sigma RMS bunch 
+        length in [s] is assumed.
+    bunch_spectrum : array, optional
+        Bunch spectrum sampled at points given in frequency.
+
+    Returns
+    -------
+    beam_spectrum : array
+
+    References
+    ----------
+    [1] Rumolo, G - Beam Instabilities - CAS - CERN Accelerator School: 
+        Advanced Accelerator Physics Course - 2014, Eq. 9
+    """
+    
+    if bunch_spectrum is None:
+        bunch_spectrum = gaussian_bunch_spectrum(frequency, sigma)
+        
+    beam_spectrum = (bunch_spectrum * np.exp(1j*np.pi*frequency * 
+                                             bunch_spacing*(M-1)) * 
+                     np.sin(M*np.pi*frequency*bunch_spacing) / 
+                     np.sin(np.pi*frequency*bunch_spacing))
+    
+    return beam_spectrum
+    
+    
+def effective_impedance(ring, imp, m, mu, sigma, M, tuneS, xi=None, 
+                        mode="Hermite"):
+    """
+    Compute the effective (longitudinal or transverse) impedance. 
+    Formulas from Eq. (1) and (2) p238 of [1].
+    
+    Parameters
+    ----------
+    ring : Synchrotron object
+    imp : Impedance object
+    mu : int
+        coupled bunch mode number, goes from 0 to (M-1) where M is the
+        number of bunches
+    m : int
+        head-tail (or azimutal/synchrotron) mode number
+    sigma : float
+        RMS bunch length in [s]
+    M : int
+        Number of bunches.
+    tuneS : float
+        Synchrotron tune.
+    xi : float, optional
+        (non-normalized) chromaticity
+    mode: str, optional
+        type of the mode taken into account for the computation:
+        -"Hermite" modes for Gaussian bunches
+
+    Returns
+    -------
+    Zeff : float 
+        effective impedance in [ohm] or in [ohm/m] depanding on the impedance
+        type.
+        
+    References
+    ----------
+    [1] : Handbook of accelerator physics and engineering, 3rd printing.
+    """
+    
+    if not isinstance(imp, Impedance):
+        raise TypeError("{} should be an Impedance object.".format(imp))
+        
+    fmin = imp.data.index.min()
+    fmax = imp.data.index.max()
+    if fmin > 0:
+        double_sided_impedance(imp)
+        
+    if mode == "Hermite":
+        def h(f):
+            return spectral_density(frequency=f, sigma=sigma, m=m,
+                                    mode="Hermite")
+    else:
+        raise NotImplementedError("Not implemanted yet.")
+    
+    pmax = fmax/(ring.f0 * M) - 1
+    pmin = fmin/(ring.f0 * M) + 1
+    
+    p = np.arange(pmin,pmax+1)
+    
+    if imp.impedance_type == "long":
+        fp = ring.f0*(p*M + mu + m*tuneS)
+        fp = fp[np.nonzero(fp)] # Avoid division by 0
+        num = np.sum( imp(fp) * h(fp) / (fp*2*np.pi) )
+        den = np.sum( h(fp) )
+        Zeff = num/den
+        
+    elif imp.impedance_type == "xdip" or imp.impedance_type == "ydip":
+        if imp.impedance_type == "xdip":
+            tuneXY = ring.tune[0]
+            if xi is None :
+                xi = ring.chro[0]
+        elif imp.impedance_type == "ydip":
+            tuneXY = ring.tune[1]
+            if xi is None:
+                xi = ring.chro[1]
+        fp = ring.f0*(p*M + mu + tuneXY + m*tuneS)
+        f_xi = xi/ring.eta*ring.f0
+        num = np.sum( imp(fp) * h(fp - f_xi) )
+        den = np.sum( h(fp - f_xi) )
+        Zeff = num/den
+    else:
+        raise TypeError("Effective impedance is only defined for long, xdip"
+                        " and ydip impedance type.")
+        
+    return Zeff
+
+
+def yokoya_elliptic(x_radius , y_radius):
+    """
+    Compute Yokoya factors for an elliptic beam pipe.
+    Function adapted from N. Mounet IW2D.
+
+    Parameters
+    ----------
+    x_radius : float
+        Horizontal semi-axis of the ellipse in [m].
+    y_radius : float
+        Vertical semi-axis of the ellipse in [m].
+
+    Returns
+    -------
+    yoklong : float
+        Yokoya factor for the longitudinal impedance.
+    yokxdip : float
+        Yokoya factor for the dipolar horizontal impedance.
+    yokydip : float
+        Yokoya factor for the dipolar vertical impedance.
+    yokxquad : float
+        Yokoya factor for the quadrupolar horizontal impedance.
+    yokyquad : float
+        Yokoya factor for the quadrupolar vertical impedance.
+    """
+    if y_radius < x_radius:
+        small_semiaxis = y_radius
+        large_semiaxis = x_radius
+    else:
+        small_semiaxis = x_radius
+        large_semiaxis = y_radius
+
+    # read Yokoya factors interpolation file
+    # BEWARE: columns are ratio, dipy, dipx, quady, quadx
+    yokoya_file = pd.read_csv("collective_effects/data/" +
+                              "Yokoya_elliptic_from_Elias_USPAS.csv")
+    ratio_col = yokoya_file["x"]
+    # compute semi-axes ratio (first column of this file)
+    ratio = (large_semiaxis - small_semiaxis)/(large_semiaxis + small_semiaxis)
+
+    # interpolate Yokoya file at the correct ratio
+    yoklong = 1
+    
+    if y_radius < x_radius:
+        yokydip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
+        yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
+        yokyquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
+        yokxquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
+    else:
+        yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
+        yokydip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
+        yokxquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
+        yokyquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])        
+
+    return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)
+
+def beam_loss_factor(impedance, frequency, spectrum, ring):
+    """
+    Compute "beam" loss factor using the beam spectrum, uses a sum instead of 
+    integral compared to loss_factor [1].
+
+    Parameters
+    ----------
+    impedance : Impedance of type "long"
+    frequency : array
+        Sample points of spectrum.
+    spectrum : array
+        Beam spectrum to consider.
+    ring : Synchrotron object
+
+    Returns
+    -------
+    kloss_beam : float
+        Beam loss factor in [V/C].
+        
+    References
+    ----------
+    [1] : Handbook of accelerator physics and engineering, 3rd printing. 
+        Eq (3) p239.
+    """
+    pmax = np.floor(impedance.data.index.max()/ring.f0)
+    pmin = np.floor(impedance.data.index.min()/ring.f0)
+    
+    if pmin >= 0:
+        double_sided_impedance(impedance)
+        pmin = -1*pmax
+    
+    p = np.arange(pmin+1,pmax)    
+    pf0 = p*ring.f0
+    ReZ = np.real(impedance(pf0))
+    spectral_density = np.abs(spectrum)**2
+    # interpolation of the spectrum is needed to avoid problems liked to 
+    # division by 0
+    # computing the spectrum directly to the frequency points gives
+    # wrong results
+    spect = interp1d(frequency, spectral_density)
+    kloss_beam = ring.f0 * np.sum(ReZ*spect(pf0))
+    
+    return kloss_beam
+
+def double_sided_impedance(impedance):
+    """
+    Add negative frequency points to single sided impedance spectrum following
+    symetries depending on impedance type.
+
+    Parameters
+    ----------
+    impedance : Impedance object
+        Single sided impedance.
+    """
+    fmin = impedance.data.index.min()
+    
+    if fmin >= 0:
+        negative_index = impedance.data.index*-1
+        negative_data = impedance.data.set_index(negative_index)
+        
+        imp_type = impedance.impedance_type
+        
+        if imp_type == "long":
+            negative_data["imag"] = -1*negative_data["imag"]
+            
+        elif (imp_type == "xdip") or (imp_type == "ydip"):
+            negative_data["real"] = -1*negative_data["real"]
+        
+        elif (imp_type == "xquad") or (imp_type == "yquad"):
+            negative_data["real"] = -1*negative_data["real"]
+            
+        else:
+            raise ValueError("Wrong impedance type")
+            
+        try:
+            negative_data = negative_data.drop(0)
+        except KeyError:
+            pass
+            
+        all_data = impedance.data.append(negative_data)
+        all_data = all_data.sort_index()
+        impedance.data = all_data
+        
diff --git a/collective_effects/wakefield.py b/collective_effects/wakefield.py
index 7854a48b66a4d130476650d607763c1ff2b1d4af..13515477a3bb0db8ae77bb29b73f241efe7c4680 100644
--- a/collective_effects/wakefield.py
+++ b/collective_effects/wakefield.py
@@ -19,6 +19,7 @@ from scipy.integrate import trapz
 from scipy.constants import c
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 from mbtrack2.tracking.element import Element
+from mbtrack2.collective_effects.utilities import beam_spectrum, gaussian_bunch_spectrum, beam_loss_factor
 
 class ComplexData:
     """
@@ -564,16 +565,16 @@ class Impedance(ComplexData):
         super().__init__(variable=freq_trun, function=imp)
         self.data.index.name = "frequency [Hz]"    
     
-class Wakefield:
+class WakeField:
     """
-    Defines a Wakefield which corresponds to a single physical element which 
+    Defines a WakeField which corresponds to a single physical element which 
     produces different types of wakes, represented by their Impedance or 
     WakeFunction objects.
     
     Parameters
     ----------
     structure_list : list, optional
-        list of Impedance/WakeFunction objects to add to the Wakefield.
+        list of Impedance/WakeFunction objects to add to the WakeField.
     name : str, optional
     
     Attributes
@@ -584,13 +585,13 @@ class Wakefield:
     Methods
     -------
     append_to_model(structure_to_add)
-        Add Impedance/WakeFunction to Wakefield.
+        Add Impedance/WakeFunction to WakeField.
     list_to_attr(structure_list)
-        Add list of Impedance/WakeFunction to Wakefield.
+        Add list of Impedance/WakeFunction to WakeField.
     add_wakefields(wake1, beta1, wake2, beta2)
-        Add two Wakefields taking into account beta functions.
+        Add two WakeFields taking into account beta functions.
     add_several_wakefields(wakefields, beta)
-        Add a list of Wakefields taking into account beta functions.
+        Add a list of WakeFields taking into account beta functions.
     """
 
     def __init__(self, structure_list=None, name=None):
@@ -599,7 +600,7 @@ class Wakefield:
 
     def append_to_model(self, structure_to_add):
         """
-        Add Impedance/WakeFunction component to Wakefield.
+        Add Impedance/WakeFunction component to WakeField.
 
         Parameters
         ----------
@@ -625,7 +626,7 @@ class Wakefield:
     
     def list_to_attr(self, structure_list):
         """
-         Add list of Impedance/WakeFunction components to Wakefield.
+         Add list of Impedance/WakeFunction components to WakeField.
 
         Parameters
         ----------
@@ -652,23 +653,23 @@ class Wakefield:
     @staticmethod
     def add_wakefields(wake1, beta1, wake2, beta2):
         """
-        Add two Wakefields taking into account beta functions.
+        Add two WakeFields taking into account beta functions.
 
         Parameters
         ----------
-        wake1 : Wakefield
-            Wakefield to add.
+        wake1 : WakeField
+            WakeField to add.
         beta1 : array of shape (2,)
             Beta functions at wake1 postion.
-        wake2 : Wakefield
-            Wakefield to add.
+        wake2 : WakeField
+            WakeField to add.
         beta2 : array of shape (2,)
             Beta functions at wake2 postion.
 
         Returns
         -------
-        wake_sum : Wakefield
-            Wakefield with summed components.
+        wake_sum : WakeField
+            WakeField with summed components.
 
         """
         wake_sum = deepcopy(wake1)
@@ -694,28 +695,28 @@ class Wakefield:
     @staticmethod
     def add_several_wakefields(wakefields, beta):
         """
-        Add a list of Wakefields taking into account beta functions.
+        Add a list of WakeFields taking into account beta functions.
         
         Parameters
         ----------
-        wakefields : list of Wakefield
-            Wakefields to add.
+        wakefields : list of WakeField
+            WakeFields to add.
         beta : array of shape (len(wakefields), 2)
-            Beta functions at the Wakefield postions.
+            Beta functions at the WakeField postions.
 
         Returns
         -------
-        wake_sum : Wakefield
-            Wakefield with summed components..
+        wake_sum : WakeField
+            WakeField with summed components..
 
         """
         if len(wakefields) == 1:
             return wakefields[0]
         elif len(wakefields) > 1:
-            wake_sum = Wakefield.add_wakefields(wakefields[0], beta[:,0],
+            wake_sum = WakeField.add_wakefields(wakefields[0], beta[:,0],
                                      wakefields[1], beta[:,1])
             for i in range(len(wakefields) - 2):
-                wake_sum = Wakefield.add_wakefields(wake_sum, [1 ,1], 
+                wake_sum = WakeField.add_wakefields(wake_sum, [1 ,1], 
                                          wakefields[i+2], beta[:,i+2])
             return wake_sum
         else:
@@ -728,25 +729,25 @@ class ImpedanceModel(Element):
     Parameters
     ----------
     ring : Synchrotron object
-    wakefield_list : list of Wakefield objects
-        Wakefields to add to the model.
+    wakefield_list : list of WakeField objects
+        WakeFields to add to the model.
     wakefiled_postions : list
         Longitudinal positions corresponding to the added Wakfields.
     
     Attributes
     ----------
-    wakefields : list of Wakefield objects
-        Wakefields in the model.
+    wakefields : list of WakeField objects
+        WakeFields in the model.
     positions : array
-        Wakefields positions.
+        WakeFields positions.
     names : array
-        Names of (unique) Wakefield objects.
-    sum : Wakefield
-        Sum of every Wakefield in the model.
-    sum_"name" : Wakefield
+        Names of (unique) WakeField objects.
+    sum : WakeField
+        Sum of every WakeField in the model.
+    sum_"name" : WakeField
         Sum of every Wakefield with the same "name".
     sum_names : array
-        Names of attributes where the wakefields are summed by name.    
+        Names of attributes where the WakeFields are summed by name.    
     
     Methods
     -------
@@ -755,17 +756,17 @@ class ImpedanceModel(Element):
     add_multiple_elements(wakefield, wakefiled_postions)
         Add the same element at different locations to the model.
     sum_elements()
-        Sum all wakefields into self.sum.
+        Sum all WakeFields into self.sum.
     update_name_list()
         Update self.names with uniques names of self.wakefields.
     find_wakefield(name)
-        Return indexes of wakefields with the same name in self.wakefields.
+        Return indexes of WakeFields with the same name in self.wakefields.
     sum_by_name(name)
         Sum the elements with the same name in the model into sum_name.
     sum_by_name_all(name)
         Sum all the elements with the same name in the model into sum_name.
     plot_area(Z_type="Zlong", component="real", sigma=None, attr_list=None)
-        Plot the contributions of different kind of wakefields.
+        Plot the contributions of different kind of WakeFields.
     """
     
     def __init__(self, ring, wakefield_list=None, wakefiled_postions=None):
@@ -788,9 +789,9 @@ class ImpedanceModel(Element):
         raise NotImplementedError
         
     def sum_elements(self):
-        """Sum all wakefields into self.sum"""
+        """Sum all WakeFields into self.sum"""
         beta = self.optics.beta(self.positions)
-        self.sum = Wakefield.add_several_wakefields(self.wakefields, beta)
+        self.sum = WakeField.add_several_wakefields(self.wakefields, beta)
         if "sum" not in self.sum_names:
             self.sum_names = np.append(self.sum_names, "sum")
     
@@ -803,7 +804,7 @@ class ImpedanceModel(Element):
                 self.names = np.append(self.names, wakefield.name)
                 
     def count_elements(self):
-        """Count number of each type of wakefield in the model."""
+        """Count number of each type of WakeField in the model."""
         self.count = np.zeros(len(self.names))
         for wakefield in self.wakefields:
             if wakefield.name is not None:
@@ -811,12 +812,12 @@ class ImpedanceModel(Element):
                 
     def find_wakefield(self, name):
         """
-        Return indexes of wakefields with the same name in self.wakefields.
+        Return indexes of WakeFields with the same name in self.wakefields.
 
         Parameters
         ----------
         name : str
-            Wakefield name.
+            WakeField name.
 
         Returns
         -------
@@ -837,7 +838,7 @@ class ImpedanceModel(Element):
         Parameters
         ----------
         name : str
-            Name of the wakefield to sum.
+            Name of the WakeField to sum.
         """
         attribute_name = "sum_" + name
         index = self.find_wakefield(name)
@@ -845,7 +846,7 @@ class ImpedanceModel(Element):
         wakes = []
         for i in index:
             wakes.append(self.wakefields[i])
-        wake_sum = Wakefield.add_several_wakefields(wakes, beta)
+        wake_sum = WakeField.add_several_wakefields(wakes, beta)
         setattr(self, attribute_name, wake_sum)
         if attribute_name not in self.sum_names:
             self.sum_names = np.append(self.sum_names, attribute_name)
@@ -863,8 +864,8 @@ class ImpedanceModel(Element):
 
         Parameters
         ----------
-        wakefield_list : list of Wakefield objects
-            Wakefields to add to the model.
+        wakefield_list : list of WakeField objects
+            WakeFields to add to the model.
         wakefiled_postions : list
             Longitudinal positions corresponding to the added Wakfields.
         """
@@ -883,8 +884,8 @@ class ImpedanceModel(Element):
 
         Parameters
         ----------
-        wakefield : Wakefield object
-            Wakefield to add to the model.
+        WakeField : WakeField object
+            WakeField to add to the model.
         wakefiled_postions : list
             Longitudinal positions corresponding to the added Wakfield.
         """
@@ -897,7 +898,7 @@ class ImpedanceModel(Element):
     def plot_area(self, Z_type="Zlong", component="real", sigma=None, 
                   attr_list=None, zoom=False):
         """
-        Plot the contributions of different kind of wakefields.
+        Plot the contributions of different kind of WakeFields.
 
         Parameters
         ----------
@@ -1128,395 +1129,3 @@ class ImpedanceModel(Element):
         summary["P (bunch) [W]"] = summary["loss factor (bunch) [V/pC]"]*1e12*Q**2/(self.ring.T0)*M
                 
         return summary
-        
-    
-def read_CST(file, impedance_type='long', divide_by=None):
-    """
-    Read CST file format into an Impedance object.
-    
-    Parameters
-    ----------
-    file : str
-        Path to the file to read.
-    impedance_type : str, optional
-        Type of the Impedance object.
-    divide_by : float, optional
-        Divide the impedance by a value. Mainly used to normalize transverse 
-        impedance by displacement.
-        
-    Returns
-    -------
-    result : Impedance object
-        Data from file.
-    """
-    df = pd.read_csv(file, comment="#", header = None, sep = "\t", 
-                    names = ["Frequency","Real","Imaginary"])
-    df["Frequency"] = df["Frequency"]*1e9 
-    if divide_by is not None:
-        df["Real"] = df["Real"]/divide_by
-        df["Imaginary"] = df["Imaginary"]/divide_by
-    if impedance_type == "long":
-        df["Real"] = np.abs(df["Real"])
-    df.set_index("Frequency", inplace = True)
-    result = Impedance(variable = df.index,
-                       function = df["Real"] + 1j*df["Imaginary"],
-                       impedance_type=impedance_type)
-    return result
-
-def read_IW2D(file, impedance_type='long'):
-    """
-    Read IW2D file format into an Impedance object.
-    
-    Parameters
-    ----------
-    file : str
-        Path to the file to read.
-    impedance_type : str, optional
-        Type of the Impedance object.
-        
-    Returns
-    -------
-    result : Impedance object
-        Data from file.
-    """
-    df = pd.read_csv(file, sep = ' ', header = None, 
-                     names = ["Frequency","Real","Imaginary"], skiprows=1)
-    df.set_index("Frequency", inplace = True)
-    result = Impedance(variable = df.index,
-                       function = df["Real"] + 1j*df["Imaginary"],
-                       impedance_type=impedance_type)
-    return result
-
-def spectral_density(frequency, sigma, m = 1, mode="Hermite"):
-    """
-    Compute the spectral density of different modes for various values of the
-    head-tail mode number, based on Table 1 p238 of [1].
-    
-    Parameters
-    ----------
-    frequency : list or numpy array
-        sample points of the spectral density in [Hz]
-    sigma : float
-        RMS bunch length in [s]
-    m : int, optional
-        head-tail (or azimutal/synchrotron) mode number
-    mode: str, optional
-        type of the mode taken into account for the computation:
-        -"Hermite" modes for Gaussian bunches
-
-    Returns
-    -------
-    numpy array
-    
-    References
-    ----------
-    [1] : Handbook of accelerator physics and engineering, 3rd printing.
-    """
-    
-    if mode == "Hermite":
-        return (2*np.pi*frequency*sigma)**(2*m)*np.exp(
-                -1*(2*np.pi*frequency*sigma)**2)
-    else:
-        raise NotImplementedError("Not implemanted yet.")
-        
-def gaussian_bunch_spectrum(frequency, sigma): 
-    """
-    Compute a Gaussian bunch spectrum [1].
-
-    Parameters
-    ----------
-    frequency : array
-        sample points of the beam spectrum in [Hz].
-    sigma : float
-        RMS bunch length in [s].
-
-    Returns
-    -------
-    bunch_spectrum : array
-        Bunch spectrum sampled at points given in frequency.
-        
-    References
-    ----------
-    [1] : Gamelin, A. (2018). Collective effects in a transient microbunching 
-    regime and ion cloud mitigation in ThomX. p86, Eq. 4.19
-    """
-    return np.exp(-1/2*(2*np.pi*frequency)**2*sigma**2)
-        
-        
-def beam_spectrum(frequency, M, bunch_spacing, sigma=None, 
-                  bunch_spectrum=None): 
-    """
-    Compute the beam spectrum assuming constant spacing between bunches [1].
-
-    Parameters
-    ----------
-    frequency : list or numpy array
-        sample points of the beam spectrum in [Hz].
-    M : int
-        Number of bunches.
-    bunch_spacing : float
-        Time between two bunches in [s].
-    sigma : float, optional
-        If bunch_spectrum is None then a Gaussian bunch with sigma RMS bunch 
-        length in [s] is assumed.
-    bunch_spectrum : array, optional
-        Bunch spectrum sampled at points given in frequency.
-
-    Returns
-    -------
-    beam_spectrum : array
-
-    References
-    ----------
-    [1] Rumolo, G - Beam Instabilities - CAS - CERN Accelerator School: 
-        Advanced Accelerator Physics Course - 2014, Eq. 9
-    """
-    
-    if bunch_spectrum is None:
-        bunch_spectrum = gaussian_bunch_spectrum(frequency, sigma)
-        
-    beam_spectrum = (bunch_spectrum * np.exp(1j*np.pi*frequency * 
-                                             bunch_spacing*(M-1)) * 
-                     np.sin(M*np.pi*frequency*bunch_spacing) / 
-                     np.sin(np.pi*frequency*bunch_spacing))
-    
-    return beam_spectrum
-    
-    
-def effective_impedance(ring, imp, m, mu, sigma, M, tuneS, xi=None, 
-                        mode="Hermite"):
-    """
-    Compute the effective (longitudinal or transverse) impedance. 
-    Formulas from Eq. (1) and (2) p238 of [1].
-    
-    Parameters
-    ----------
-    ring : Synchrotron object
-    imp : Impedance object
-    mu : int
-        coupled bunch mode number, goes from 0 to (M-1) where M is the
-        number of bunches
-    m : int
-        head-tail (or azimutal/synchrotron) mode number
-    sigma : float
-        RMS bunch length in [s]
-    M : int
-        Number of bunches.
-    tuneS : float
-        Synchrotron tune.
-    xi : float, optional
-        (non-normalized) chromaticity
-    mode: str, optional
-        type of the mode taken into account for the computation:
-        -"Hermite" modes for Gaussian bunches
-
-    Returns
-    -------
-    Zeff : float 
-        effective impedance in [ohm] or in [ohm/m] depanding on the impedance
-        type.
-        
-    References
-    ----------
-    [1] : Handbook of accelerator physics and engineering, 3rd printing.
-    """
-    
-    if not isinstance(imp, Impedance):
-        raise TypeError("{} should be an Impedance object.".format(imp))
-        
-    fmin = imp.data.index.min()
-    fmax = imp.data.index.max()
-    if fmin > 0:
-        double_sided_impedance(imp)
-        
-    if mode == "Hermite":
-        def h(f):
-            return spectral_density(frequency=f, sigma=sigma, m=m,
-                                    mode="Hermite")
-    else:
-        raise NotImplementedError("Not implemanted yet.")
-    
-    pmax = fmax/(ring.f0 * M) - 1
-    pmin = fmin/(ring.f0 * M) + 1
-    
-    p = np.arange(pmin,pmax+1)
-    
-    if imp.impedance_type == "long":
-        fp = ring.f0*(p*M + mu + m*tuneS)
-        fp = fp[np.nonzero(fp)] # Avoid division by 0
-        num = np.sum( imp(fp) * h(fp) / (fp*2*np.pi) )
-        den = np.sum( h(fp) )
-        Zeff = num/den
-        
-    elif imp.impedance_type == "xdip" or imp.impedance_type == "ydip":
-        if imp.impedance_type == "xdip":
-            tuneXY = ring.tune[0]
-            if xi is None :
-                xi = ring.chro[0]
-        elif imp.impedance_type == "ydip":
-            tuneXY = ring.tune[1]
-            if xi is None:
-                xi = ring.chro[1]
-        fp = ring.f0*(p*M + mu + tuneXY + m*tuneS)
-        f_xi = xi/ring.eta*ring.f0
-        num = np.sum( imp(fp) * h(fp - f_xi) )
-        den = np.sum( h(fp - f_xi) )
-        Zeff = num/den
-    else:
-        raise TypeError("Effective impedance is only defined for long, xdip"
-                        " and ydip impedance type.")
-        
-    return Zeff
-
-
-def yokoya_elliptic(x_radius , y_radius):
-    """
-    Compute Yokoya factors for an elliptic beam pipe.
-    Function adapted from N. Mounet IW2D.
-
-    Parameters
-    ----------
-    x_radius : float
-        Horizontal semi-axis of the ellipse in [m].
-    y_radius : float
-        Vertical semi-axis of the ellipse in [m].
-
-    Returns
-    -------
-    yoklong : float
-        Yokoya factor for the longitudinal impedance.
-    yokxdip : float
-        Yokoya factor for the dipolar horizontal impedance.
-    yokydip : float
-        Yokoya factor for the dipolar vertical impedance.
-    yokxquad : float
-        Yokoya factor for the quadrupolar horizontal impedance.
-    yokyquad : float
-        Yokoya factor for the quadrupolar vertical impedance.
-    """
-    if y_radius < x_radius:
-        small_semiaxis = y_radius
-        large_semiaxis = x_radius
-    else:
-        small_semiaxis = x_radius
-        large_semiaxis = y_radius
-
-    # read Yokoya factors interpolation file
-    # BEWARE: columns are ratio, dipy, dipx, quady, quadx
-    yokoya_file = pd.read_csv("collective_effects/data/" +
-                              "Yokoya_elliptic_from_Elias_USPAS.csv")
-    ratio_col = yokoya_file["x"]
-    # compute semi-axes ratio (first column of this file)
-    ratio = (large_semiaxis - small_semiaxis)/(large_semiaxis + small_semiaxis)
-
-    # interpolate Yokoya file at the correct ratio
-    yoklong = 1
-    
-    if y_radius < x_radius:
-        yokydip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
-        yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
-        yokyquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
-        yokxquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
-    else:
-        yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
-        yokydip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
-        yokxquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
-        yokyquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])        
-
-    return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)
-
-def beam_loss_factor(impedance, frequency, spectrum, ring):
-    """
-    Compute "beam" loss factor using the beam spectrum, uses a sum instead of 
-    integral compared to loss_factor [1].
-
-    Parameters
-    ----------
-    impedance : Impedance of type "long"
-    frequency : array
-        Sample points of spectrum.
-    spectrum : array
-        Beam spectrum to consider.
-    ring : Synchrotron object
-
-    Returns
-    -------
-    kloss_beam : float
-        Beam loss factor in [V/C].
-        
-    References
-    ----------
-    [1] : Handbook of accelerator physics and engineering, 3rd printing. 
-        Eq (3) p239.
-    """
-    pmax = np.floor(impedance.data.index.max()/ring.f0)
-    pmin = np.floor(impedance.data.index.min()/ring.f0)
-    
-    if pmin >= 0:
-        double_sided_impedance(impedance)
-        pmin = -1*pmax
-    
-    p = np.arange(pmin+1,pmax)    
-    pf0 = p*ring.f0
-    ReZ = np.real(impedance(pf0))
-    spectral_density = np.abs(spectrum)**2
-    # interpolation of the spectrum is needed to avoid problems liked to 
-    # division by 0
-    # computing the spectrum directly to the frequency points gives
-    # wrong results
-    spect = interp1d(frequency, spectral_density)
-    kloss_beam = ring.f0 * np.sum(ReZ*spect(pf0))
-    
-    return kloss_beam
-
-def double_sided_impedance(impedance):
-    """
-    Add negative frequency points to single sided impedance spectrum following
-    symetries depending on impedance type.
-
-    Parameters
-    ----------
-    impedance : Impedance object
-        Single sided impedance.
-    """
-    fmin = impedance.data.index.min()
-    
-    if fmin >= 0:
-        negative_index = impedance.data.index*-1
-        negative_data = impedance.data.set_index(negative_index)
-        
-        imp_type = impedance.impedance_type
-        
-        if imp_type == "long":
-            negative_data["imag"] = -1*negative_data["imag"]
-            
-        elif (imp_type == "xdip") or (imp_type == "ydip"):
-            negative_data["real"] = -1*negative_data["real"]
-        
-        elif (imp_type == "xquad") or (imp_type == "yquad"):
-            negative_data["real"] = -1*negative_data["real"]
-            
-        else:
-            raise ValueError("Wrong impedance type")
-            
-        try:
-            negative_data = negative_data.drop(0)
-        except KeyError:
-            pass
-            
-        all_data = impedance.data.append(negative_data)
-        all_data = all_data.sort_index()
-        impedance.data = all_data
-        
-
-        
-        
-        
-        
-        
-        
-        
-        
-        
-        
-