diff --git a/gprMax/exceptions.py b/gprMax/exceptions.py
index 4f502289..75a34d9c 100644
--- a/gprMax/exceptions.py
+++ b/gprMax/exceptions.py
@@ -18,7 +18,9 @@
 
 import sys
 
-from colorama import init, Fore
+from colorama import init
+from colorama import Fore
+
 init()
 
 sys.tracebacklimit = None
diff --git a/gprMax/fields_updates.pyx b/gprMax/fields_updates.pyx
index 4af4c94c..5ff58514 100644
--- a/gprMax/fields_updates.pyx
+++ b/gprMax/fields_updates.pyx
@@ -20,7 +20,8 @@ import numpy as np
 cimport numpy as np
 from cython.parallel import prange
 
-from gprMax.constants cimport floattype_t, complextype_t
+from gprMax.constants cimport floattype_t
+from gprMax.constants cimport complextype_t
 
 
 ###############################################
@@ -28,13 +29,13 @@ from gprMax.constants cimport floattype_t, complextype_t
 ###############################################
 cpdef void update_electric(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
     """This function updates the electric field components.
-        
+
     Args:
         nx, ny, nz (int): Grid size in cells
         nthreads (int): Number of threads to use
         updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int materialEx, materialEy, materialEz
 
@@ -98,14 +99,14 @@ cpdef void update_electric(int nx, int ny, int nz, int nthreads, floattype_t[:,
 #################################################
 cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
     """This function updates the electric field components when dispersive materials (with multiple poles) are present.
-        
+
     Args:
         nx, ny, nz (int): Grid size in cells
         maxpoles (int): Maximum number of poles
         nthreads (int): Number of threads to use
         updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
     """
-    
+
     cdef Py_ssize_t i, j, k, pole
     cdef int material
     cdef float phi = 0
@@ -150,14 +151,14 @@ cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nt
 
 cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez):
     """This function updates a temporary dispersive material array when disperisive materials (with multiple poles) are present.
-        
+
     Args:
         nx, ny, nz (int): Grid size in cells
         maxpoles (int): Maximum number of poles
         nthreads (int): Number of threads to use
         updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
     """
-    
+
     cdef Py_ssize_t i, j, k, pole
     cdef int material
 
@@ -191,13 +192,13 @@ cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nt
 
 cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
     """This function updates the electric field components when dispersive materials (with 1 pole) are present.
-        
+
     Args:
         nx, ny, nz (int): Grid size in cells
         nthreads (int): Number of threads to use
         updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int material
     cdef float phi = 0
@@ -235,13 +236,13 @@ cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthrea
 
 cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez):
     """This function updates a temporary dispersive material array when disperisive materials (with 1 pole) are present.
-        
+
     Args:
         nx, ny, nz (int): Grid size in cells
         nthreads (int): Number of threads to use
         updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int material
 
@@ -275,13 +276,13 @@ cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthrea
 ##########################
 cpdef void update_magnetic(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
     """This function updates the magnetic field components.
-        
+
     Args:
         nx, ny, nz (int): Grid size in cells
         nthreads (int): Number of threads to use
         updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int materialHx, materialHy, materialHz
 
diff --git a/gprMax/fractals.py b/gprMax/fractals.py
index 3eca2f3e..d4123bd8 100644
--- a/gprMax/fractals.py
+++ b/gprMax/fractals.py
@@ -18,7 +18,8 @@
 
 import numpy as np
 
-from gprMax.constants import floattype, complextype
+from gprMax.constants import floattype
+from gprMax.constants import complextype
 from gprMax.utilities import round_value
 
 np.seterr(divide='raise')
diff --git a/gprMax/geometry_primitives.pyx b/gprMax/geometry_primitives.pyx
index ac2a5659..6f2cc73f 100644
--- a/gprMax/geometry_primitives.pyx
+++ b/gprMax/geometry_primitives.pyx
@@ -20,27 +20,36 @@ import numpy as np
 cimport numpy as np
 
 from gprMax.utilities import round_value
-from gprMax.yee_cell_setget_rigid cimport set_rigid_Ex, set_rigid_Ey, set_rigid_Ez, set_rigid_Hx, set_rigid_Hy, set_rigid_Hz, set_rigid_E, unset_rigid_E, set_rigid_H, unset_rigid_H
+from gprMax.yee_cell_setget_rigid cimport set_rigid_Ex
+from gprMax.yee_cell_setget_rigid cimport set_rigid_Ey
+from gprMax.yee_cell_setget_rigid cimport set_rigid_Ez
+from gprMax.yee_cell_setget_rigid cimport set_rigid_Hx
+from gprMax.yee_cell_setget_rigid cimport set_rigid_Hy
+from gprMax.yee_cell_setget_rigid cimport set_rigid_Hz
+from gprMax.yee_cell_setget_rigid cimport set_rigid_E
+from gprMax.yee_cell_setget_rigid cimport unset_rigid_E
+from gprMax.yee_cell_setget_rigid cimport set_rigid_H
+from gprMax.yee_cell_setget_rigid cimport unset_rigid_H
 
 np.seterr(divide='raise')
 
 
 cpdef bint are_clockwise(float v1x, float v1y, float v2x, float v2y):
     """Find if vector 2 is clockwise relative to vector 1.
-        
+
     Args:
         v1x, v1y, v2x, v2y (float): Coordinates of vectors.
 
     Returns:
         (boolean)
     """
-    
+
     return -v1x*v2y + v1y*v2x > 0
 
 
 cpdef bint is_within_radius(float vx, float vy, float radius):
     """Check if the point is within a given radius of the centre of the circle.
-        
+
     Args:
         vx, vy (float): Coordinates of vector.
         radius (float): Radius.
@@ -48,7 +57,7 @@ cpdef bint is_within_radius(float vx, float vy, float radius):
     Returns:
         (boolean)
     """
-    
+
     return vx*vx + vy*vy <= radius*radius
 
 
@@ -59,7 +68,7 @@ cpdef bint is_inside_sector(float px, float py, float ctrx, float ctry, float se
         It has to be closer to the center of the circle than the sectors radius.
         Assumes sector start is always clockwise from sector end,
         i.e. sector defined in an anti-clockwise direction
-        
+
     Args:
         px, py (float): Coordinates of point.
         ctrx, ctry (float): Coordinates of centre of circle.
@@ -70,26 +79,26 @@ cpdef bint is_inside_sector(float px, float py, float ctrx, float ctry, float se
     Returns:
         (boolean)
     """
-    
+
     cdef float sectorstart1, sectorstart2, sectorend1, sectorend2, relpoint1, relpoint2
-    
+
     sectorstart1 = radius * np.cos(sectorstartangle)
     sectorstart2 = radius * np.sin(sectorstartangle)
     sectorend1 = radius * np.cos(sectorstartangle + sectorangle)
     sectorend2 = radius * np.sin(sectorstartangle + sectorangle)
     relpoint1 = px - ctrx
     relpoint2 = py - ctry
-    
+
     return not are_clockwise(sectorstart1, sectorstart2, relpoint1, relpoint2) and are_clockwise(sectorend1, sectorend2, relpoint1, relpoint2) and is_within_radius(relpoint1, relpoint2, radius)
 
 
 cpdef bint point_in_polygon(float px, float py, list polycoords):
     """Calculates, using a ray casting algorithm, whether a point lies within a polygon.
-            
+
     Args:
         px, py (float): Coordinates of point to test.
         polycoords (list): x, y tuples of coordinates that define the polygon.
-        
+
     Returns:
         inside (boolean)
     """
@@ -114,7 +123,7 @@ cpdef bint point_in_polygon(float px, float py, list polycoords):
             p2x, p2y = polycoords[i]
         if p1y == p2y and p1y == py and px > min(p1x, p2x) and px < max(p1x, p2x):
             return True
-      
+
     inside = False
 
     p1x, p1y = polycoords[0]
@@ -134,52 +143,52 @@ cpdef bint point_in_polygon(float px, float py, list polycoords):
 
 cpdef void build_edge_x(int i, int j, int k, int numIDx, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Set x-orientated edges in the rigid and ID arrays for a Yee voxel.
-        
+
     Args:
         i, j, k (int): Cell coordinates of edge.
         numIDz (int): Numeric ID of material.
         rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
     """
-    
+
     set_rigid_Ex(i, j, k, rigidE)
     ID[0, i, j, k] = numIDx
 
 
 cpdef void build_edge_y(int i, int j, int k, int numIDy, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Set y-orientated edges in the rigid and ID arrays for a Yee voxel.
-        
+
     Args:
         i, j, k (int): Cell coordinates of edge.
         numIDz (int): Numeric ID of material.
         rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
     """
-    
+
     set_rigid_Ey(i, j, k, rigidE)
     ID[1, i, j, k] = numIDy
 
 
 cpdef void build_edge_z(int i, int j, int k, int numIDz, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Set z-orientated edges in the rigid and ID arrays for a Yee voxel.
-        
+
     Args:
         i, j, k (int): Cell coordinates of edge.
         numIDz (int): Numeric ID of material.
         rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
     """
-    
+
     set_rigid_Ez(i, j, k, rigidE)
     ID[2, i, j, k] = numIDz
 
 
 cpdef void build_face_yz(int i, int j, int k, int numIDy, int numIDz, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Set the edges of the yz-plane face of a Yell cell in the rigid and ID arrays.
-        
+
     Args:
         i, j, k (int): Cell coordinates of the face.
         numIDx, numIDy (int): Numeric ID of material.
         rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
     """
-    
+
     set_rigid_Ey(i, j, k, rigidE)
     set_rigid_Ez(i, j, k, rigidE)
     set_rigid_Ey(i, j, k + 1, rigidE)
@@ -200,13 +209,13 @@ cpdef void build_face_yz(int i, int j, int k, int numIDy, int numIDz, np.int8_t[
 
 cpdef void build_face_xz(int i, int j, int k, int numIDx, int numIDz, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Set the edges of the xz-plane face of a Yell cell in the rigid and ID arrays.
-        
+
     Args:
         i, j, k (int): Cell coordinates of the face.
         numIDx, numIDy (int): Numeric ID of material.
         rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
     """
-    
+
     set_rigid_Ex(i, j, k, rigidE)
     set_rigid_Ez(i, j, k, rigidE)
     set_rigid_Ex(i, j, k + 1, rigidE)
@@ -227,13 +236,13 @@ cpdef void build_face_xz(int i, int j, int k, int numIDx, int numIDz, np.int8_t[
 
 cpdef void build_face_xy(int i, int j, int k, int numIDx, int numIDy, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Set the edges of the xy-plane face of a Yell cell in the rigid and ID arrays.
-        
+
     Args:
         i, j, k (int): Cell coordinates of the face.
         numIDx, numIDy (int): Numeric ID of material.
         rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
     """
-    
+
     set_rigid_Ex(i, j, k, rigidE)
     set_rigid_Ey(i, j, k, rigidE)
     set_rigid_Ex(i, j + 1, k, rigidE)
@@ -254,14 +263,14 @@ cpdef void build_face_xy(int i, int j, int k, int numIDx, int numIDy, np.int8_t[
 
 cpdef void build_voxel(int i, int j, int k, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Set values in the solid, rigid and ID arrays for a Yee voxel.
-        
+
     Args:
         i, j, k (int): Cell coordinates of voxel.
         numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
         averaging (bint): Whether material property averaging will occur for the object.
         solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
     """
-    
+
     if averaging:
         solid[i, j, k] = numID
         unset_rigid_E(i, j, k, rigidE)
@@ -276,12 +285,12 @@ cpdef void build_voxel(int i, int j, int k, int numID, int numIDx, int numIDy, i
         ID[0, i, j + 1, k + 1] = numIDx
         ID[0, i, j + 1, k] = numIDx
         ID[0, i, j, k + 1] = numIDx
-        
+
         ID[1, i, j, k] = numIDy
         ID[1, i + 1, j, k + 1] = numIDy
         ID[1, i + 1, j, k] = numIDy
         ID[1, i, j, k + 1] = numIDy
-        
+
         ID[2, i, j, k] = numIDz
         ID[2, i + 1, j + 1, k] = numIDz
         ID[2, i + 1, j, k] = numIDz
@@ -291,12 +300,12 @@ cpdef void build_voxel(int i, int j, int k, int numID, int numIDx, int numIDy, i
         ID[3, i, j + 1, k + 1] = numIDx
         ID[3, i, j + 1, k] = numIDx
         ID[3, i, j, k + 1] = numIDx
-        
+
         ID[4, i, j, k] = numIDy
         ID[4, i + 1, j, k + 1] = numIDy
         ID[4, i + 1, j, k] = numIDy
         ID[4, i, j, k + 1] = numIDy
-        
+
         ID[5, i, j, k] = numIDz
         ID[5, i + 1, j + 1, k] = numIDz
         ID[5, i + 1, j, k] = numIDz
@@ -305,7 +314,7 @@ cpdef void build_voxel(int i, int j, int k, int numID, int numIDx, int numIDy, i
 
 cpdef void build_triangle(float x1, float y1, float z1, float x2, float y2, float z2, float x3, float y3, float z3, str normal, float thickness, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Builds #triangle and #triangular_prism commands which sets values in the solid, rigid and ID arrays for a Yee voxel.
-        
+
     Args:
         x1, y1, z1, x2, y2, z2, x3, y3, z3 (float): Coordinates of the vertices of the triangular prism.
         normal (char): Normal direction to the plane of the triangular prism.
@@ -315,11 +324,11 @@ cpdef void build_triangle(float x1, float y1, float z1, float x2, float y2, floa
         averaging (bint): Whether material property averaging will occur for the object.
         solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int i1, i2, j1, j2, sign, level, thicknesscells
     cdef float area, s, t
-    
+
     # Calculate a bounding box for the triangle
     if normal == 'x':
         area = 0.5 * (-z2 * y3 + z1 * (-y2 + y3) + y1 * (z2 - z3) + y2 * z3)
@@ -350,7 +359,7 @@ cpdef void build_triangle(float x1, float y1, float z1, float x2, float y2, floa
 
     for i in range(i1, i2):
         for j in range(j1, j2):
-            
+
             # Calculate the areas of the 3 triangles defined by the 3 vertices of the main triangle and the point under test
             if normal == 'x':
                 ir = (i + 0.5) * dy
@@ -389,7 +398,7 @@ cpdef void build_triangle(float x1, float y1, float z1, float x2, float y2, floa
 
 cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sectorstartangle, float sectorangle, float radius, str normal, float thickness, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Builds #cylindrical_sector commands which sets values in the solid, rigid and ID arrays for a Yee voxel. It defines a sector of cylinder given by the direction of the axis of the coordinates of the cylinder face centre, depth coordinates, sector start point, sector angle, and sector radius. N.B Assumes sector start is always clockwise from sector end, i.e. sector defined in an anti-clockwise direction.
-        
+
     Args:
         ctr1, ctr2 (float): Coordinates of centre of circle.
         level (int): Third dimensional coordinate.
@@ -403,10 +412,10 @@ cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sec
         averaging (bint): Whether material property averaging will occur for the object.
         solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
     """
-    
+
     cdef Py_ssize_t x, y, z
     cdef int x1, x2, y1, y2, z1, z2, thicknesscells
-    
+
     if normal == 'x':
         # Angles are defined from zero degrees on the positive y-axis going towards positive z-axis
         y1 = round_value((ctr1 - radius)/dy)
@@ -414,7 +423,7 @@ cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sec
         z1 = round_value((ctr2 - radius)/dz)
         z2 = round_value((ctr2 + radius)/dz)
         thicknesscells = round_value(thickness/dx)
-        
+
         for y in range(y1, y2):
             for z in range(z1, z2):
                 if is_inside_sector(y * dy + 0.5 * dy, z * dz + 0.5 * dz, ctr1, ctr2, sectorstartangle, sectorangle, radius):
@@ -431,7 +440,7 @@ cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sec
         z1 = round_value((ctr2 - radius)/dz)
         z2 = round_value((ctr2 + radius)/dz)
         thicknesscells = round_value(thickness/dy)
-            
+
         for x in range(x1, x2):
             for z in range(z1, z2):
                 if is_inside_sector(x * dx + 0.5 * dx, z * dz + 0.5 * dz, ctr1, ctr2, sectorstartangle, sectorangle, radius):
@@ -448,7 +457,7 @@ cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sec
         y1 = round_value((ctr2 - radius)/dy)
         y2 = round_value((ctr2 + radius)/dy)
         thicknesscells = round_value(thickness/dz)
-            
+
         for x in range(x1, x2):
             for y in range(y1, y2):
                 if is_inside_sector(x * dx + 0.5 * dx, y * dy + 0.5 * dy, ctr1, ctr2, sectorstartangle, sectorangle, radius):
@@ -461,14 +470,14 @@ cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sec
 
 cpdef void build_box(int xs, int xf, int ys, int yf, int zs, int zf, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Builds #box commands which sets values in the solid, rigid and ID arrays.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box.
         numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
         averaging (bint): Whether material property averaging will occur for the object.
         solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
     """
-    
+
     cdef Py_ssize_t i, j, k
 
     if averaging:
@@ -505,12 +514,12 @@ cpdef void build_box(int xs, int xf, int ys, int yf, int zs, int zf, int numID,
             for j in range(ys, yf):
                 for k in range(zs, zf):
                     ID[3, i, j, k] = numIDx
-                        
+
         for i in range(xs, xf):
             for j in range(ys, yf + 1):
                 for k in range(zs, zf):
                     ID[4, i, j, k] = numIDy
-                        
+
         for i in range(xs, xf):
             for j in range(ys, yf):
                 for k in range(zs, zf + 1):
@@ -519,7 +528,7 @@ cpdef void build_box(int xs, int xf, int ys, int yf, int zs, int zf, int numID,
 
 cpdef void build_cylinder(float x1, float y1, float z1, float x2, float y2, float z2, float r, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Builds #cylinder commands which sets values in the solid, rigid and ID arrays for a Yee voxel.
-        
+
     Args:
         x1, y1, z1, x2, y2, z2 (float): Coordinates of the centres of cylinder faces.
         r (float): Radius of the cylinder.
@@ -528,7 +537,7 @@ cpdef void build_cylinder(float x1, float y1, float z1, float x2, float y2, floa
         averaging (bint): Whether material property averaging will occur for the object.
         solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int xs, xf, ys, yf, zs, zf
     cdef float f1f2mag, f2f1mag, f1ptmag, f2ptmag, dot1, dot2, factor1, factor2, theta1, theta2, distance1, distance2
@@ -611,14 +620,14 @@ cpdef void build_cylinder(float x1, float y1, float z1, float x2, float y2, floa
                     distance2 = f2ptmag * np.sin(theta2)
                     if (distance1 <= r or distance2 <= r) and theta1 <= np.pi/2 and theta2 <= np.pi/2:
                         build = 1
-                
+
                 if build:
                     build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
 
 
 cpdef void build_sphere(int xc, int yc, int zc, float r, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Builds #sphere commands which sets values in the solid, rigid and ID arrays for a Yee voxel.
-        
+
     Args:
         xc, yc, zc (int): Cell coordinates of the centre of the sphere.
         r (float): Radius of the sphere.
@@ -627,10 +636,10 @@ cpdef void build_sphere(int xc, int yc, int zc, float r, float dx, float dy, flo
         averaging (bint): Whether material property averaging will occur for the object.
         solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int xs, xf, ys, yf, zs, zf
-    
+
     # Calculate a bounding box for sphere
     xs = round_value(((xc * dx) - r) / dx) - 1
     xf = round_value(((xc * dx) + r) / dx) + 1
@@ -638,7 +647,7 @@ cpdef void build_sphere(int xc, int yc, int zc, float r, float dx, float dy, flo
     yf = round_value(((yc * dy) + r) / dy) + 1
     zs = round_value(((zc * dz) - r) / dz) - 1
     zf = round_value(((zc * dz) + r) / dz) + 1
-    
+
     # Set bounds to domain if they outside
     if xs < 0:
         xs = 0
@@ -662,7 +671,7 @@ cpdef void build_sphere(int xc, int yc, int zc, float r, float dx, float dy, flo
 
 cpdef void build_voxels_from_array(int xs, int ys, int zs, int numexistmaterials, bint averaging, np.int16_t[:, :, ::1] data, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Builds Yee voxels by reading integers from an array.
-        
+
     Args:
         xs, ys, zs (int): Cell coordinates of position of start of array in domain.
         numexistmaterials (int): Number of existing materials in model prior to building voxels.
@@ -673,7 +682,7 @@ cpdef void build_voxels_from_array(int xs, int ys, int zs, int numexistmaterials
 
     cdef Py_ssize_t i, j, k
     cdef int xf, yf, zf, numID
-    
+
     # Set bounds to domain if they outside
     if xs < 0:
         xs = 0
@@ -707,7 +716,7 @@ cpdef void build_voxels_from_array(int xs, int ys, int zs, int numexistmaterials
 
 cpdef void build_voxels_from_array_mask(int xs, int ys, int zs, int waternumID, int grassnumID, bint averaging, np.int8_t[:, :, ::1] mask, np.int16_t[:, :, ::1] data, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID):
     """Builds Yee voxels by reading integers from an array.
-        
+
     Args:
         xs, ys, zs (int): Cell coordinates of position of start of array in domain.
         waternumID, grassnumID (int): Numeric ID of water and grass materials.
@@ -719,7 +728,7 @@ cpdef void build_voxels_from_array_mask(int xs, int ys, int zs, int waternumID,
 
     cdef Py_ssize_t i, j, k
     cdef int xf, yf, zf, numID, numIDx, numIDy, numIDz
-    
+
     # Set upper bounds
     xf = xs + data.shape[0]
     yf = ys + data.shape[1]
diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py
index 258464cc..464013bc 100644
--- a/gprMax/gprMax.py
+++ b/gprMax/gprMax.py
@@ -20,19 +20,27 @@
 
 import argparse
 import datetime
-from enum import Enum
 import os
-from time import perf_counter
 import sys
 
+from enum import Enum
+from time import perf_counter
+
 import h5py
 import numpy as np
 
 from gprMax._version import __version__, codename
-from gprMax.constants import c, e0, m0, z0
+from gprMax.constants import c
+from gprMax.constants import e0
+from gprMax.constants import m0
+from gprMax.constants import z0
 from gprMax.exceptions import GeneralError
 from gprMax.model_build_run import run_model
-from gprMax.utilities import get_host_info, get_terminal_width, human_size, logo, open_path_file
+from gprMax.utilities import get_host_info
+from gprMax.utilities import get_terminal_width
+from gprMax.utilities import human_size
+from gprMax.utilities import logo
+from gprMax.utilities import open_path_file
 
 
 def main():
@@ -90,7 +98,7 @@ def run_main(args):
     Args:
         args (dict): Namespace with input arguments from command line or api.
     """
-    
+
     with open_path_file(args.inputfile) as inputfile:
 
         # Get information about host machine
@@ -130,7 +138,7 @@ def run_main(args):
                 if args.task:
                     raise GeneralError('MPI cannot be combined with job array mode')
                 run_mpi_sim(args, inputfile, usernamespace)
-        
+
             # Standard behaviour - models run serially with each model parallelised with OpenMP (CPU) or CUDA (GPU)
             else:
                 if args.task and args.restart:
@@ -160,7 +168,7 @@ def run_std_sim(args, inputfile, usernamespace, optparams=None):
         modelstart = 1
         modelend = modelstart + args.n
     numbermodelruns = args.n
-    
+
     tsimstart = perf_counter()
     for currentmodelrun in range(modelstart, modelend):
         if optparams:  # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace
@@ -211,7 +219,7 @@ def run_benchmark_sim(args, inputfile, usernamespace):
 
     numbermodelruns = len(cputhreads)
     modelend = numbermodelruns + 1
-                
+
     usernamespace['number_model_runs'] = numbermodelruns
 
     for currentmodelrun in range(1, modelend):
@@ -246,7 +254,7 @@ def run_mpi_sim(args, inputfile, usernamespace, optparams=None):
     """
 
     from mpi4py import MPI
-    
+
     # Get name of processor/host
     name = MPI.Get_processor_name()
 
@@ -254,18 +262,18 @@ def run_mpi_sim(args, inputfile, usernamespace, optparams=None):
     modelstart = args.restart if args.restart else 1
     modelend = modelstart + args.n
     numbermodelruns = args.n
-    
+
     # Number of workers and command line flag to indicate a spawned worker
     worker = '--mpi-worker'
     numberworkers = args.mpi - 1
 
     # Master process
     if worker not in sys.argv:
-        
+
         tsimstart = perf_counter()
-        
+
         print('MPI master rank (PID {}) on {} using {} workers'.format(os.getpid(), name, numberworkers))
-        
+
         # Create a list of work
         worklist = []
         for model in range(modelstart, modelend):
@@ -306,10 +314,10 @@ def run_mpi_sim(args, inputfile, usernamespace, optparams=None):
         # Ask for work until stop sentinel
         for work in iter(lambda: comm.sendrecv(0, dest=0), StopIteration):
             currentmodelrun = work['currentmodelrun']
-            
+
             gpuinfo = ''
             print('MPI worker rank {} (PID {}) starting model {}/{}{} on {}'.format(rank, os.getpid(), currentmodelrun, numbermodelruns, gpuinfo, name))
-            
+
             # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace
             if 'optparams' in work:
                 tmp = {}
diff --git a/gprMax/grid.py b/gprMax/grid.py
index a9b9ef4c..df92b6d3 100644
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -18,12 +18,16 @@
 
 from collections import OrderedDict
 
-from colorama import init, Fore, Style
+from colorama import init
+from colorama import Fore
+from colorama import Style
 init()
 import numpy as np
 np.seterr(invalid='raise')
 
-from gprMax.constants import c, floattype, complextype
+from gprMax.constants import c
+from gprMax.constants import floattype
+from gprMax.constants import complextype
 from gprMax.materials import Material
 from gprMax.pml import PML
 from gprMax.utilities import round_value
diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py
index cfa48b30..4aef1091 100644
--- a/gprMax/input_cmds_geometry.py
+++ b/gprMax/input_cmds_geometry.py
@@ -27,10 +27,25 @@ from gprMax.constants import floattype
 from gprMax.input_cmds_file import check_cmd_names
 from gprMax.input_cmds_multiuse import process_multicmds
 from gprMax.exceptions import CmdInputError
-from gprMax.fractals import FractalSurface, FractalVolume, Grass
-from gprMax.geometry_primitives import build_edge_x, build_edge_y, build_edge_z, build_face_yz, build_face_xz, build_face_xy, build_triangle, build_box, build_cylinder, build_cylindrical_sector, build_sphere, build_voxels_from_array, build_voxels_from_array_mask
+from gprMax.fractals import FractalSurface
+from gprMax.fractals import FractalVolume
+from gprMax.fractals import Grass
+from gprMax.geometry_primitives import build_edge_x
+from gprMax.geometry_primitives import build_edge_y
+from gprMax.geometry_primitives import build_edge_z
+from gprMax.geometry_primitives import build_face_yz
+from gprMax.geometry_primitives import build_face_xz
+from gprMax.geometry_primitives import build_face_xy
+from gprMax.geometry_primitives import build_triangle
+from gprMax.geometry_primitives import build_box
+from gprMax.geometry_primitives import build_cylinder
+from gprMax.geometry_primitives import build_cylindrical_sector
+from gprMax.geometry_primitives import build_sphere
+from gprMax.geometry_primitives import build_voxels_from_array
+from gprMax.geometry_primitives import build_voxels_from_array_mask
 from gprMax.materials import Material
-from gprMax.utilities import round_value, get_terminal_width
+from gprMax.utilities import round_value
+from gprMax.utilities import get_terminal_width
 
 
 def process_geometrycmds(geometry, G):
diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py
index 92f2aec0..d4df19b9 100644
--- a/gprMax/input_cmds_multiuse.py
+++ b/gprMax/input_cmds_multiuse.py
@@ -16,19 +16,28 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-from colorama import init, Fore, Style
+from colorama import init
+from colorama import Fore
+from colorama import Style
 init()
 import numpy as np
 from tqdm import tqdm
 
-from gprMax.constants import z0, floattype
+from gprMax.constants import z0
+from gprMax.constants import floattype
 from gprMax.exceptions import CmdInputError
-from gprMax.geometry_outputs import GeometryView, GeometryObjects
-from gprMax.materials import Material, PeplinskiSoil
-from gprMax.pml import CFSParameter, CFS
+from gprMax.geometry_outputs import GeometryView
+from gprMax.geometry_outputs import GeometryObjects
+from gprMax.materials import Material
+from gprMax.materials import PeplinskiSoil
+from gprMax.pml import CFSParameter
+from gprMax.pml import CFS
 from gprMax.receivers import Rx
 from gprMax.snapshots import Snapshot
-from gprMax.sources import VoltageSource, HertzianDipole, MagneticDipole, TransmissionLine
+from gprMax.sources import VoltageSource
+from gprMax.sources import HertzianDipole
+from gprMax.sources import MagneticDipole
+from gprMax.sources import TransmissionLine
 from gprMax.utilities import round_value
 from gprMax.waveforms import Waveform
 
diff --git a/gprMax/input_cmds_singleuse.py b/gprMax/input_cmds_singleuse.py
index 747ae146..382dabab 100644
--- a/gprMax/input_cmds_singleuse.py
+++ b/gprMax/input_cmds_singleuse.py
@@ -24,9 +24,14 @@ from colorama import init, Fore, Style
 init()
 import numpy as np
 
-from gprMax.constants import c, floattype
-from gprMax.exceptions import CmdInputError, GeneralError
-from gprMax.utilities import get_host_info, human_size, memory_usage, round_value
+from gprMax.constants import c
+from gprMax.constants import floattype
+from gprMax.exceptions import CmdInputError
+from gprMax.exceptions import GeneralError
+from gprMax.utilities import get_host_info
+from gprMax.utilities import human_size
+from gprMax.utilities import memory_usage
+from gprMax.utilities import round_value
 from gprMax.waveforms import Waveform
 
 
diff --git a/gprMax/materials.py b/gprMax/materials.py
index 7f37e534..2232419b 100644
--- a/gprMax/materials.py
+++ b/gprMax/materials.py
@@ -18,7 +18,9 @@
 
 import numpy as np
 
-from gprMax.constants import e0, m0, complextype
+from gprMax.constants import e0
+from gprMax.constants import m0
+from gprMax.constants import complextype
 
 
 class Material(object):
diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py
index e48309ad..f8bb3069 100644
--- a/gprMax/model_build_run.py
+++ b/gprMax/model_build_run.py
@@ -31,17 +31,33 @@ from tqdm import tqdm
 
 from gprMax.constants import floattype, cfloattype, ccomplextype
 from gprMax.exceptions import GeneralError
-from gprMax.fields_outputs import store_outputs, write_hdf5_outputfile
-from gprMax.fields_updates import update_electric, update_magnetic, update_electric_dispersive_multipole_A, update_electric_dispersive_multipole_B, update_electric_dispersive_1pole_A, update_electric_dispersive_1pole_B
-from gprMax.grid import FDTDGrid, dispersion_analysis
+
+from gprMax.fields_outputs import store_outputs
+from gprMax.fields_outputs import write_hdf5_outputfile
+
+from gprMax.fields_updates import update_electric
+from gprMax.fields_updates import update_magnetic
+from gprMax.fields_updates import update_electric_dispersive_multipole_A
+from gprMax.fields_updates import update_electric_dispersive_multipole_B
+from gprMax.fields_updates import update_electric_dispersive_1pole_A
+from gprMax.fields_updates import update_electric_dispersive_1pole_B
+
+from gprMax.grid import FDTDGrid
+from gprMax.grid import dispersion_analysis
 from gprMax.input_cmds_geometry import process_geometrycmds
-from gprMax.input_cmds_file import process_python_include_code, write_processed_file, check_cmd_names
+from gprMax.input_cmds_file import process_python_include_code
+from gprMax.input_cmds_file import write_processed_file
+from gprMax.input_cmds_file import check_cmd_names
 from gprMax.input_cmds_multiuse import process_multicmds
 from gprMax.input_cmds_singleuse import process_singlecmds
 from gprMax.materials import Material, process_materials
-from gprMax.pml import PML, build_pmls
-from gprMax.utilities import get_terminal_width, human_size, open_path_file
-from gprMax.yee_cell_build import build_electric_components, build_magnetic_components
+from gprMax.pml import PML
+from gprMax.pml import build_pmls
+from gprMax.utilities import get_terminal_width
+from gprMax.utilities import human_size
+from gprMax.utilities import open_path_file
+from gprMax.yee_cell_build import build_electric_components
+from gprMax.yee_cell_build import build_magnetic_components
 
 
 def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usernamespace):
@@ -67,7 +83,7 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
 
     # Used for naming geometry and output files
     appendmodelnumber = '' if numbermodelruns == 1 and not args.task and not args.restart else str(currentmodelrun)
-    
+
     # Normal model reading/building process; bypassed if geometry information to be reused
     if 'G' not in globals():
 
diff --git a/gprMax/optimisation_taguchi.py b/gprMax/optimisation_taguchi.py
index 474bbd28..29ecf8e9 100644
--- a/gprMax/optimisation_taguchi.py
+++ b/gprMax/optimisation_taguchi.py
@@ -30,8 +30,10 @@ import numpy as np
 
 from gprMax.constants import floattype
 from gprMax.exceptions import CmdInputError
-from gprMax.gprMax import run_std_sim, run_mpi_sim
-from gprMax.utilities import get_terminal_width, open_path_file
+from gprMax.gprMax import run_std_sim
+from gprMax.gprMax import run_mpi_sim
+from gprMax.utilities import get_terminal_width
+from gprMax.utilities import  open_path_file
 
 
 def run_opt_sim(args, inputfile, usernamespace):
@@ -190,7 +192,7 @@ def taguchi_code_blocks(inputfile, taguchinamespace):
 
     # Strip out any newline characters and comments that must begin with double hashes
     inputlines = [line.rstrip() for line in inputfile if(not line.startswith('##') and line.rstrip('\n'))]
-    
+
     # Rewind input file in preparation for passing to standard command reading function
     inputfile.seek(0)
 
diff --git a/gprMax/pml.py b/gprMax/pml.py
index 95ec1107..974df3db 100644
--- a/gprMax/pml.py
+++ b/gprMax/pml.py
@@ -21,7 +21,9 @@ from importlib import import_module
 import numpy as np
 from tqdm import tqdm
 
-from gprMax.constants import e0, z0, floattype
+from gprMax.constants import e0
+from gprMax.constants import z0
+from gprMax.constants import floattype
 
 
 class CFSParameter(object):
diff --git a/gprMax/pml_1order_update.pyx b/gprMax/pml_1order_update.pyx
index 231acef6..f81d06d0 100644
--- a/gprMax/pml_1order_update.pyx
+++ b/gprMax/pml_1order_update.pyx
@@ -20,7 +20,8 @@ import numpy as np
 cimport numpy as np
 from cython.parallel import prange
 
-from gprMax.constants cimport floattype_t, complextype_t
+from gprMax.constants cimport floattype_t
+from gprMax.constants cimport complextype_t
 
 
 ############################################
@@ -28,7 +29,7 @@ from gprMax.constants cimport floattype_t, complextype_t
 ############################################
 cpdef void update_pml_1order_electric_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ey and Ez field components for the xminus slab.
-        
+
         Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -36,7 +37,7 @@ cpdef void update_pml_1order_electric_xminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEy, materialEz
     cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0
@@ -44,7 +45,7 @@ cpdef void update_pml_1order_electric_xminus(int xs, int xf, int ys, int yf, int
     nx = xf - xs
     ny = yf - ys
     nz = zf - zs
-    
+
     for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
         RA0 = (ERA[0, i] - 1)
         RB0 = ERB[0, i]
@@ -72,7 +73,7 @@ cpdef void update_pml_1order_electric_xminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_1order_electric_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ey and Ez field components for the xplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -80,7 +81,7 @@ cpdef void update_pml_1order_electric_xplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEy, materialEz
     cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0
@@ -116,7 +117,7 @@ cpdef void update_pml_1order_electric_xplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_1order_electric_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ez field components for the yminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -124,7 +125,7 @@ cpdef void update_pml_1order_electric_yminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEz
     cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0
@@ -160,7 +161,7 @@ cpdef void update_pml_1order_electric_yminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_1order_electric_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ez field components for the yplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -168,7 +169,7 @@ cpdef void update_pml_1order_electric_yplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEz
     cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0
@@ -204,7 +205,7 @@ cpdef void update_pml_1order_electric_yplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_1order_electric_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ey field components for the zminus slab.
-        
+
         Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -212,7 +213,7 @@ cpdef void update_pml_1order_electric_zminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
         """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEy
     cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0
@@ -220,7 +221,7 @@ cpdef void update_pml_1order_electric_zminus(int xs, int xf, int ys, int yf, int
     nx = xf - xs
     ny = yf - ys
     nz = zf - zs
-    
+
     for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
         ii = i + xs
         for j in range(0, ny):
@@ -248,7 +249,7 @@ cpdef void update_pml_1order_electric_zminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_1order_electric_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ey field components for the zplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -256,7 +257,7 @@ cpdef void update_pml_1order_electric_zplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEy
     cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0
@@ -292,7 +293,7 @@ cpdef void update_pml_1order_electric_zplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_1order_magnetic_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hy and Hz field components for the xminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -300,7 +301,7 @@ cpdef void update_pml_1order_magnetic_xminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, ERE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHy, materialHz
     cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0
@@ -336,7 +337,7 @@ cpdef void update_pml_1order_magnetic_xminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_1order_magnetic_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hy and Hz field components for the xplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -344,7 +345,7 @@ cpdef void update_pml_1order_magnetic_xplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHy, materialHz
     cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0
@@ -380,7 +381,7 @@ cpdef void update_pml_1order_magnetic_xplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_1order_magnetic_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hz field components for the yminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -388,7 +389,7 @@ cpdef void update_pml_1order_magnetic_yminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHz
     cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0
@@ -424,7 +425,7 @@ cpdef void update_pml_1order_magnetic_yminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_1order_magnetic_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hz field components for the yplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -432,7 +433,7 @@ cpdef void update_pml_1order_magnetic_yplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHz
     cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0
@@ -468,7 +469,7 @@ cpdef void update_pml_1order_magnetic_yplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_1order_magnetic_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hy field components for the zminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -476,7 +477,7 @@ cpdef void update_pml_1order_magnetic_zminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHy
     cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0
@@ -512,7 +513,7 @@ cpdef void update_pml_1order_magnetic_zminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_1order_magnetic_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hy field components for the zplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -520,7 +521,7 @@ cpdef void update_pml_1order_magnetic_zplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHy
     cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0
diff --git a/gprMax/pml_2order_update.pyx b/gprMax/pml_2order_update.pyx
index 5cab5255..d6da0a75 100644
--- a/gprMax/pml_2order_update.pyx
+++ b/gprMax/pml_2order_update.pyx
@@ -20,7 +20,8 @@ import numpy as np
 cimport numpy as np
 from cython.parallel import prange
 
-from gprMax.constants cimport floattype_t, complextype_t
+from gprMax.constants cimport floattype_t
+from gprMax.constants cimport complextype_t
 
 
 ############################################
@@ -28,7 +29,7 @@ from gprMax.constants cimport floattype_t, complextype_t
 ############################################
 cpdef void update_pml_2order_electric_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ey and Ez field components for the xminus slab.
-        
+
         Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -36,7 +37,7 @@ cpdef void update_pml_2order_electric_xminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEy, materialEz
     cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -44,7 +45,7 @@ cpdef void update_pml_2order_electric_xminus(int xs, int xf, int ys, int yf, int
     nx = xf - xs
     ny = yf - ys
     nz = zf - zs
-    
+
     for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
         RA0 = ERA[0, i]
         RB0 = ERB[0, i]
@@ -79,7 +80,7 @@ cpdef void update_pml_2order_electric_xminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_2order_electric_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ey and Ez field components for the xplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -87,7 +88,7 @@ cpdef void update_pml_2order_electric_xplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEy, materialEz
     cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -130,7 +131,7 @@ cpdef void update_pml_2order_electric_xplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_2order_electric_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ez field components for the yminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -138,7 +139,7 @@ cpdef void update_pml_2order_electric_yminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEz
     cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -181,7 +182,7 @@ cpdef void update_pml_2order_electric_yminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_2order_electric_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ez field components for the yplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -189,7 +190,7 @@ cpdef void update_pml_2order_electric_yplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEz
     cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -232,7 +233,7 @@ cpdef void update_pml_2order_electric_yplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_2order_electric_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ey field components for the zminus slab.
-        
+
         Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -240,7 +241,7 @@ cpdef void update_pml_2order_electric_zminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
         """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEy
     cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -248,7 +249,7 @@ cpdef void update_pml_2order_electric_zminus(int xs, int xf, int ys, int yf, int
     nx = xf - xs
     ny = yf - ys
     nz = zf - zs
-    
+
     for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
         ii = i + xs
         for j in range(0, ny):
@@ -283,7 +284,7 @@ cpdef void update_pml_2order_electric_zminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_2order_electric_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d):
     """This function updates the Ex and Ey field components for the zplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -291,7 +292,7 @@ cpdef void update_pml_2order_electric_zplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialEx, materialEy
     cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -334,7 +335,7 @@ cpdef void update_pml_2order_electric_zplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_2order_magnetic_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hy and Hz field components for the xminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -342,7 +343,7 @@ cpdef void update_pml_2order_magnetic_xminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, ERE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHy, materialHz
     cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -385,7 +386,7 @@ cpdef void update_pml_2order_magnetic_xminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_2order_magnetic_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hy and Hz field components for the xplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -393,7 +394,7 @@ cpdef void update_pml_2order_magnetic_xplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHy, materialHz
     cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -436,7 +437,7 @@ cpdef void update_pml_2order_magnetic_xplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_2order_magnetic_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hz field components for the yminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -444,7 +445,7 @@ cpdef void update_pml_2order_magnetic_yminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHz
     cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -487,7 +488,7 @@ cpdef void update_pml_2order_magnetic_yminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_2order_magnetic_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hz field components for the yplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -495,7 +496,7 @@ cpdef void update_pml_2order_magnetic_yplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHz
     cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -538,7 +539,7 @@ cpdef void update_pml_2order_magnetic_yplus(int xs, int xf, int ys, int yf, int
 ############################################
 cpdef void update_pml_2order_magnetic_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hy field components for the zminus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -546,7 +547,7 @@ cpdef void update_pml_2order_magnetic_zminus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHy
     cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
@@ -589,7 +590,7 @@ cpdef void update_pml_2order_magnetic_zminus(int xs, int xf, int ys, int yf, int
 ###########################################
 cpdef void update_pml_2order_magnetic_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d):
     """This function updates the Hx and Hy field components for the zplus slab.
-        
+
     Args:
         xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box
         nthreads (int): Number of threads to use
@@ -597,7 +598,7 @@ cpdef void update_pml_2order_magnetic_zplus(int xs, int xf, int ys, int yf, int
         EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays
         d (float): Spatial discretisation, e.g. dx, dy or dz
     """
-    
+
     cdef Py_ssize_t i, j, k, ii, jj, kk
     cdef int nx, ny, nz, materialHx, materialHy
     cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01
diff --git a/gprMax/snapshots.py b/gprMax/snapshots.py
index 34431135..73b976b0 100644
--- a/gprMax/snapshots.py
+++ b/gprMax/snapshots.py
@@ -23,7 +23,9 @@ from struct import pack
 import numpy as np
 
 from gprMax.constants import floattype
-from gprMax.grid import Ix, Iy, Iz
+from gprMax.grid import Ix
+from gprMax.grid import Iy
+from gprMax.grid import Iz
 from gprMax.utilities import round_value
 
 
diff --git a/gprMax/sources.py b/gprMax/sources.py
index c4def072..b9709ed6 100644
--- a/gprMax/sources.py
+++ b/gprMax/sources.py
@@ -20,8 +20,11 @@ from copy import deepcopy
 
 import numpy as np
 
-from gprMax.constants import c, floattype
-from gprMax.grid import Ix, Iy, Iz
+from gprMax.constants import c
+from gprMax.constants import floattype
+from gprMax.grid import Ix
+from gprMax.grid import Iy
+from gprMax.grid import Iz
 from gprMax.utilities import round_value
 
 
diff --git a/gprMax/utilities.py b/gprMax/utilities.py
index 794574b9..28c85fae 100644
--- a/gprMax/utilities.py
+++ b/gprMax/utilities.py
@@ -26,7 +26,9 @@ from shutil import get_terminal_size
 import sys
 import textwrap
 
-from colorama import init, Fore, Style
+from colorama import init
+from colorama import Fore
+from colorama import Style
 init()
 import numpy as np
 
@@ -157,10 +159,10 @@ def get_host_info():
     Returns:
         hostinfo (dict): Manufacturer and model of machine; description of CPU type, speed, cores; RAM; name and version of operating system.
     """
-    
+
     # Default to 'unknown' if any of the detection fails
     manufacturer = model = cpuID = sockets = threadspercore = 'unknown'
-    
+
     # Windows
     if sys.platform == 'win32':
         # Manufacturer/model
diff --git a/gprMax/yee_cell_build.pyx b/gprMax/yee_cell_build.pyx
index 5f09fe29..27383463 100644
--- a/gprMax/yee_cell_build.pyx
+++ b/gprMax/yee_cell_build.pyx
@@ -20,12 +20,17 @@ import numpy as np
 cimport numpy as np
 
 from gprMax.materials import Material
-from gprMax.yee_cell_setget_rigid cimport get_rigid_Ex, get_rigid_Ey, get_rigid_Ez, get_rigid_Hx, get_rigid_Hy, get_rigid_Hz
+from gprMax.yee_cell_setget_rigid cimport get_rigid_Ex
+from gprMax.yee_cell_setget_rigid cimport get_rigid_Ey
+from gprMax.yee_cell_setget_rigid cimport get_rigid_Ez
+from gprMax.yee_cell_setget_rigid cimport get_rigid_Hx
+from gprMax.yee_cell_setget_rigid cimport get_rigid_Hy
+from gprMax.yee_cell_setget_rigid cimport get_rigid_Hz
 
 
 cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2, int numID3, int numID4, int componentID, G):
     """This function creates a new material by averaging the dielectric properties of the surrounding cells.
-        
+
     Args:
         i, j, k (int): Cell coordinates.
         numID1, numID2, numID3, numID4 (int): Numeric IDs for materials in surrounding cells.
@@ -35,7 +40,7 @@ cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2,
 
     # Make an ID composed of the names of the four materials that will be averaged
     requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID + '+' + G.materials[numID3].ID + '+' + G.materials[numID4].ID
-    
+
     # Check if this material already exists
     tmp = requiredID.split('+')
     material = [x for x in G.materials if
@@ -43,7 +48,7 @@ cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2,
              x.ID.count(tmp[1]) == requiredID.count(tmp[1]) and
              x.ID.count(tmp[2]) == requiredID.count(tmp[2]) and
              x.ID.count(tmp[3]) == requiredID.count(tmp[3])]
-             
+
     if material:
         G.ID[componentID, i, j, k] = material[0].numID
     else:
@@ -59,13 +64,13 @@ cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2,
 
         # Append the new material object to the materials list
         G.materials.append(m)
-        
+
         G.ID[componentID, i, j, k] = newNumID
 
 
 cpdef void create_magnetic_average(int i, int j, int k, int numID1, int numID2, int componentID, G):
     """This function creates a new material by averaging the dielectric properties of the surrounding cells.
-        
+
     Args:
         i, j, k (int): Cell coordinates.
         numID1, numID2 (int): Numeric IDs for materials in surrounding cells.
@@ -75,14 +80,14 @@ cpdef void create_magnetic_average(int i, int j, int k, int numID1, int numID2,
 
     # Make an ID composed of the names of the two materials that will be averaged
     requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID
-    
+
     # Check if this material already exists
     tmp = requiredID.split('+')
     material = [x for x in G.materials if
                 (x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
                  x.ID.count(tmp[1]) == requiredID.count(tmp[1])) or
                 (x.ID.count(tmp[0]) % 2 == 0 and x.ID.count(tmp[1]) % 2 == 0)]
-    
+
     if material:
         G.ID[componentID, i, j, k] = material[0].numID
     else:
@@ -95,21 +100,21 @@ cpdef void create_magnetic_average(int i, int j, int k, int numID1, int numID2,
         m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se), axis=0)
         m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr), axis=0)
         m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm), axis=0)
-        
+
         # Append the new material object to the materials list
         G.materials.append(m)
-        
+
         G.ID[componentID, i, j, k] = newNumID
 
 
 cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.uint32_t[:, :, :, ::1] ID, G):
     """This function builds the electric field components in the ID array.
-        
+
     Args:
         solid, rigid, ID (memoryviews): Access to solid, rigid and ID arrays
         G (class): Grid class instance - holds essential parameters describing the model.
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int numID1, numID2, numID3, numID4, componentID
 
@@ -118,7 +123,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     for i in range(0, G.nx):
         for j in range(1, G.ny):
             for k in range(1, G.nz):
-                
+
                 # If rigid is True do not average
                 if get_rigid_Ex(i, j, k, rigidE):
                     pass
@@ -127,7 +132,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                     numID2 = solid[i, j - 1, k]
                     numID3 = solid[i, j - 1, k - 1]
                     numID4 = solid[i, j, k - 1]
-                    
+
                     # If all values are the same no need to average
                     if numID1 == numID2 and numID1 == numID3 and numID1 == numID4:
                         ID[componentID, i, j, k] = numID1
@@ -140,7 +145,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     for i in range(1, G.nx):
         for j in range(0, G.ny):
             for k in range(1, G.nz):
-                
+
                 # If rigid is True do not average
                 if get_rigid_Ey(i, j, k, rigidE):
                     pass
@@ -149,7 +154,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                     numID2 = solid[i - 1, j, k]
                     numID3 = solid[i - 1, j, k - 1]
                     numID4 = solid[i, j, k - 1]
-                    
+
                     # If all values are the same no need to average
                     if numID1 == numID2 and numID1 == numID3 and numID1 == numID4:
                         ID[componentID, i, j, k] = numID1
@@ -162,7 +167,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     for i in range(1, G.nx):
         for j in range(1, G.ny):
             for k in range(0, G.nz):
-                
+
                 # If rigid is True do not average
                 if get_rigid_Ez(i, j, k, rigidE):
                     pass
@@ -171,7 +176,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                     numID2 = solid[i - 1, j, k]
                     numID3 = solid[i - 1, j - 1, k]
                     numID4 = solid[i, j - 1, k]
-                    
+
                     # If all values are the same no need to average
                     if numID1 == numID2 and numID1 == numID3 and numID1 == numID4:
                         ID[componentID, i, j, k] = numID1
@@ -182,12 +187,12 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
 
 cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID, G):
     """This function builds the magnetic field components in the ID array.
-        
+
     Args:
         solid, rigid, ID (memoryviews): Access to solid, rigid and ID arrays
         G (class): Grid class instance - holds essential parameters describing the model.
     """
-    
+
     cdef Py_ssize_t i, j, k
     cdef int numID1, numID2, componentID
 
@@ -196,14 +201,14 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     for i in range(1, G.nx):
         for j in range(0, G.ny):
             for k in range(0, G.nz):
-                
+
                 # If rigid is True do not average
                 if get_rigid_Hx(i, j, k, rigidH):
                     pass
                 else:
                     numID1 = solid[i, j, k]
                     numID2 = solid[i - 1, j, k]
-                    
+
                     # If all values are the same no need to average
                     if numID1 == numID2:
                         ID[componentID, i, j, k] = numID1
@@ -216,14 +221,14 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     for i in range(0, G.nx):
         for j in range(1, G.ny):
             for k in range(0, G.nz):
-                
+
                 # If rigid is True do not average
                 if get_rigid_Hy(i, j, k, rigidH):
                     pass
                 else:
                     numID1 = solid[i, j, k]
                     numID2 = solid[i, j - 1, k]
-                    
+
                     # If all values are the same no need to average
                     if numID1 == numID2:
                         ID[4, i, j, k] = numID1
@@ -236,14 +241,14 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     for i in range(0, G.nx):
         for j in range(0, G.ny):
             for k in range(1, G.nz):
-                
+
                 # If rigid is True do not average
                 if get_rigid_Hz(i, j, k, rigidH):
                     pass
                 else:
                     numID1 = solid[i, j, k]
                     numID2 = solid[i, j, k - 1]
-                    
+
                     # If all values are the same no need to average
                     if numID1 == numID2:
                         ID[5, i, j, k] = numID1