diff --git a/collective_effects/wakefield.py b/collective_effects/wakefield.py
index ea3fe51dc102b035efaa385d5c4c1ec40e0cfb42..1a1ebb3344d6a57b16d3d3cac83540c8c53ac6f5 100644
--- a/collective_effects/wakefield.py
+++ b/collective_effects/wakefield.py
@@ -71,14 +71,13 @@ class ComplexData:
         # from the two input data
 
         if isinstance(structure_to_add, (int, float, complex)):
-            structure_to_add = ComplexData(
-                                variable=self.data.index,
-                                function=structure_to_add * np.ones(len(self.data.index)))
-
-        initial_data = pd.concat(
-                        [self.data,
-                         structure_to_add.data.index.to_frame().set_index(index_name)],
-                        sort=True)
+            structure_to_add = ComplexData(variable=self.data.index,
+                                           function=(structure_to_add * 
+                                                     np.ones(len(self.data.index))))
+                                
+        data_to_concat = structure_to_add.data.index.to_frame().set_index(index_name)
+        
+        initial_data = pd.concat([self.data, data_to_concat], sort=True)
         initial_data = initial_data[~initial_data.index.duplicated(keep='first')]
 
         data_to_add = pd.concat(
@@ -325,16 +324,41 @@ class Impedance(ComplexData):
     
 class Wakefield:
     """
-    Define a general Wakefield machine element.
-    It contains Imepdance or WakeFunction objects which defines the wakefield
-    and other informations: beta functions.
+    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.
+    name : str, optional
+    
+    Attributes
+    ----------
+    impedance_components : array of str
+        Impedance component names present for this element.
+        
+    Methods
+    -------
+    append_to_model(structure_to_add)
+        Add Impedance/WakeFunction to Wakefield.
+    list_to_attr(structure_list)
+        Add list of Impedance/WakeFunction to Wakefield.
     """
 
-    def __init__(self, structure_list=[], name=None):
+    def __init__(self, structure_list=None, name=None):
         self.list_to_attr(structure_list)
         self.name = name
 
     def append_to_model(self, structure_to_add):
+        """
+        Add Impedance/WakeFunction component to Wakefield.
+
+        Parameters
+        ----------
+        structure_to_add : Impedance or WakeFunction object
+        """
         list_of_attributes = dir(self)
         if isinstance(structure_to_add, Impedance):
             attribute_name = "Z" + structure_to_add.impedance_type
@@ -354,26 +378,46 @@ class Wakefield:
             raise ValueError("{} is not an Impedance nor a WakeFunction.".format(structure_to_add))
     
     def list_to_attr(self, structure_list):
-        for component in structure_list:
-            self.append_to_model(component)
+        """
+         Add list of Impedance/WakeFunction components to Wakefield.
+
+        Parameters
+        ----------
+        structure_list : list of Impedance or WakeFunction objects.
+        """
+        if structure_list is not None:
+            for component in structure_list:
+                self.append_to_model(component)
     
     @property
-    def list_impedance_components(self):
+    def impedance_components(self):
         """
-        Return the list of impedance components for the element,
-        based on the attributes names.
+        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)])
     
     
 class ImpedanceModel(Element):
+    """
+    Define the impedance model of the machine.
+    
+    Parameters
+    ----------
+    ring : Synchrotron object
+    wakefield_list : list of Wakefield objects
+        Wakefields to add to the model.
+    wakefiled_postions : list
+        Longitudinal positions corresponding to the added Wakfields.
+        
+    
+    """
     
-    def __init__(self, ring, wakefields_to_add=None, wakefileds_postions=None):
+    def __init__(self, ring, wakefield_list=None, wakefiled_postions=None):
         self.ring = ring
         self.optics = self.ring.optics
         self.wakefields = []
         self.positions = np.array([])
-        self.add(wakefields_to_add, wakefileds_postions)
+        self.add(wakefield_list, wakefiled_postions)
         
     def track(self, beam):
         """
@@ -386,7 +430,9 @@ class ImpedanceModel(Element):
         raise NotImplementedError
         
     def sum_elements(self):
-        
+        """
+        Sum all the elements in the model into sum_wakefield.
+        """
         beta = self.optics.beta(self.positions)
         list_impedance_components = ["Zlong","Zxdip","Zydip","Zxquad","Zyquad"]
         sum_wakefield = Wakefield()
@@ -405,11 +451,21 @@ class ImpedanceModel(Element):
             sum_wakefield.append_to_model(sum_imp)
         self.sum_wakefield = sum_wakefield
         
-    def add(self, wakefield_list, position_list):
+    def add(self, wakefield_list, wakefiled_postions):
+        """
+        Add elements to the model.
+
+        Parameters
+        ----------
+        wakefield_list : list of Wakefield objects
+            Wakefields to add to the model.
+        wakefiled_postions : list
+            Longitudinal positions corresponding to the added Wakfields.
+        """
         
-        if (wakefield_list is not None) and (position_list is not None):
+        if (wakefield_list is not None) and (wakefiled_postions is not None):
             for wakefield in wakefield_list:
                 self.wakefields.append(wakefield)
                 
-            for position in position_list:
+            for position in wakefiled_postions:
                 self.positions = np.append(self.positions, position)
\ No newline at end of file
diff --git a/tracking/optics.py b/tracking/optics.py
index 2482a373478ba532fb17b22db7cf4286324d0606..2e3483657181ebab0073f5d9eb135ac7808d07dd 100644
--- a/tracking/optics.py
+++ b/tracking/optics.py
@@ -1,8 +1,10 @@
 # -*- coding: utf-8 -*-
 """
-Created on Wed Mar 11 18:00:28 2020
+Module where the class to store the optic functions and the lattice physical 
+parameters are defined.
 
-@author: gamelina
+@author: Alexis Gamelin
+@date: 10/04/2020
 """
 
 import at
@@ -12,7 +14,49 @@ from scipy.interpolate import interp1d
 
 
 class Optics:
-    """Handle optic functions"""
+    """
+    Class used to handle optic functions.
+    
+    Parameters
+    ----------
+    lattice_file : str, optional if local_beta, local_alpha and 
+        local_dispersion are specified.
+        AT lattice file path.
+    local_beta : array of shape (2,), optional if lattice_file is specified.
+        Beta function at the location of the tracking. Default is mean beta if
+        lattice has been loaded.
+    local_alpha : array of shape (2,), optional if lattice_file is specified.
+        Alpha function at the location of the tracking. Default is mean alpha 
+        if lattice has been loaded.
+    local_dispersion : array of shape (4,), optional if lattice_file is 
+        specified.
+        Dispersion function and derivative at the location of the tracking. 
+        Default is zero if lattice has been loaded.  
+        
+    Attributes
+    ----------
+    use_local_values : bool
+        True if no lattice has been loaded.
+    local_gamma : array of shape (2,)
+        Gamma function at the location of the tracking.
+    lattice : AT lattice
+        
+    Methods
+    -------
+    load_from_AT(lattice_file, **kwargs)
+        Load a lattice from accelerator toolbox (AT).
+    setup_interpolation()
+        Setup interpolation of the optic functions.
+    beta(position)
+        Return beta functions at specific locations given by position.
+    alpha(position)
+        Return alpha functions at specific locations given by position.
+    gamma(position)
+        Return gamma functions at specific locations given by position.
+    dispersion(position)
+        Return dispersion functions at specific locations given by position.
+    """
+    
     def __init__(self, lattice_file=None, local_beta=None, local_alpha=None, 
                  local_dispersion=None, **kwargs):
         
@@ -41,7 +85,19 @@ class Optics:
             self.local_dispersion = local_dispersion
 
     def load_from_AT(self, lattice_file, **kwargs):
-        
+        """
+        Load a lattice from accelerator toolbox (AT).
+
+        Parameters
+        ----------
+        lattice_file : str
+            AT lattice file path.
+        n_points : int or float, optional
+            Minimum number of points to use for the optic function arrays.
+        periodicity : int, optional
+            Lattice periodicity, if not specified the AT lattice periodicity is
+            used.
+        """
         self.n_points = int(kwargs.get("n_points", 1e3))
         periodicity = kwargs.get("periodicity")
         
@@ -62,8 +118,17 @@ class Optics:
         self.position = pos
         self.beta_array = np.tile(twiss.beta.T, self.periodicity)
         self.alpha_array = np.tile(twiss.alpha.T, self.periodicity)
-        self.gamma_array = (1 + self.alpha_array**2)/self.beta_array
         self.dispersion_array = np.tile(twiss.dispersion.T, self.periodicity)
+        
+        self.position = np.append(self.position, self.lattice.circumference)
+        self.beta_array = np.append(self.beta_array, self.beta_array[:,0:1],
+                                    axis=1)
+        self.alpha_array = np.append(self.alpha_array, self.alpha_array[:,0:1],
+                                     axis=1)
+        self.dispersion_array = np.append(self.dispersion_array,
+                                          self.dispersion_array[:,0:1], axis=1)
+        
+        self.gamma_array = (1 + self.alpha_array**2)/self.beta_array
         self.tune = tune * self.periodicity
         self.chro = chrom * self.periodicity
         self.ac = at.get_mcf(self.lattice)
@@ -71,7 +136,8 @@ class Optics:
         self.setup_interpolation()
         
         
-    def setup_interpolation(self):        
+    def setup_interpolation(self):      
+        """Setup interpolation of the optic functions."""
         self.betaX = interp1d(self.position, self.beta_array[0,:],
                               kind='linear')
         self.betaY = interp1d(self.position, self.beta_array[1,:],
@@ -94,6 +160,20 @@ class Optics:
                                kind='linear')
         
     def beta(self, position):
+        """
+        Return beta functions at specific locations given by position. If no
+        lattice has been loaded, local values are returned.
+
+        Parameters
+        ----------
+        position : array or float
+            Longitudinal position at which the beta functions are returned.
+
+        Returns
+        -------
+        beta : array
+            Beta functions.
+        """
         if self.use_local_values:
             return np.outer(self.local_beta, np.ones_like(position))
         else:
@@ -101,6 +181,20 @@ class Optics:
             return np.array(beta)
     
     def alpha(self, position):
+        """
+        Return alpha functions at specific locations given by position. If no
+        lattice has been loaded, local values are returned.
+
+        Parameters
+        ----------
+        position : array or float
+            Longitudinal position at which the alpha functions are returned.
+
+        Returns
+        -------
+        alpha : array
+            Alpha functions.
+        """
         if self.use_local_values:
             return np.outer(self.local_alpha, np.ones_like(position))
         else:
@@ -108,6 +202,20 @@ class Optics:
             return np.array(alpha)
     
     def gamma(self, position):
+        """
+        Return gamma functions at specific locations given by position. If no
+        lattice has been loaded, local values are returned.
+
+        Parameters
+        ----------
+        position : array or float
+            Longitudinal position at which the gamma functions are returned.
+
+        Returns
+        -------
+        gamma : array
+            Gamma functions.
+        """
         if self.use_local_values:
             return np.outer(self.local_gamma, np.ones_like(position))
         else:
@@ -115,6 +223,21 @@ class Optics:
             return np.array(gamma)
     
     def dispersion(self, position):
+        """
+        Return dispersion functions at specific locations given by position. 
+        If no lattice has been loaded, local values are returned.
+
+        Parameters
+        ----------
+        position : array or float
+            Longitudinal position at which the dispersion functions are 
+            returned.
+
+        Returns
+        -------
+        dispersion : array
+            Dispersion functions.
+        """
         if self.use_local_values:
             return np.outer(self.local_dispersion, np.ones_like(position))
         else:
@@ -123,70 +246,271 @@ class Optics:
             return np.array(dispersion)
     
 class PhyisicalModel:
-    """Store lattice physical parameters"""
-    def __init__(self, ring, a0=5e-3, b0=5e-3, shape='circ', rho=1.7e-8,
-                 n_points=1e4):
+    """
+    Store the lattice physical parameters such as apperture and resistivity.
+    
+    Parameters
+    ----------
+    ring : Synchrotron object
+    x_right : float
+        Default value for the right horizontal aperture in [m].
+    y_top : float
+        Default value for the top vertical aperture in [m].
+    shape : str
+        Default value for the shape of the chamber cross section 
+        (elli/circ/rect).
+    rho : float
+        Default value for the resistivity of the chamber material [ohm.m].
+    x_left : float, optional
+        Default value for the left horizontal aperture in [m].
+    y_bottom : float, optional
+        Default value for the bottom vertical aperture in [m]. 
+    n_points : int or float, optional
+        Number of points to use in class arrays
+        
+    Attributes
+    ----------
+    position : array of shape (n_points,)
+        Longitudinal position in [m].
+    x_right : array of shape (n_points,)
+        Right horizontal aperture in [m].
+    y_top : array of shape (n_points,)
+        Top vertical aperture in [m].
+    shape : array of shape (n_points - 1,)
+        Shape of the chamber cross section (elli/circ/rect).
+    rho : array of shape (n_points - 1,)
+        Resistivity of the chamber material in [ohm.m].
+    x_left : array of shape (n_points,)
+        Left horizontal aperture in [m].
+    y_bottom : array of shape (n_points,)
+        Bottom vertical aperture in [m].
+    length : array of shape (n_points - 1,)
+        Length of each segments between two longitudinal positions in [m].
+    center : array of shape (n_points - 1,)
+        Center of each segments between two longitudinal positions in [m].
+        
+    Methods
+    -------
+    change_values(start_position, end_position, x_right, y_top, shape, rho)
+        Change the physical parameters between start_position and end_position.
+    taper(start_position, end_position, x_right_start, x_right_end, 
+          y_top_start, y_top_end, shape, rho)
+        Change the physical parameters to have a tapered transition between 
+        start_position and end_position.
+    plot_aperture()
+        Plot horizontal and vertical apertures.
+    resistive_wall_effective_radius(optics)
+        Return the effective radius of the chamber for resistive wall 
+        calculations.
+    """
+    def __init__(self, ring, x_right, y_top, shape, rho, x_left=None, 
+                 y_bottom=None, n_points=1e4):
         
         self.n_points = int(n_points)
         self.position = np.linspace(0, ring.L, self.n_points)
-        self.x0 = np.ones_like(self.position)*a0 # Horizontal half aperture [m]
-        self.y0 = np.ones_like(self.position)*b0 # Vertical half aperture [m]
+        self.x_right = np.ones_like(self.position)*x_right 
+        self.y_top = np.ones_like(self.position)*y_top
+        
+        if x_left is None:
+            self.x_left = np.ones_like(self.position)*-1*x_right
+        else:
+            self.x_left = np.ones_like(self.position)*x_left
+        
+        if y_bottom is None:
+            self.y_bottom = np.ones_like(self.position)*-1*y_top
+        else:
+            self.y_bottom = np.ones_like(self.position)*y_bottom
         
         self.length = self.position[1:] - self.position[:-1]
         self.center = (self.position[1:] + self.position[:-1])/2
-        self.rho = np.ones_like(self.center)*rho # Resistivity of the chamber material [ohm*m]
+        self.rho = np.ones_like(self.center)*rho
         
-        self.shape = np.repeat(np.array([shape]), self.n_points-1) # Shape of the chamber cross section (elli/circ/rect)
+        self.shape = np.repeat(np.array([shape]), self.n_points-1)
+
+    def resistive_wall_effective_radius(self, optics, x_right=True, 
+                                        y_top=True):
+        """
+        Return the effective radius of the chamber for resistive wall 
+        calculations as defined in Eq. 27 of [1].
+
+        Parameters
+        ----------
+        optics : Optics object
+        x_right : bool, optional
+            If True, x_right is used, if Fasle, x_left is used.
+        y_top : TYPE, optional
+            If True, y_top is used, if Fasle, y_bottom is used.
 
-    def resistive_wall_coef(self):
-        L = self.end_position[-1]
+        Returns
+        -------
+        rho_star : float
+            Effective resistivity of the chamber material in [ohm.m].
+        a1_L : float
+            Effective longitudinal radius of the chamber of the first power in
+            [m].
+        a2_L : float
+            Effective longitudinal radius of the chamber of the second power 
+            in [m].
+        a3_H : float
+            Effective horizontal radius of the chamber of the third power in
+            [m].
+        a4_H : float
+            Effective horizontal radius of the chamber of the fourth power in
+            [m].
+        a3_V : float
+            Effective vertical radius of the chamber of the third power in [m].
+        a4_V : float
+            Effective vertical radius of the chamber of the fourth power in 
+            [m].
+
+        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.         
+        """
+        
+        if x_right is True:
+            a0 = (self.x_right[1:] + self.x_right[:-1])/2
+        else:
+            a0 = np.abs((self.x_left[1:] + self.x_left[:-1])/2)
+            
+        if y_top is True:
+            b0 = (self.y_top[1:] + self.y_top[:-1])/2
+        else:
+            b0 = np.abs((self.y_bottom[1:] + self.y_bottom[:-1])/2)
+        
+        beta = optics.beta(self.center)
+        L = self.position[-1]
         sigma = 1/self.rho
-        beta_H_star = 1/L*(self.length*self.beta[0,:]).sum()
-        beta_V_star = 1/L*(self.length*self.beta[1,:]).sum()
+        beta_H_star = 1/L*(self.length*beta[0,:]).sum()
+        beta_V_star = 1/L*(self.length*beta[1,:]).sum()
         sigma_star = 1/L*(self.length*sigma).sum()
         
-        a1_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**1)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/1)
-        a2_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**2)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/2)
+        a1_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**1)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/1)
+        a2_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**2)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/2)
 
-        a1_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.b0)**1)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/1)
-        a2_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.b0)**2)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/2)
+        a1_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**1)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/1)
+        a2_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**2)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/2)
         
-        a3_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**3)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/3)
-        a4_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**4)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/4)
+        a3_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**3)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/3)
+        a4_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**4)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/4)
         
-        a3_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.b0)**3)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/3)
-        a4_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.b0)**4)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/4)
+        a3_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**3)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/3)
+        a4_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**4)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/4)
         
         a1_L = min((a1_H,a1_V))
         a2_L = min((a2_H,a2_V))
         
-        print("rhorw = " + str(1/sigma_star))
-        print("beffL1 = " + str(a1_L))
-        print("beffL2 = " + str(a2_L))
-        print("beffH3 = " + str(a3_H))
-        print("beffH4 = " + str(a4_H))
-        print("beffV3 = " + str(a3_V))
-        print("beffV4 = " + str(a4_V))
-        
-    def change_values(self, pos1, pos2, a0, b0, shape, rho):
-        ind = (self.position > pos1) & (self.position < pos2)
-        self.x0[ind] = a0
-        self.y0[ind] = b0
-        ind2 = (self.position[:-1] > pos1) & (self.position[1:] < pos2)
+        return (1/sigma_star, a1_L, a2_L, a3_H, a4_H, a3_V, a4_V)
+        
+    def change_values(self, start_position, end_position, x_right, y_top, 
+                      shape, rho, x_left=None, y_bottom=None):
+        """
+        Change the physical parameters between start_position and end_position.
+
+        Parameters
+        ----------
+        start_position : float
+        end_position : float
+        x_right : float
+            Right horizontal aperture in [m].
+        y_top : float
+            Top vertical aperture in [m].
+        shape : str
+            Shape of the chamber cross section (elli/circ/rect).
+        rho : float
+            Resistivity of the chamber material [ohm.m].
+        x_left : float, optional
+            Left horizontal aperture in [m].
+        y_bottom : float, optional
+            Bottom vertical aperture in [m].
+        """
+        ind = (self.position > start_position) & (self.position < end_position)
+        self.x_right[ind] = x_right
+        self.y_top[ind] = y_top
+        
+        if x_left is None:
+            self.x_left[ind] = -1*x_right
+        else:
+            self.x_left[ind] = x_left
+        
+        if y_bottom is None:
+            self.y_bottom[ind] = -1*y_top
+        else:
+            self.y_bottom[ind] = y_bottom
+        
+        ind2 = ((self.position[:-1] > start_position) & 
+                (self.position[1:] < end_position))
         self.rho[ind2] = rho
         self.shape[ind2] = shape
         
-    def taper(self, pos1, pos2, x0_start, y0_start, x0_end, y0_end, shape, rho):
-        ind = (self.position > pos1) & (self.position < pos2)
-        self.x0[ind] = np.linspace(x0_start, x0_end, sum(ind))
-        self.y0[ind] = np.linspace(y0_start, y0_end, sum(ind))
-        ind2 = (self.position[:-1] > pos1) & (self.position[1:] < pos2)
+    def taper(self, start_position, end_position, x_right_start, x_right_end,
+              y_top_start, y_top_end, shape, rho, x_left_start=None, 
+              x_left_end=None, y_bottom_start=None, y_bottom_end=None):
+        """
+        Change the physical parameters to have a tapered transition between 
+        start_position and end_position.
+
+        Parameters
+        ----------
+        start_position : float
+        end_position : float
+        x_right_start : float
+            Right horizontal aperture at the start of the taper in [m].
+        x_right_end : float
+            Right horizontal aperture at the end of the taper in [m].
+        y_top_start : float
+            Top vertical aperture at the start of the taper in [m].
+        y_top_end : float
+            Top vertical aperture at the end of the taper in [m].
+        shape : str
+            Shape of the chamber cross section (elli/circ/rect).
+        rho : float
+            Resistivity of the chamber material [ohm.m].
+        x_left_start : float, optional
+            Left horizontal aperture at the start of the taper in [m].
+        x_left_end : float, optional
+            Left horizontal aperture at the end of the taper in [m].
+        y_bottom_start : float, optional
+            Bottom vertical aperture at the start of the taper in [m].
+        y_bottom_end : float, optional
+            Bottom vertical aperture at the end of the taper in [m].
+        """
+        ind = (self.position > start_position) & (self.position < end_position)
+        self.x_right[ind] = np.linspace(x_right_start, x_right_end, sum(ind))
+        self.y_top[ind] = np.linspace(y_top_start, y_top_end, sum(ind))
+        
+        if (x_left_start is not None) and (x_left_end is not None):
+            self.x_left[ind] = np.linspace(x_left_start, x_left_end, sum(ind))
+        else:
+            self.x_left[ind] = -1*np.linspace(x_right_start, x_right_end, 
+                                              sum(ind))
+            
+        if (y_bottom_start is not None) and (y_bottom_end is not None):
+            self.y_bottom[ind] = np.linspace(y_bottom_start, y_bottom_end, 
+                                             sum(ind))
+        else:
+            self.y_bottom[ind] = -1*np.linspace(y_top_start, y_top_end, 
+                                                sum(ind))
+        
+        ind2 = ((self.position[:-1] > start_position) 
+                & (self.position[1:] < end_position))
         self.rho[ind2] = rho
         self.shape[ind2] = shape
         
     def plot_aperture(self):
-        plt.plot(self.position,self.x0*1e3)
-        plt.plot(self.position,self.y0*1e3)
-        plt.xlabel("Longitudinal position [m]")
-        plt.ylabel("H/V aperture [mm]")
-        plt.legend(["H","V"])
+        """Plot horizontal and vertical apertures."""
+        fig, axs = plt.subplots(2)
+        axs[0].plot(self.position,self.x_right*1e3)
+        axs[0].plot(self.position,self.x_left*1e3)
+        axs[0].set(xlabel="Longitudinal position [m]", 
+                   ylabel="Horizontal aperture [mm]")
+        axs[0].legend(["Right","Left"])
+        
+        axs[1].plot(self.position,self.y_top*1e3)
+        axs[1].plot(self.position,self.y_bottom*1e3)
+        axs[1].set(xlabel="Longitudinal position [m]", 
+                   ylabel="Vertical aperture [mm]")
+        axs[1].legend(["Top","Bottom"])
diff --git a/tracking/particles.py b/tracking/particles.py
index 4b1b783d8da6e84a7cdd1905390f4d255d396a17..94d525a47503c5bcb10e3276491bf01952b1d16c 100644
--- a/tracking/particles.py
+++ b/tracking/particles.py
@@ -12,7 +12,8 @@ from tracking.parallel import Mpi
 from scipy.constants import c, m_e, m_p, e
 
 class Particle:
-    """ Define a particle object
+    """
+    Define a particle object
 
     Attributes
     ----------
diff --git a/tracking/synchrotron.py b/tracking/synchrotron.py
index b32c3bc64a5eff43ac96f57ca1272e3fcb6c76cd..038248540226a1ad46854473f32f2a10665de7c0 100644
--- a/tracking/synchrotron.py
+++ b/tracking/synchrotron.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 """
-General elements
+Module where the synchrotron class is defined.
 
 @author: Alexis Gamelin
 @date: 15/01/2020
@@ -8,11 +8,66 @@ General elements
 
 import numpy as np
 from scipy.constants import c, e
-
         
 class Synchrotron:
-    """ Synchrotron class to store main properties. """
+    """
+    Synchrotron class to store main properties.
+    
+    Optional parameters are optional only if the Optics object passed to the 
+    class uses a loaded lattice.
     
+    Parameters
+    ----------
+    h : int
+        Harmonic number of the accelerator.
+    optics : Optics object
+        Object where the optic functions are stored.
+    particle : Particle object
+        Particle which is accelerated.
+    tau : array of shape (3,)
+        Horizontal, vertical and longitudinal damping times in [s].
+    sigma_delta : float
+        Equilibrium energy spread.
+    sigma_0 : float
+        Natural bunch length in [s].
+    emit : array of shape (2,)
+        Horizontal and vertical equilibrium emittance in [m.rad].
+    L : float, optional
+        Ring circumference in [m].
+    E0 : float, optional
+        Nominal (total) energy of the ring in [eV].
+    ac : float, optional
+        Momentum compaction factor.
+    tune : array of shape (2,), optional
+        Horizontal and vertical tunes.
+    chro : array of shape (2,), optional
+        Horizontal and vertical (non-normalized) chromaticities.
+    U0 : float, optional
+        Energy loss per turn in [eV].
+        
+    Attributes
+    ----------
+    T0 : float
+        Revolution time in [s].
+    f0 : float
+        Revolution frequency in [Hz].
+    omega0 : float
+        Angular revolution frequency in [Hz.rad]
+    f1 : float
+        Fundamental RF frequency in [Hz].
+    omega1 : float
+        Fundamental RF angular frequency in [Hz.rad].
+    k1 : float
+        Fundamental RF wave number in [m**-1].
+    gamma : float
+        Relativistic Lorentz gamma.
+    beta : float
+        Relativistic Lorentz beta.
+    eta : float
+        Momentum compaction.
+    sigma : array of shape (4,)
+        RMS beam size at equilibrium in [m].
+    """
     def __init__(self, h, optics, particle, **kwargs):
         self._h = h
         self.particle = particle
@@ -36,9 +91,7 @@ class Synchrotron:
         self.tau = kwargs.get('tau') # X/Y/S damping times [s]
         self.sigma_delta = kwargs.get('sigma_delta') # Equilibrium energy spread
         self.sigma_0 = kwargs.get('sigma_0') # Natural bunch length [s]
-        self.tune = kwargs.get('tune') # X/Y/S tunes
         self.emit = kwargs.get('emit') # X/Y emittances in [m.rad]
-        self.chro = kwargs.get('chro') # X/Y (non-normalized) chromaticities
                 
     @property
     def h(self):