From 69e3fe55c54266089a847bc3284eab5ad8e63c9a Mon Sep 17 00:00:00 2001
From: Craig Warren <craig.warren@northumbria.ac.uk>
Date: Fri, 4 Nov 2022 16:17:20 +0000
Subject: [PATCH] Updates to docstrings and formatting.

---
 gprMax/cmds_geometry/add_grass.py             |  34 +-
 gprMax/cmds_geometry/add_surface_roughness.py |  26 +-
 gprMax/cmds_geometry/add_surface_water.py     |  16 +-
 gprMax/cmds_geometry/box.py                   |  29 +-
 gprMax/cmds_geometry/cmds_geometry.py         |  42 +-
 gprMax/cmds_geometry/cylinder.py              |  24 +-
 gprMax/cmds_geometry/cylindrical_sector.py    |  37 +-
 gprMax/cmds_geometry/edge.py                  |  12 +-
 gprMax/cmds_geometry/fractal_box.py           |  37 +-
 gprMax/cmds_geometry/fractal_box_builder.py   |  18 +-
 gprMax/cmds_geometry/geometry_objects_read.py |   5 +-
 gprMax/cmds_geometry/plate.py                 |  12 +-
 gprMax/cmds_geometry/sphere.py                |  29 +-
 gprMax/cmds_geometry/triangle.py              |  34 +-
 .../fields_updates_dispersive_template.jinja  | 144 ++--
 gprMax/cython/fields_updates_hsg.pyx          | 449 +++++------
 gprMax/cython/fields_updates_normal.pyx       | 136 ++--
 gprMax/cython/fractals_generate.pyx           |  61 +-
 gprMax/cython/geometry_outputs.pyx            |  33 +-
 gprMax/cython/geometry_primitives.pyx         | 712 ++++++++++--------
 gprMax/cython/snapshots.pyx                   |  44 +-
 gprMax/cython/yee_cell_build.pyx              | 103 ++-
 gprMax/cython/yee_cell_setget_rigid.pyx       | 127 +++-
 gprMax/subgrids/base.py                       |   3 +-
 gprMax/subgrids/multi.py                      |  46 +-
 gprMax/subgrids/precursor_nodes.py            |   5 +-
 gprMax/subgrids/subgrid_hsg.py                |   6 +-
 gprMax/subgrids/updates.py                    |  30 +-
 gprMax/subgrids/user_objects.py               |  58 +-
 29 files changed, 1301 insertions(+), 1011 deletions(-)

diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py
index a3b4b628..64d27034 100644
--- a/gprMax/cmds_geometry/add_grass.py
+++ b/gprMax/cmds_geometry/add_grass.py
@@ -29,20 +29,20 @@ logger = logging.getLogger(__name__)
 
 
 class AddGrass(UserObjectGeometry):
-    """Allows you to add grass with roots to a FractalBox class in the model.
+    """Adds grass with roots to a FractalBox class in the model.
 
     Attributes:
-        p1: a list of the lower left (x,y,z) coordinates of a surface on a 
+        p1: list of the lower left (x,y,z) coordinates of a surface on a 
             FractalBox class.
-        p2: a list of the upper right (x,y,z) coordinates of a surface on a 
+        p2: list of the upper right (x,y,z) coordinates of a surface on a 
             FractalBox class.
-        frac_dim: a float for the fractal dimension which, for an orthogonal 
+        frac_dim:float for the fractal dimension which, for an orthogonal 
                     parallelepiped, should take values between zero and three.
-        limits: a list to define lower and upper limits for a range over which 
+        limits: list to define lower and upper limits for a range over which 
                     the height of the blades of grass can vary.
-        n_blades: an int for the number of blades of grass that should be 
+        n_blades:int for the number of blades of grass that should be 
                     applied to the surface area.
-        fractal_box_id: a string identifier for the FractalBox class that the 
+        fractal_box_id: string identifier for the FractalBox class that the 
                         grass should be applied to.
     """
 
@@ -191,7 +191,8 @@ class AddGrass(UserObjectGeometry):
                              'enough for the number of grass blades/roots specified')
             raise ValueError
 
-        # Scale the distribution so that the summation is equal to one, i.e. a probability distribution
+        # Scale the distribution so that the summation is equal to one, 
+        # i.e. a probability distribution
         surface.fractalsurface = surface.fractalsurface / np.sum(surface.fractalsurface)
 
         # Set location of grass blades using probability distribution
@@ -202,16 +203,23 @@ class AddGrass(UserObjectGeometry):
         R = np.random.RandomState(surface.seed)
         A = R.random_sample(n_blades)
 
-        # Locate the random numbers in the bins created by the 1D vector of probability values, and convert the 1D index back into a x, y index for the original surface.
-        bladesindex = np.unravel_index(np.digitize(A, probability1D), (surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))
+        # Locate the random numbers in the bins created by the 1D vector of 
+        # probability values, and convert the 1D index back into a x, y index 
+        # for the original surface.
+        bladesindex = np.unravel_index(np.digitize(A, probability1D), 
+                                       (surface.fractalsurface.shape[0], 
+                                        surface.fractalsurface.shape[1]))
 
         # Set the fractal range to minimum and maximum heights of the grass blades
         surface.fractalrange = fractalrange
 
-        # Set the fractal surface using the pre-calculated spatial distribution and a random height
-        surface.fractalsurface = np.zeros((surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))
+        # Set the fractal surface using the pre-calculated spatial distribution 
+        # and a random height
+        surface.fractalsurface = np.zeros((surface.fractalsurface.shape[0], 
+                                           surface.fractalsurface.shape[1]))
         for i in range(len(bladesindex[0])):
-            surface.fractalsurface[bladesindex[0][i], bladesindex[1][i]] = R.randint(surface.fractalrange[0], surface.fractalrange[1], size=1)
+            surface.fractalsurface[bladesindex[0][i], bladesindex[1][i]] = R.randint(surface.fractalrange[0], 
+                                                                                     surface.fractalrange[1], size=1)
 
         # Create grass geometry parameters
         g = Grass(n_blades)
diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py
index 50f77024..9def72ba 100644
--- a/gprMax/cmds_geometry/add_surface_roughness.py
+++ b/gprMax/cmds_geometry/add_surface_roughness.py
@@ -28,20 +28,20 @@ logger = logging.getLogger(__name__)
 
 
 class AddSurfaceRoughness(UserObjectGeometry):
-    """Allows you to add grass with roots to a FractalBox class in the model.
+    """Adds surface roughness to a FractalBox class in the model.
 
     Attributes:
-        p1: a list of the lower left (x,y,z) coordinates of a surface on a 
+        p1: list of the lower left (x,y,z) coordinates of a surface on a 
             FractalBox class.
-        p2: a list of the upper right (x,y,z) coordinates of a surface on a 
+        p2: list of the upper right (x,y,z) coordinates of a surface on a 
             FractalBox class.
-        frac_dim: a float for the fractal dimension which, for an orthogonal 
+        frac_dim: float for the fractal dimension which, for an orthogonal 
                     parallelepiped, should take values between zero and three.
-        weighting: a list with weightings in the first and second direction of 
+        weighting: list with weightings in the first and second direction of 
                     the surface.
-        limits: a list to define lower and upper limits for a range over which 
+        limits: ist to define lower and upper limits for a range over which 
                     the surface roughness can vary.
-        fractal_box_id: a string identifier for the FractalBox class
+        fractal_box_id: string identifier for the FractalBox class
                         that the surface roughness should be applied to.
         seed: (optional) float parameter which controls the seeding of the random 
                 number generator used to create the fractals.
@@ -119,7 +119,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
                 logger.exception(self.__str__() + ' can only be used on the external ' + 
                                  'surfaces of a fractal box')
                 raise ValueError
-            fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
+            fractalrange = (round_value(limits[0] / grid.dx), 
+                            round_value(limits[1] / grid.dx))
             # xminus surface
             if xs == volume.xs:
                 if fractalrange[0] < 0 or fractalrange[1] > volume.xf:
@@ -147,7 +148,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
                 logger.exception(self.__str__() + ' can only be used on the external ' + 
                                  'surfaces of a fractal box')
                 raise ValueError
-            fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
+            fractalrange = (round_value(limits[0] / grid.dy), 
+                            round_value(limits[1] / grid.dy))
             # yminus surface
             if ys == volume.ys:
                 if fractalrange[0] < 0 or fractalrange[1] > volume.yf:
@@ -175,7 +177,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
                 logger.exception(self.__str__() + ' can only be used on the external ' +
                                  'surfaces of a fractal box')
                 raise ValueError
-            fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
+            fractalrange = (round_value(limits[0] / grid.dz), 
+                            round_value(limits[1] / grid.dz))
             # zminus surface
             if zs == volume.zs:
                 if fractalrange[0] < 0 or fractalrange[1] > volume.zf:
@@ -212,7 +215,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
         # List of existing surfaces IDs
         existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
         if surface.surfaceID in existingsurfaceIDs:
-            logger.exception(self.__str__() + f' has already been used on the {surface.surfaceID} surface')
+            logger.exception(self.__str__() + f' has already been used on the ' +
+                             f'{surface.surfaceID} surface')
             raise ValueError
 
         surface.generate_fractal_surface()
diff --git a/gprMax/cmds_geometry/add_surface_water.py b/gprMax/cmds_geometry/add_surface_water.py
index 6718a604..10bb3a7a 100644
--- a/gprMax/cmds_geometry/add_surface_water.py
+++ b/gprMax/cmds_geometry/add_surface_water.py
@@ -28,17 +28,17 @@ logger = logging.getLogger(__name__)
 
 
 class AddSurfaceWater(UserObjectGeometry):
-    """Allows you to add surface water to a FractalBox class in the model.
+    """Adds surface water to a FractalBox class in the model.
 
     Attributes:
-        p1: a list of the lower left (x,y,z) coordinates of a surface on a 
+        p1: list of the lower left (x,y,z) coordinates of a surface on a 
             FractalBox class.
-        p2: a list of the upper right (x,y,z) coordinates of a surface on a 
+        p2: list of the upper right (x,y,z) coordinates of a surface on a 
             FractalBox class.
-        depth: a float that defines the depth of the water, which should be 
+        depth: float that defines the depth of the water, which should be 
                 specified relative to the dimensions of the #fractal_box that 
                 the surface water is being applied.
-        fractal_box_id: a string identifier for the FractalBox class that the 
+        fractal_box_id: string identifier for the FractalBox class that the 
                         surface water should be applied to.
     """
 
@@ -87,7 +87,8 @@ class AddSurfaceWater(UserObjectGeometry):
         xf, yf, zf = p2
 
         if depth <= 0:
-            logger.exception(self.__str__() + ' requires a positive value for the depth of water')
+            logger.exception(self.__str__() + ' requires a positive value for ' +
+                             'the depth of water')
             raise ValueError
 
         # Check for valid orientations
@@ -148,7 +149,8 @@ class AddSurfaceWater(UserObjectGeometry):
 
         surface = next((x for x in volume.fractalsurfaces if x.surfaceID == requestedsurface), None)
         if not surface:
-            logger.exception(self.__str__() + f' specified surface {requestedsurface} does not have a rough surface applied')
+            logger.exception(self.__str__() + f' specified surface {requestedsurface} ' +
+                             'does not have a rough surface applied')
             raise ValueError
 
         surface.filldepth = filldepthcells
diff --git a/gprMax/cmds_geometry/box.py b/gprMax/cmds_geometry/box.py
index 5fa2be69..bf0b19f8 100644
--- a/gprMax/cmds_geometry/box.py
+++ b/gprMax/cmds_geometry/box.py
@@ -28,16 +28,16 @@ logger = logging.getLogger(__name__)
 
 
 class Box(UserObjectGeometry):
-    """Allows you to introduce an orthogonal parallelepiped with specific 
-        properties into the model.
+    """Introduces an orthogonal parallelepiped with specific properties into 
+        the model.
 
     Attributes:
-        p1: a list of the lower left (x,y,z) coordinates of the parallelepiped.
-        p2: a list of the upper right (x,y,z) coordinates of the parallelepiped.
-        material_id: a string for the material identifier that must correspond 
+        p1: list of the lower left (x,y,z) coordinates of the parallelepiped.
+        p2: list of the upper right (x,y,z) coordinates of the parallelepiped.
+        material_id: string for the material identifier that must correspond 
                         to material that has already been defined.
-        material_ids: a list of material identifiers in the x, y, z directions.
-        averaging: a string (y or n) used to switch on and off dielectric smoothing.
+        material_ids: list of material identifiers in the x, y, z directions.
+        averaging: string (y or n) used to switch on and off dielectric smoothing.
     """
 
     def __init__(self, **kwargs):
@@ -124,15 +124,20 @@ class Box(UserObjectGeometry):
                 m = Material(numID, requiredID)
                 m.type = 'dielectric-smoothed'
                 # Create dielectric-smoothed constituents for material
-                m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
-                m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
-                m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
-                m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                m.er = np.mean((materials[0].er, materials[1].er, 
+                                materials[2].er), axis=0)
+                m.se = np.mean((materials[0].se, materials[1].se, 
+                                materials[2].se), axis=0)
+                m.mr = np.mean((materials[0].mr, materials[1].mr, 
+                                materials[2].mr), axis=0)
+                m.sm = np.mean((materials[0].mr, materials[1].mr, 
+                                materials[2].mr), axis=0)
 
                 # Append the new material object to the materials list
                 grid.materials.append(m)
 
-        build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
+        build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz, 
+                  averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
 
         dielectricsmoothing = 'on' if averaging else 'off'
 
diff --git a/gprMax/cmds_geometry/cmds_geometry.py b/gprMax/cmds_geometry/cmds_geometry.py
index 803d4959..5330342c 100644
--- a/gprMax/cmds_geometry/cmds_geometry.py
+++ b/gprMax/cmds_geometry/cmds_geometry.py
@@ -45,15 +45,15 @@ class UserObjectGeometry:
         return f'{self.hash}: {s[:-1]}'
 
     def create(self, grid, uip):
-        """Create object and add it to the grid."""
+        """Creates object and adds it to the grid."""
         pass
 
     def rotate(self, axis, angle, origin=None):
-        """Rotate object - specialised for each object."""
+        """Rotates object - specialised for each object."""
         pass
 
     def grid_name(self, grid):
-        """Return subgrid name for use with logging info. Return an empty
+        """Returns subgrid name for use with logging info. Returns an empty
             string if the grid is the main grid.
         """
 
@@ -64,16 +64,16 @@ class UserObjectGeometry:
 
 
 def rotate_point(p, axis, angle, origin=(0, 0, 0)):
-    """Rotate a point.
+    """Rotates a point.
     
     Args:
-        p: an array of coordinates of point (x, y, z).
-        axis: a string which defines the axis about which to perform rotation (x, y, or z).
-        angle: an int specifying the angle of rotation (degrees).
-        origin: a tuple defining the point about which to perform rotation (x, y, z).
+        p: array of coordinates of point (x, y, z).
+        axis: string which defines the axis about which to perform rotation (x, y, or z).
+        angle: int specifying the angle of rotation (degrees).
+        origin: tuple defining the point about which to perform rotation (x, y, z).
 
     Returns:
-        p: an array of coordinates of rotated point (x, y, z)
+        p: array of coordinates of rotated point (x, y, z)
     """
 
     origin = np.array(origin)
@@ -97,13 +97,13 @@ def rotate_2point_object(pts, axis, angle, origin=None):
     """Rotate a geometry object that is defined by 2 points.
     
     Args:
-        pts: an array ofcoordinates of points of object to be rotated.
-        axis: a string which defines the axis about which to perform rotation (x, y, or z).
-        angle: an int specifying the angle of rotation (degrees).
-        origin: a tuple defining the point about which to perform rotation (x, y, z).
+        pts: array ofcoordinates of points of object to be rotated.
+        axis: string which defines the axis about which to perform rotation (x, y, or z).
+        angle: int specifying the angle of rotation (degrees).
+        origin: tuple defining the point about which to perform rotation (x, y, z).
 
     Returns:
-        new_pts: an array of coordinates of points of rotated object.
+        new_pts: array of coordinates of points of rotated object.
     """
     
     # Use origin at centre of object if not given
@@ -152,18 +152,18 @@ def rotate_2point_object(pts, axis, angle, origin=None):
 
 
 def rotate_polarisation(p, polarisation, axis, angle, G):
-    """Rotate a geometry object that is defined by a point and polarisation.
+    """Rotates a geometry object that is defined by a point and polarisation.
     
     Args:
-        p: an array of coordinates of point (x, y, z).
-        polarisation: a string defining the current polarisation (x, y, or z).
-        axis: a string which defines the axis about which to perform rotation (x, y, or z).
-        angle: an int specifying the angle of rotation (degrees).
+        p: array of coordinates of point (x, y, z).
+        polarisation: string defining the current polarisation (x, y, or z).
+        axis: string which defines the axis about which to perform rotation (x, y, or z).
+        angle: int specifying the angle of rotation (degrees).
         G: FDTDGrid class describing a grid in a model.
 
     Returns:
-        pts: an array of coordinates of points of rotated object.
-        new_polarisation: a string defining the new polarisation (x, y, or z).
+        pts: array of coordinates of points of rotated object.
+        new_polarisation: string defining the new polarisation (x, y, or z).
     """
 
     if polarisation.lower() == 'x':
diff --git a/gprMax/cmds_geometry/cylinder.py b/gprMax/cmds_geometry/cylinder.py
index e1541a46..65814855 100644
--- a/gprMax/cmds_geometry/cylinder.py
+++ b/gprMax/cmds_geometry/cylinder.py
@@ -28,18 +28,18 @@ logger = logging.getLogger(__name__)
 
 
 class Cylinder(UserObjectGeometry):
-    """Allows you to introduce a circular cylinder into the model.
+    """Introduces a circular cylinder into the model.
 
     Attributes:
-        p1: a list of the coordinates (x,y,z) of the centre of the first face 
+        p1: list of the coordinates (x,y,z) of the centre of the first face 
             of the cylinder.
-        p2: a list of the coordinates (x,y,z) of the centre of the second face 
+        p2: list of the coordinates (x,y,z) of the centre of the second face 
             of the cylinder.
         r: float of the radius of the cylinder.
-        material_id: a string for the material identifier that must correspond 
+        material_id: string for the material identifier that must correspond 
                         to material that has already been defined.
-        material_ids: a list of material identifiers in the x, y, z directions.
-        averaging: a string (y or n) used to switch on and off dielectric smoothing.
+        material_ids: list of material identifiers in the x, y, z directions.
+        averaging: string (y or n) used to switch on and off dielectric smoothing.
     """
 
     def __init__(self, **kwargs):
@@ -113,10 +113,14 @@ class Cylinder(UserObjectGeometry):
                 m = Material(numID, requiredID)
                 m.type = 'dielectric-smoothed'
                 # Create dielectric-smoothed constituents for material
-                m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
-                m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
-                m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
-                m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                m.er = np.mean((materials[0].er, materials[1].er, 
+                                materials[2].er), axis=0)
+                m.se = np.mean((materials[0].se, materials[1].se, 
+                                materials[2].se), axis=0)
+                m.mr = np.mean((materials[0].mr, materials[1].mr, 
+                                materials[2].mr), axis=0)
+                m.sm = np.mean((materials[0].mr, materials[1].mr, 
+                                materials[2].mr), axis=0)
 
                 # Append the new material object to the materials list
                 grid.materials.append(m)
diff --git a/gprMax/cmds_geometry/cylindrical_sector.py b/gprMax/cmds_geometry/cylindrical_sector.py
index 88f6dce2..d1d78957 100644
--- a/gprMax/cmds_geometry/cylindrical_sector.py
+++ b/gprMax/cmds_geometry/cylindrical_sector.py
@@ -28,28 +28,27 @@ logger = logging.getLogger(__name__)
 
 
 class CylindricalSector(UserObjectGeometry):
-    """Allows you to introduce a cylindrical sector (shaped like a slice of pie) 
-        into the model.
+    """Introduces a cylindrical sector (shaped like a slice of pie) into the model.
 
     Attributes:
-        normal: a string for the direction of the axis of the cylinder from which 
+        normal: string for the direction of the axis of the cylinder from which 
                 the sector is defined and can be x, y, or z.
-        ctr1: a float for the first coordinate of the centre of the cylindrical 
+        ctr1: float for the first coordinate of the centre of the cylindrical 
                 sector.
-        ctr2: a float for the second coordinate of the centre of the cylindrical 
+        ctr2: float for the second coordinate of the centre of the cylindrical 
                 sector.
-        extent1: a float for the first thickness from the centre of the 
+        extent1: float for the first thickness from the centre of the 
                     cylindrical sector.
-        extent2: a float for the second thickness from the centre of the 
+        extent2: float for the second thickness from the centre of the 
                     cylindrical sector.
-        r: a float for the radius of the cylindrical sector.
-        start: a float for the starting angle (in degrees) for the cylindrical 
+        r: float for the radius of the cylindrical sector.
+        start: float for the starting angle (in degrees) for the cylindrical 
                 sector.
-        end: a float for the angle (in degrees) swept by the cylindrical sector.
-        material_id: a string for the material identifier that must correspond 
+        end: float for the angle (in degrees) swept by the cylindrical sector.
+        material_id: string for the material identifier that must correspond 
                         to material that has already been defined.
-        material_ids: a list of material identifiers in the x, y, z directions.
-        averaging: a string (y or n) used to switch on and off dielectric smoothing.
+        material_ids: list of material identifiers in the x, y, z directions.
+        averaging: string (y or n) used to switch on and off dielectric smoothing.
     """
 
     def __init__(self, **kwargs):
@@ -138,10 +137,14 @@ class CylindricalSector(UserObjectGeometry):
                     m = Material(numID, requiredID)
                     m.type = 'dielectric-smoothed'
                     # Create dielectric-smoothed constituents for material
-                    m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
-                    m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
-                    m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
-                    m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                    m.er = np.mean((materials[0].er, materials[1].er, 
+                                    materials[2].er), axis=0)
+                    m.se = np.mean((materials[0].se, materials[1].se, 
+                                    materials[2].se), axis=0)
+                    m.mr = np.mean((materials[0].mr, materials[1].mr, 
+                                    materials[2].mr), axis=0)
+                    m.sm = np.mean((materials[0].mr, materials[1].mr, 
+                                    materials[2].mr), axis=0)
 
                     # Append the new material object to the materials list
                     grid.materials.append(m)
diff --git a/gprMax/cmds_geometry/edge.py b/gprMax/cmds_geometry/edge.py
index 5f8940da..f858782c 100644
--- a/gprMax/cmds_geometry/edge.py
+++ b/gprMax/cmds_geometry/edge.py
@@ -28,12 +28,12 @@ logger = logging.getLogger(__name__)
 
 
 class Edge(UserObjectGeometry):
-    """Allows you to introduce a wire with specific properties into the model.
+    """Introduces a wire with specific properties into the model.
 
     Attributes:
-        p1: a list of the coordinates (x,y,z) of the starting point of the edge.
-        p2: a list of the coordinates (x,y,z) of the ending point of the edge.
-        material_id: a string for the material identifier that must correspond 
+        p1: list of the coordinates (x,y,z) of the starting point of the edge.
+        p2: list of the coordinates (x,y,z) of the ending point of the edge.
+        material_id: string for the material identifier that must correspond 
                         to material that has already been defined.
     """
 
@@ -49,14 +49,14 @@ class Edge(UserObjectGeometry):
         self.dorotate = True
 
     def __dorotate(self):
-        """Perform rotation."""
+        """Performs rotation."""
         pts = np.array([self.kwargs['p1'], self.kwargs['p2']])
         rot_pts = rotate_2point_object(pts, self.axis, self.angle, self.origin)
         self.kwargs['p1'] = tuple(rot_pts[0, :])
         self.kwargs['p2'] = tuple(rot_pts[1, :])
         
     def create(self, grid, uip):
-        """Create edge and add it to the grid."""
+        """Creates edge and adds it to the grid."""
         try:
             p1 = self.kwargs['p1']
             p2 = self.kwargs['p2']
diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py
index f03725bc..4f1ae5ac 100644
--- a/gprMax/cmds_geometry/fractal_box.py
+++ b/gprMax/cmds_geometry/fractal_box.py
@@ -27,27 +27,27 @@ logger = logging.getLogger(__name__)
 
 
 class FractalBox(UserObjectGeometry):
-    """Allows you to introduce an orthogonal parallelepiped with fractal 
-        distributed properties which are related to a mixing model or 
-        normal material into the model.
+    """Introduces an orthogonal parallelepiped with fractal distributed 
+        properties which are related to a mixing model or normal material into 
+        the model.
 
     Attributes:
-        p1: a list of the lower left (x,y,z) coordinates of the parallelepiped.
-        p2: a list of the upper right (x,y,z) coordinates of the parallelepiped.
-        frac_dim: a float for the fractal dimension which, for an orthogonal 
+        p1: list of the lower left (x,y,z) coordinates of the parallelepiped.
+        p2: list of the upper right (x,y,z) coordinates of the parallelepiped.
+        frac_dim: float for the fractal dimension which, for an orthogonal 
                     parallelepiped, should take values between zero and three.
-        weighting: a list of the weightings in the x, y, z direction of the 
+        weighting: list of the weightings in the x, y, z direction of the 
                     parallelepiped.
-        n_materials: an int of the number of materials to use for the fractal 
+        n_materials: int of the number of materials to use for the fractal 
                         distribution (defined according to the associated 
                         mixing model). This should be set to one if using a 
                         normal material instead of a mixing model.
-        mixing_model_id: a string identifier for the associated mixing model or 
+        mixing_model_id: string identifier for the associated mixing model or 
                             material.
-        id: a string identifier for the fractal box itself.
-        seed: seed: (optional) float parameter which controls the seeding of the random 
-                number generator used to create the fractals.
-        averaging: a string (y or n) used to switch on and off dielectric smoothing.
+        id: string identifier for the fractal box itself.
+        seed: (optional) float parameter which controls the seeding of the 
+                random number generator used to create the fractals.
+        averaging: string (y or n) used to switch on and off dielectric smoothing.
     """
 
     def __init__(self, **kwargs):
@@ -62,7 +62,7 @@ class FractalBox(UserObjectGeometry):
         self.dorotate = True
 
     def __dorotate(self):
-        """Perform rotation."""
+        """Performs rotation."""
         pts = np.array([self.kwargs['p1'], self.kwargs['p2']])
         rot_pts = rotate_2point_object(pts, self.axis, self.angle, self.origin)
         self.kwargs['p1'] = tuple(rot_pts[0, :])
@@ -94,7 +94,8 @@ class FractalBox(UserObjectGeometry):
             # Go with user specified averaging
             averagefractalbox = self.kwargs['averaging']
         except KeyError:
-            # If they havent specfied - default is no dielectric smoothing for a fractal box
+            # If they havent specified - default is no dielectric smoothing for 
+            # a fractal box.
             averagefractalbox = False
 
         p3 = uip.round_to_grid_static_point(p1)
@@ -124,7 +125,8 @@ class FractalBox(UserObjectGeometry):
                              'number of bins')
             raise ValueError
 
-        # Find materials to use to build fractal volume, either from mixing models or normal materials
+        # Find materials to use to build fractal volume, either from mixing 
+        # models or normal materials.
         mixingmodel = next((x for x in grid.mixingmodels if x.ID == mixing_model_id), None)
         material = next((x for x in grid.materials if x.ID == mixing_model_id), None)
         nbins = n_materials
@@ -134,7 +136,8 @@ class FractalBox(UserObjectGeometry):
                 logger.exception(self.__str__() + ' must be used with more than ' +
                                  'one material from the mixing model.')
                 raise ValueError
-            # Create materials from mixing model as number of bins now known from fractal_box command
+            # Create materials from mixing model as number of bins now known 
+            # from fractal_box command.
             mixingmodel.calculate_debye_properties(nbins, grid)
         elif not material:
             logger.exception(self.__str__() + f' mixing model or material with ' +
diff --git a/gprMax/cmds_geometry/fractal_box_builder.py b/gprMax/cmds_geometry/fractal_box_builder.py
index a8f10b42..4114ebd1 100644
--- a/gprMax/cmds_geometry/fractal_box_builder.py
+++ b/gprMax/cmds_geometry/fractal_box_builder.py
@@ -124,7 +124,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                                 # Add y, z coordinates to existing location
                                                 yy = int(j - volume.ys + y)
                                                 zz = int(k - volume.zs + z)
-                                                # If these coordinates are outwith fractal volume stop building the blade, otherwise set the mask for grass
+                                                # If these coordinates are outwith fractal volume stop building the blade, 
+                                                # otherwise set the mask for grass.
                                                 if yy < 0 or yy >= volume.mask.shape[1] or zz < 0 or zz >= volume.mask.shape[2]:
                                                     break
                                                 else:
@@ -146,7 +147,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                                 # Add y, z coordinates to existing location
                                                 yy = int(j - volume.ys + y)
                                                 zz = int(k - volume.zs + z)
-                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
+                                                # If these coordinates are outwith the fractal volume stop building the root, 
+                                                # otherwise set the mask for grass.
                                                 if yy < 0 or yy >= volume.mask.shape[1] or zz < 0 or zz >= volume.mask.shape[2]:
                                                     break
                                                 else:
@@ -192,7 +194,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                                 # Add x, z coordinates to existing location
                                                 xx = int(i - volume.xs + x)
                                                 zz = int(k - volume.zs + z)
-                                                # If these coordinates are outwith fractal volume stop building the blade, otherwise set the mask for grass
+                                                # If these coordinates are outwith fractal volume stop building the blade, 
+                                                # otherwise set the mask for grass.
                                                 if xx < 0 or xx >= volume.mask.shape[0] or zz < 0 or zz >= volume.mask.shape[2]:
                                                     break
                                                 else:
@@ -214,7 +217,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                                 # Add x, z coordinates to existing location
                                                 xx = int(i - volume.xs + x)
                                                 zz = int(k - volume.zs + z)
-                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
+                                                # If these coordinates are outwith the fractal volume stop building the root, 
+                                                # otherwise set the mask for grass.
                                                 if xx < 0 or xx >= volume.mask.shape[0] or zz < 0 or zz >= volume.mask.shape[2]:
                                                     break
                                                 else:
@@ -260,7 +264,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                                 # Add x, y coordinates to existing location
                                                 xx = int(i - volume.xs + x)
                                                 yy = int(j - volume.ys + y)
-                                                # If these coordinates are outwith the fractal volume stop building the blade, otherwise set the mask for grass
+                                                # If these coordinates are outwith the fractal volume stop building the blade, 
+                                                # otherwise set the mask for grass.
                                                 if xx < 0 or xx >= volume.mask.shape[0] or yy < 0 or yy >= volume.mask.shape[1]:
                                                     break
                                                 else:
@@ -282,7 +287,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                                 # Add x, y coordinates to existing location
                                                 xx = int(i - volume.xs + x)
                                                 yy = int(j - volume.ys + y)
-                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
+                                                # If these coordinates are outwith the fractal volume stop building the root, 
+                                                # otherwise set the mask for grass.
                                                 if xx < 0 or xx >= volume.mask.shape[0] or yy < 0 or yy >= volume.mask.shape[1]:
                                                     break
                                                 else:
diff --git a/gprMax/cmds_geometry/geometry_objects_read.py b/gprMax/cmds_geometry/geometry_objects_read.py
index 9e6865d3..d1277aa7 100644
--- a/gprMax/cmds_geometry/geometry_objects_read.py
+++ b/gprMax/cmds_geometry/geometry_objects_read.py
@@ -40,7 +40,7 @@ class GeometryObjectsRead(UserObjectGeometry):
         pass
 
     def create(self, grid, uip):
-        """Create the object and add it to the grid."""
+        """Creates the object and adds it to the grid."""
         try:
             p1 = self.kwargs['p1']
             geofile = self.kwargs['geofile']
@@ -88,7 +88,8 @@ class GeometryObjectsRead(UserObjectGeometry):
                 else:
                     material.type = 'imported'
 
-        # See if geometry object file exists at specified path and if not try input file directory
+        # See if geometry object file exists at specified path and if not try 
+        # input file directory.
         geofile = Path(geofile)
         if not geofile.exists():
             geofile = Path(config.sim_config.input_file_path.parent, geofile)
diff --git a/gprMax/cmds_geometry/plate.py b/gprMax/cmds_geometry/plate.py
index 81e226ff..13555424 100644
--- a/gprMax/cmds_geometry/plate.py
+++ b/gprMax/cmds_geometry/plate.py
@@ -28,14 +28,14 @@ logger = logging.getLogger(__name__)
 
 
 class Plate(UserObjectGeometry):
-    """Allows you to introduce a plate with specific properties into the model.
+    """Introduces a plate with specific properties into the model.
 
     Attributes:
-        p1: a list of the lower left (x,y,z) coordinates of the plate.
-        p2: a list of the upper right (x,y,z) coordinates of the plate.
-        material_id: a string for the material identifier that must correspond 
+        p1: list of the lower left (x,y,z) coordinates of the plate.
+        p2: list of the upper right (x,y,z) coordinates of the plate.
+        material_id: string for the material identifier that must correspond 
                         to material that has already been defined.
-        material_ids: a list of material identifiers in the x, y, z directions.
+        material_ids: list of material identifiers in the x, y, z directions.
     """
 
     def __init__(self, **kwargs):
@@ -50,7 +50,7 @@ class Plate(UserObjectGeometry):
         self.dorotate = True
 
     def __dorotate(self):
-        """Perform rotation."""
+        """Performs rotation."""
         pts = np.array([self.kwargs['p1'], self.kwargs['p2']])
         rot_pts = rotate_2point_object(pts, self.axis, self.angle, self.origin)
         self.kwargs['p1'] = tuple(rot_pts[0, :])
diff --git a/gprMax/cmds_geometry/sphere.py b/gprMax/cmds_geometry/sphere.py
index e2303b8f..ccb64240 100644
--- a/gprMax/cmds_geometry/sphere.py
+++ b/gprMax/cmds_geometry/sphere.py
@@ -28,16 +28,15 @@ logger = logging.getLogger(__name__)
 
 
 class Sphere(UserObjectGeometry):
-    """Allows you to introduce a spherical object with specific parameters 
-        into the model.
+    """Introduces a spherical object with specific parameters into the model.
 
     Attributes:
-        p1: a list of the coordinates (x,y,z) of the centre of the sphere.
-        r: a float of radius of the sphere.
-        material_id: a string for the material identifier that must correspond 
+        p1: list of the coordinates (x,y,z) of the centre of the sphere.
+        r: float of radius of the sphere.
+        material_id: string for the material identifier that must correspond 
                         to material that has already been defined.
-        material_ids: a list of material identifiers in the x, y, z directions.
-        averaging: a string (y or n) used to switch on and off dielectric smoothing.
+        material_ids: list of material identifiers in the x, y, z directions.
+        averaging: string (y or n) used to switch on and off dielectric smoothing.
     """
 
     def __init__(self, **kwargs):
@@ -104,15 +103,21 @@ class Sphere(UserObjectGeometry):
                 m = Material(numID, requiredID)
                 m.type = 'dielectric-smoothed'
                 # Create dielectric-smoothed constituents for material
-                m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
-                m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
-                m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
-                m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                m.er = np.mean((materials[0].er, materials[1].er, 
+                                materials[2].er), axis=0)
+                m.se = np.mean((materials[0].se, materials[1].se, 
+                                materials[2].se), axis=0)
+                m.mr = np.mean((materials[0].mr, materials[1].mr, 
+                                materials[2].mr), axis=0)
+                m.sm = np.mean((materials[0].mr, materials[1].mr, 
+                                materials[2].mr), axis=0)
 
                 # Append the new material object to the materials list
                 grid.materials.append(m)
 
-        build_sphere(xc, yc, zc, r, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
+        build_sphere(xc, yc, zc, r, grid.dx, grid.dy, grid.dz, numID, 
+                     numIDx, numIDy, numIDz, averaging, grid.solid, 
+                     grid.rigidE, grid.rigidH, grid.ID)
 
         dielectricsmoothing = 'on' if averaging else 'off'
         logger.info(self.grid_name(grid) + f"Sphere with centre {p2[0]:g}m, " +
diff --git a/gprMax/cmds_geometry/triangle.py b/gprMax/cmds_geometry/triangle.py
index 141e6d7f..f72e80f3 100644
--- a/gprMax/cmds_geometry/triangle.py
+++ b/gprMax/cmds_geometry/triangle.py
@@ -28,19 +28,19 @@ logger = logging.getLogger(__name__)
 
 
 class Triangle(UserObjectGeometry):
-    """Allows you to introduce a triangular patch or a triangular prism with 
-        specific properties into the model.
+    """Introduces a triangular patch or a triangular prism with specific 
+        properties into the model.
 
     Attributes:
-        p1: a list of the coordinates (x,y,z) of the first apex of the triangle.
-        p2: a list of the coordinates (x,y,z) of the second apex of the triangle.
-        p3: a list of the coordinates (x,y,z) of the third apex of the triangle.
-        thickness: a float for the thickness of the triangular prism. If the 
+        p1: list of the coordinates (x,y,z) of the first apex of the triangle.
+        p2: list of the coordinates (x,y,z) of the second apex of the triangle.
+        p3: list of the coordinates (x,y,z) of the third apex of the triangle.
+        thickness: float for the thickness of the triangular prism. If the 
                     thickness is zero then a triangular patch is created.
-        material_id: a string for the material identifier that must correspond 
+        material_id: string for the material identifier that must correspond 
                         to material that has already been defined.
-        material_ids: a list of material identifiers in the x, y, z directions.
-        averaging: a string (y or n) used to switch on and off dielectric smoothing.
+        material_ids: list of material identifiers in the x, y, z directions.
+        averaging: string (y or n) used to switch on and off dielectric smoothing.
     """
 
     def __init__(self, **kwargs):
@@ -48,14 +48,14 @@ class Triangle(UserObjectGeometry):
         self.hash = '#triangle'
 
     def rotate(self, axis, angle, origin=None):
-        """Set parameters for rotation."""
+        """Sets parameters for rotation."""
         self.axis = axis
         self.angle = angle
         self.origin = origin
         self.dorotate = True
 
     def __dorotate(self):
-        """Perform rotation."""
+        """Performs rotation."""
         p1 = rotate_point(self.kwargs['p1'], self.axis, self.angle, self.origin)
         p2 = rotate_point(self.kwargs['p2'], self.axis, self.angle, self.origin)
         p3 = rotate_point(self.kwargs['p3'], self.axis, self.angle, self.origin)
@@ -154,10 +154,14 @@ class Triangle(UserObjectGeometry):
                     m = Material(numID, requiredID)
                     m.type = 'dielectric-smoothed'
                     # Create dielectric-smoothed constituents for material
-                    m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
-                    m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
-                    m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
-                    m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                    m.er = np.mean((materials[0].er, materials[1].er, 
+                                    materials[2].er), axis=0)
+                    m.se = np.mean((materials[0].se, materials[1].se, 
+                                    materials[2].se), axis=0)
+                    m.mr = np.mean((materials[0].mr, materials[1].mr, 
+                                    materials[2].mr), axis=0)
+                    m.sm = np.mean((materials[0].mr, materials[1].mr, 
+                                    materials[2].mr), axis=0)
 
                     # Append the new material object to the materials list
                     grid.materials.append(m)
diff --git a/gprMax/cython/fields_updates_dispersive_template.jinja b/gprMax/cython/fields_updates_dispersive_template.jinja
index 8c49d69c..48e506bd 100644
--- a/gprMax/cython/fields_updates_dispersive_template.jinja
+++ b/gprMax/cython/fields_updates_dispersive_template.jinja
@@ -39,26 +39,26 @@ cdef extern from "complex.h" nogil:
 
 {% for item in functions %}
 cpdef void {{ item.name_a }}(
-                    int nx,
-                    int ny,
-                    int nz,
-                    int nthreads,
-                    int maxpoles,
-                    {{ item.field_type }}[:, ::1] updatecoeffsE,
-                    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
-                    np.uint32_t[:, :, :, ::1] ID,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
-                    {{ item.field_type }}[:, :, ::1] Ex,
-                    {{ item.field_type }}[:, :, ::1] Ey,
-                    {{ item.field_type }}[:, :, ::1] Ez,
-                    {{ item.field_type }}[:, :, ::1] Hx,
-                    {{ item.field_type }}[:, :, ::1] Hy,
-                    {{ item.field_type }}[:, :, ::1] Hz
-            ):
-    """This function updates the electric field components when dispersive 
-        materials (with multiple poles) are present.
+    int nx,
+    int ny,
+    int nz,
+    int nthreads,
+    int maxpoles,
+    {{ item.field_type }}[:, ::1] updatecoeffsE,
+    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
+    np.uint32_t[:, :, :, ::1] ID,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
+    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
+    {{ item.field_type }}[:, :, ::1] Ex,
+    {{ item.field_type }}[:, :, ::1] Ey,
+    {{ item.field_type }}[:, :, ::1] Ez,
+    {{ item.field_type }}[:, :, ::1] Hx,
+    {{ item.field_type }}[:, :, ::1] Hy,
+    {{ item.field_type }}[:, :, ::1] Hz
+):
+    """Updates the electric field components when dispersive materials 
+        (with multiple poles) are present.
 
     Args:
         nx, ny, nz: int for grid size in cells.
@@ -158,22 +158,22 @@ cpdef void {{ item.name_a }}(
 
 {% for item in functions %}
 cpdef void {{ item.name_b }}(
-                    int nx,
-                    int ny,
-                    int nz,
-                    int nthreads,
-                    int maxpoles,
-                    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
-                    np.uint32_t[:, :, :, ::1] ID,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
-                    {{ item.field_type }}[:, :, ::1] Ex,
-                    {{ item.field_type }}[:, :, ::1] Ey,
-                    {{ item.field_type }}[:, :, ::1] Ez
-            ):
-    """This function updates a temporary dispersive material array when 
-        disperisive materials (with multiple poles) are present.
+    int nx,
+    int ny,
+    int nz,
+    int nthreads,
+    int maxpoles,
+    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
+    np.uint32_t[:, :, :, ::1] ID,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
+    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
+    {{ item.field_type }}[:, :, ::1] Ex,
+    {{ item.field_type }}[:, :, ::1] Ey,
+    {{ item.field_type }}[:, :, ::1] Ez
+):
+    """Updates a temporary dispersive material array when disperisive materials 
+        (with multiple poles) are present.
 
     Args:
         nx, ny, nz: int for grid size in cells.
@@ -228,26 +228,26 @@ cpdef void {{ item.name_b }}(
 
 {% for item in functions %}
 cpdef void {{ item.name_a_1 }}(
-                    int nx,
-                    int ny,
-                    int nz,
-                    int nthreads,
-                    int maxpoles,
-                    {{ item.field_type }}[:, ::1] updatecoeffsE,
-                    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
-                    np.uint32_t[:, :, :, ::1] ID,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
-                    {{ item.field_type }}[:, :, ::1] Ex,
-                    {{ item.field_type }}[:, :, ::1] Ey,
-                    {{ item.field_type }}[:, :, ::1] Ez,
-                    {{ item.field_type }}[:, :, ::1] Hx,
-                    {{ item.field_type }}[:, :, ::1] Hy,
-                    {{ item.field_type }}[:, :, ::1] Hz
-            ):
-    """This function updates the electric field components when dispersive 
-        materials (with 1 pole) are present.
+    int nx,
+    int ny,
+    int nz,
+    int nthreads,
+    int maxpoles,
+    {{ item.field_type }}[:, ::1] updatecoeffsE,
+    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
+    np.uint32_t[:, :, :, ::1] ID,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
+    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
+    {{ item.field_type }}[:, :, ::1] Ex,
+    {{ item.field_type }}[:, :, ::1] Ey,
+    {{ item.field_type }}[:, :, ::1] Ez,
+    {{ item.field_type }}[:, :, ::1] Hx,
+    {{ item.field_type }}[:, :, ::1] Hy,
+    {{ item.field_type }}[:, :, ::1] Hz
+):
+    """Updates the electric field components when dispersive materials 
+        (with 1 pole) are present.
 
     Args:
         nx, ny, nz: int for grid size in cells.
@@ -339,22 +339,22 @@ cpdef void {{ item.name_a_1 }}(
 
 {% for item in functions %}
 cpdef void {{ item.name_b_1 }}(
-                    int nx,
-                    int ny,
-                    int nz,
-                    int nthreads,
-                    int maxpoles,
-                    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
-                    np.uint32_t[:, :, :, ::1] ID,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
-                    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
-                    {{ item.field_type }}[:, :, ::1] Ex,
-                    {{ item.field_type }}[:, :, ::1] Ey,
-                    {{ item.field_type }}[:, :, ::1] Ez
-            ):
-    """This function updates a temporary dispersive material array when 
-        disperisive materials (with 1 pole) are present.
+    int nx,
+    int ny,
+    int nz,
+    int nthreads,
+    int maxpoles,
+    {{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
+    np.uint32_t[:, :, :, ::1] ID,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tx,
+    {{ item.dispersive_type }}[:, :, :, ::1] Ty,
+    {{ item.dispersive_type }}[:, :, :, ::1] Tz,
+    {{ item.field_type }}[:, :, ::1] Ex,
+    {{ item.field_type }}[:, :, ::1] Ey,
+    {{ item.field_type }}[:, :, ::1] Ez
+):
+    """Updates a temporary dispersive material array when disperisive materials 
+        (with 1 pole) are present.
 
     Args:
         nx, ny, nz: int for grid size in cells.
diff --git a/gprMax/cython/fields_updates_hsg.pyx b/gprMax/cython/fields_updates_hsg.pyx
index 2d3d982f..4f03d77d 100644
--- a/gprMax/cython/fields_updates_hsg.pyx
+++ b/gprMax/cython/fields_updates_hsg.pyx
@@ -22,262 +22,267 @@ from cython.parallel import prange
 
 
 cpdef void cython_update_electric_os(
-        np.float64_t[:, :] updatecoeffsE,
-        np.uint32_t[:, :, :, :] ID,
-        int face,
-        int l_l,
-        int l_u,
-        int m_l,
-        int m_u,
-        size_t n_l,
-        size_t n_u,
-        int nwn,
-        size_t lookup_id,
-        np.float64_t[:, :, :] field,
-        np.float64_t[:, :, :] inc_field,
-        size_t co,
-        int sign_n,
-        int sign_f,
-        int mid,
-        int r,
-        int s,
-        int nb,
-        int nthreads
+    np.float64_t[:, :] updatecoeffsE,
+    np.uint32_t[:, :, :, :] ID,
+    int face,
+    int l_l,
+    int l_u,
+    int m_l,
+    int m_u,
+    size_t n_l,
+    size_t n_u,
+    int nwn,
+    size_t lookup_id,
+    np.float64_t[:, :, :] field,
+    np.float64_t[:, :, :] inc_field,
+    size_t co,
+    int sign_n,
+    int sign_f,
+    int mid,
+    int r,
+    int s,
+    int nb,
+    int nthreads
 ):
-        """
-            Args:
+    """
+    Args:
+        subgrid: Subgrid class instance.
+        n: string for the normal to the face to update.
+        nwn: int for the number of working cell in the normal direction
+                to the face.
+        lookup_id: int id of the H component to update at each node.
+        field: nparray of main grid field to be updated.
+        inc_field: nparray of incident sub_grid field.
+        co: int coefficient used by gprMax update equations which
+                is specific to the field component being updated.
+        sign_n: int 1 or -1 sign of the incident field on the near face.
+        sign_f: int 1 or -1 sign of the incident field on the far face.
+        mid: boolean for the H node midway along the lower edge?
+        r: int for subgrid ratio.
+        s: int for separation of inner and outer surfaces.
+        nb: int for number of boundary cells.
+    """
+    # Comments here as as per left and right face
 
-            subgrid: (Subgrid)
-            n: (String) the normal to the face to update
-            nwn: (Int) number of working cell in the normal direction
-            to the face
-            lookup_id: (Int) id of the H component we wish to update at
-            each node
-            field: (Numpy array) main grid field to be updated
-            inc_field: (Numpy array) incident sub_grid field
-            co: (Int) Coefficient used by gprMax update equations which
-            is specific to the field component being updated.
-            sign_n: (Int) 1 or -1 sign of the incident field on the near face.
-            sign_f: (Int) 1 or -1 sign of the incident field on the far face.
-            mid: (Bool) is the H node midway along the lower edge?
-            r = self.ratio
-            s = self.is_os_sep
-            nb = self.n_boundary_cells
-        """
-        # Comments here as as per left and right face
+    cdef Py_ssize_t l, m, l_s, m_s, n_s_l, n_s_r, material_e_l, material_e_r, i0, j0, k0, i1, j1, k1, i2, j2, k2, i3, j3, k3
+    cdef int os
+    cdef double inc_n, inc_f
 
-        cdef Py_ssize_t l, m, l_s, m_s, n_s_l, n_s_r, material_e_l, material_e_r, i0, j0, k0, i1, j1, k1, i2, j2, k2, i3, j3, k3
-        cdef int os
-        cdef double inc_n, inc_f
+    # Surface normal index for the subgrid near face h nodes (left i index)
+    n_s_l = nb - s * r - r + r // 2
+    # Surface normal index for the subgrid far face h nodes (right i index)
+    n_s_r = nb + nwn + s * r + r // 2
+    # OS at the left face
+    os = nb - r * s
 
-        # surface normal index for the subgrid near face h nodes (left i index)
-        n_s_l = nb - s * r - r + r // 2
-        # surface normal index for the subgrid far face h nodes (right i index)
-        n_s_r = nb + nwn + s * r + r // 2
-        # OS at the left face
-        os = nb - r * s
+    # Iterate over a slice of the main grid using dummy indices
+    for l in prange(l_l, l_u, nogil=True, schedule='static', num_threads=nthreads):
 
-        # Iterate over a slice of the main grid using dummy indices
-        for l in prange(l_l, l_u, nogil=True, schedule='static', num_threads=nthreads):
+        # Calculate the subgrid j component of the H nodes
+        # i.e. Hz node of the left or right face
+        if mid == 1:
+            l_s = os + (l - l_l) * r + r // 2
+        # i.e. the Hy node of the left or right face
+        else:
+            l_s = os + (l - l_l) * r
 
-            # Calculate the subgrid j component of the H nodes
-            # i.e. Hz node of the left or right face
+        for m in range(m_l, m_u):
+
+            # Calculate the subgrid k component of the H nodes
             if mid == 1:
-                l_s = os + (l - l_l) * r + r // 2
-            # i.e. the Hy node of the left or right face
+                m_s = os + (m - m_l) * r
             else:
-                l_s = os + (l - l_l) * r
+                m_s = os + (m - m_l) * r + r // 2
 
-            for m in range(m_l, m_u):
+            # left and right
+            if face == 2:
+                # main grid index
+                i0, j0, k0 = n_l, l, m
+                # equivalent subgrid index
+                i1, j1, k1 = n_s_l, l_s, m_s
+                i2, j2, k2 = n_u, l, m
+                i3, j3, k3 = n_s_r, l_s, m_s
+            # front and back
+            if face == 3:
+                i0, j0, k0 = l, n_l, m
+                i1, j1, k1 = l_s, n_s_l, m_s
+                i2, j2, k2 = l, n_u, m
+                i3, j3, k3 = l_s, n_s_r, m_s
+            # top bottom
+            if face == 1:
+                i0, j0, k0 = l, m, n_l
+                i1, j1, k1 = l_s, m_s, n_s_l
+                i2, j2, k2 = l, m, n_u
+                i3, j3, k3 = l_s, m_s, n_s_r
+            # Update the left face
 
-                # Calculate the subgrid k component of the H nodes
-                if mid == 1:
-                    m_s = os + (m - m_l) * r
-                else:
-                    m_s = os + (m - m_l) * r + r // 2
+            # Get the material at main grid index
+            material_e_l = ID[lookup_id, i0, j0, k0]
+            # Get the associated indident field from the subgrid
+            inc_n = inc_field[i1, j1, k1] * sign_n
+            # Update the main grid E field with the corrected H field
+            field[i0, j0, k0] += updatecoeffsE[material_e_l, co] * inc_n
 
-                # left and right
-                if face == 2:
-                    # main grid index
-                    i0, j0, k0 = n_l, l, m
-                    # equivalent subgrid index
-                    i1, j1, k1 = n_s_l, l_s, m_s
-                    i2, j2, k2 = n_u, l, m
-                    i3, j3, k3 = n_s_r, l_s, m_s
-                # front and back
-                if face == 3:
-                    i0, j0, k0 = l, n_l, m
-                    i1, j1, k1 = l_s, n_s_l, m_s
-                    i2, j2, k2 = l, n_u, m
-                    i3, j3, k3 = l_s, n_s_r, m_s
-                # top bottom
-                if face == 1:
-                    i0, j0, k0 = l, m, n_l
-                    i1, j1, k1 = l_s, m_s, n_s_l
-                    i2, j2, k2 = l, m, n_u
-                    i3, j3, k3 = l_s, m_s, n_s_r
-                # Update the left face
+            # Update the right face
+            material_e_r = ID[lookup_id, i2, j2, k2]
+            inc_f = inc_field[i3, j3, k3] * sign_f
+            field[i2, j2, k2] += updatecoeffsE[material_e_r, co] * inc_f
 
-                # Get the material at main grid index
-                material_e_l = ID[lookup_id, i0, j0, k0]
-                # Get the associated indident field from the subgrid
-                inc_n = inc_field[i1, j1, k1] * sign_n
-                # Update the main grid E field with the corrected H field
-                field[i0, j0, k0] += updatecoeffsE[material_e_l, co] * inc_n
-
-                # Update the right face
-                material_e_r = ID[lookup_id, i2, j2, k2]
-                inc_f = inc_field[i3, j3, k3] * sign_f
-                field[i2, j2, k2] += updatecoeffsE[material_e_r, co] * inc_f
 
 cpdef void cython_update_magnetic_os(
-        np.float64_t[:, :] updatecoeffsH,
-        np.uint32_t[:, :, :, :] ID,
-        int face,
-        int l_l,
-        int l_u,
-        int m_l,
-        int m_u,
-        size_t n_l,
-        size_t n_u,
-        int nwn,
-        size_t lookup_id,
-        np.float64_t[:, :, :] field,
-        np.float64_t[:, :, :] inc_field,
-        size_t co,
-        int sign_n,
-        int sign_f,
-        int mid,
-        int r,
-        int s,
-        int nb,
-        int nthreads
+    np.float64_t[:, :] updatecoeffsH,
+    np.uint32_t[:, :, :, :] ID,
+    int face,
+    int l_l,
+    int l_u,
+    int m_l,
+    int m_u,
+    size_t n_l,
+    size_t n_u,
+    int nwn,
+    size_t lookup_id,
+    np.float64_t[:, :, :] field,
+    np.float64_t[:, :, :] inc_field,
+    size_t co,
+    int sign_n,
+    int sign_f,
+    int mid,
+    int r,
+    int s,
+    int nb,
+    int nthreads
 ):
-        """
-        int r ratio,
-        int s is_os_sep,
-        int nb n_boundary_cells
-        """
+    """
+    Args:
+        r: int for subgrid ratio.
+        s: int for separation of inner and outer surfaces.
+        nb: int for number of boundary cells.
+    """
 
-        cdef Py_ssize_t l, m, l_s, m_s, n_s_l, n_s_r, material_e_l, material_e_r, i0, j0, k0, i1, j1, k1, i2, j2, k2, i3, j3, k3
-        cdef int os
-        cdef double inc_n, inc_f
+    cdef Py_ssize_t l, m, l_s, m_s, n_s_l, n_s_r, material_e_l, material_e_r, i0, j0, k0, i1, j1, k1, i2, j2, k2, i3, j3, k3
+    cdef int os
+    cdef double inc_n, inc_f
 
-        # i index (normal to os) for the subgrid near face e node
-        n_s_l = nb - r * s
-        # Normal index for the subgrid far face e node
-        n_s_r = nb + nwn + s * r
+    # i index (normal to os) for the subgrid near face e node
+    n_s_l = nb - r * s
+    # Normal index for the subgrid far face e node
+    n_s_r = nb + nwn + s * r
 
-        # os inner index for the sub grid
-        os = nb - r * s
+    # os inner index for the sub grid
+    os = nb - r * s
 
-        for l in prange(l_l, l_u, nogil=True, schedule='static', num_threads=nthreads):
+    for l in prange(l_l, l_u, nogil=True, schedule='static', num_threads=nthreads):
 
-            # y coord of the Ex field component
+        # y coord of the Ex field component
+        if mid == 1:
+            l_s = os + (l - l_l) * r + r // 2
+        # y coord of the Ez field component
+        else:
+            l_s = os + (l - l_l) * r
+
+        for m in range(m_l, m_u):
+
+            # z coordinate of the Ex node in the subgrid
             if mid == 1:
-                l_s = os + (l - l_l) * r + r // 2
-            # y coord of the Ez field component
+                m_s = os + (m - m_l) * r
             else:
-                l_s = os + (l - l_l) * r
+                m_s = os + (m - m_l) * r + r // 2
 
-            for m in range(m_l, m_u):
+            # associate the given indices with their i, j, k values
 
-                # z coordinate of the Ex node in the subgrid
-                if mid == 1:
-                    m_s = os + (m - m_l) * r
-                else:
-                    m_s = os + (m - m_l) * r + r // 2
+            # left and right
+            if face == 2:
+                # main grid index
+                i0, j0, k0 = n_l, l, m
+                # equivalent subgrid index
+                i1, j1, k1 = n_s_l, l_s, m_s
+                i2, j2, k2 = n_u, l, m
+                i3, j3, k3 = n_s_r, l_s, m_s
+            # front and back
+            if face == 3:
+                i0, j0, k0 = l, n_l, m
+                i1, j1, k1 = l_s, n_s_l, m_s
+                i2, j2, k2 = l, n_u, m
+                i3, j3, k3 = l_s, n_s_r, m_s
+            # top bottom
+            if face == 1:
+                i0, j0, k0 = l, m, n_l
+                i1, j1, k1 = l_s, m_s, n_s_l
+                i2, j2, k2 = l, m, n_u
+                i3, j3, k3 = l_s, m_s, n_s_r
 
-                # associate the given indices with their i, j, k values
+            material_e_l = ID[lookup_id, i0, j0, k0]
+            inc_n = inc_field[i1, j1, k1] * sign_n
 
-                # left and right
-                if face == 2:
-                    # main grid index
-                    i0, j0, k0 = n_l, l, m
-                    # equivalent subgrid index
-                    i1, j1, k1 = n_s_l, l_s, m_s
-                    i2, j2, k2 = n_u, l, m
-                    i3, j3, k3 = n_s_r, l_s, m_s
-                # front and back
-                if face == 3:
-                    i0, j0, k0 = l, n_l, m
-                    i1, j1, k1 = l_s, n_s_l, m_s
-                    i2, j2, k2 = l, n_u, m
-                    i3, j3, k3 = l_s, n_s_r, m_s
-                # top bottom
-                if face == 1:
-                    i0, j0, k0 = l, m, n_l
-                    i1, j1, k1 = l_s, m_s, n_s_l
-                    i2, j2, k2 = l, m, n_u
-                    i3, j3, k3 = l_s, m_s, n_s_r
+            # make sure these are the correct grid
+            field[i0, j0, k0] += updatecoeffsH[material_e_l, co] * inc_n
 
-                material_e_l = ID[lookup_id, i0, j0, k0]
-                inc_n = inc_field[i1, j1, k1] * sign_n
+            # Far face
+            material_e_r = ID[lookup_id, i2, j2, k2]
+            inc_f = inc_field[i3, j3, k3] * sign_f
+            field[i2, j2, k2] += updatecoeffsH[material_e_r, co] * inc_f
 
-                # make sure these are the correct grid
-                field[i0, j0, k0] += updatecoeffsH[material_e_l, co] * inc_n
-
-                # Far face
-                material_e_r = ID[lookup_id, i2, j2, k2]
-                inc_f = inc_field[i3, j3, k3] * sign_f
-                field[i2, j2, k2] += updatecoeffsH[material_e_r, co] * inc_f
 
 cpdef void cython_update_is(
-        int nwx,
-        int nwy,
-        int nwz,
-        np.float64_t[:, :] updatecoeffsE,
-        np.uint32_t[:, :, :, :] ID,
-        int n,
-        int offset,
-        int nwl,
-        int nwm,
-        int nwn,
-        int face,
-        np.float64_t[:, :, :] field,
-        np.float64_t[:, :] inc_field_l,
-        np.float64_t[:, :] inc_field_u,
-        Py_ssize_t lookup_id,
-        int sign_l,
-        int sign_u,
-        Py_ssize_t co,
-        int nthreads
+    int nwx,
+    int nwy,
+    int nwz,
+    np.float64_t[:, :] updatecoeffsE,
+    np.uint32_t[:, :, :, :] ID,
+    int n,
+    int offset,
+    int nwl,
+    int nwm,
+    int nwn,
+    int face,
+    np.float64_t[:, :, :] field,
+    np.float64_t[:, :] inc_field_l,
+    np.float64_t[:, :] inc_field_u,
+    Py_ssize_t lookup_id,
+    int sign_l,
+    int sign_u,
+    Py_ssize_t co,
+    int nthreads
     ):
 
-        cdef Py_ssize_t l, m, i1, j1, k1, i2, j2, k2, field_material_l, field_material_u, inc_i, inc_j
-        cdef double inc_l, inc_u, f_l, f_u
-        # for inner faces H nodes are 1 cell before n boundary cells
-        cdef int n_o = n + offset
+    """
+    Args:
+    """
 
-        for l in prange(n, nwl + n, nogil=True, schedule='static', num_threads=nthreads):
-            for m in range(n, nwm + n):
+    cdef Py_ssize_t l, m, i1, j1, k1, i2, j2, k2, field_material_l, field_material_u, inc_i, inc_j
+    cdef double inc_l, inc_u, f_l, f_u
+    # For inner faces H nodes are 1 cell before n boundary cells
+    cdef int n_o = n + offset
 
-                # bottom and top
-                if face == 1:
-                    i1, j1, k1 = l, m, n_o
-                    i2, j2, k2 = l, m, n + nwz
-                # left and right
-                if face == 2:
-                    i1, j1, k1 = n_o, l, m
-                    i2, j2, k2 = n + nwx, l, m
-                # front and back
-                if face == 3:
-                    i1, j1, k1 = l, n_o, m
-                    i2, j2, k2 = l, n + nwy, m
+    for l in prange(n, nwl + n, nogil=True, schedule='static', num_threads=nthreads):
+        for m in range(n, nwm + n):
 
-                inc_i = l - n
-                inc_j = m - n
+            # bottom and top
+            if face == 1:
+                i1, j1, k1 = l, m, n_o
+                i2, j2, k2 = l, m, n + nwz
+            # left and right
+            if face == 2:
+                i1, j1, k1 = n_o, l, m
+                i2, j2, k2 = n + nwx, l, m
+            # front and back
+            if face == 3:
+                i1, j1, k1 = l, n_o, m
+                i2, j2, k2 = l, n + nwy, m
 
-                field_material_l = ID[lookup_id, i1, j1, k1]
-                inc_l = inc_field_l[inc_i, inc_j]
-                # Additional field at i, j, k
-                f_l = updatecoeffsE[field_material_l, co] * inc_l * sign_l
-                # Set the new value
-                field[i1, j1, k1] += f_l
+            inc_i = l - n
+            inc_j = m - n
 
-                field_material_u = ID[lookup_id, i2, j2, k2]
-                inc_u = inc_field_u[inc_i, inc_j]
-                # Additional field at i, j, k
-                f_u = updatecoeffsE[field_material_u, co] * inc_u * sign_u
-                # Set the new value
-                field[i2, j2, k2] += f_u
+            field_material_l = ID[lookup_id, i1, j1, k1]
+            inc_l = inc_field_l[inc_i, inc_j]
+            # Additional field at i, j, k
+            f_l = updatecoeffsE[field_material_l, co] * inc_l * sign_l
+            # Set the new value
+            field[i1, j1, k1] += f_l
+
+            field_material_u = ID[lookup_id, i2, j2, k2]
+            inc_u = inc_field_u[inc_i, inc_j]
+            # Additional field at i, j, k
+            f_u = updatecoeffsE[field_material_u, co] * inc_u * sign_u
+            # Set the new value
+            field[i2, j2, k2] += f_u
diff --git a/gprMax/cython/fields_updates_normal.pyx b/gprMax/cython/fields_updates_normal.pyx
index b6921417..681fa9d1 100644
--- a/gprMax/cython/fields_updates_normal.pyx
+++ b/gprMax/cython/fields_updates_normal.pyx
@@ -23,29 +23,27 @@ from cython.parallel import prange
 from gprMax.config cimport float_or_double
 
 
-###############################################
-# Electric field updates - standard materials #
-###############################################
 cpdef void update_electric(
-                    int nx,
-                    int ny,
-                    int nz,
-                    int nthreads,
-                    float_or_double[:, ::1] updatecoeffsE,
-                    np.uint32_t[:, :, :, ::1] ID,
-                    float_or_double[:, :, ::1] Ex,
-                    float_or_double[:, :, ::1] Ey,
-                    float_or_double[:, :, ::1] Ez,
-                    float_or_double[:, :, ::1] Hx,
-                    float_or_double[:, :, ::1] Hy,
-                    float_or_double[:, :, ::1] Hz
-            ):
-    """This function updates the electric field components.
+    int nx,
+    int ny,
+    int nz,
+    int nthreads,
+    float_or_double[:, ::1] updatecoeffsE,
+    np.uint32_t[:, :, :, ::1] ID,
+    float_or_double[:, :, ::1] Ex,
+    float_or_double[:, :, ::1] Ey,
+    float_or_double[:, :, ::1] Ez,
+    float_or_double[:, :, ::1] Hx,
+    float_or_double[:, :, ::1] Hy,
+    float_or_double[:, :, ::1] Hz
+):
+    """Updates the electric field components (for standard materials).
 
     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
+        nx, ny, nz: ints for grid size in cells.
+        nthreads: int for number of threads to use.
+        updatecoeffs, ID, E, H: memoryviews to access update coefficients, 
+                                ID and field component arrays
     """
 
     cdef Py_ssize_t i, j, k
@@ -57,7 +55,9 @@ cpdef void update_electric(
             for j in range(1, ny):
                 for k in range(1, nz):
                     materialEx = ID[0, i, j, k]
-                    Ex[i, j, k] = updatecoeffsE[materialEx, 0] * Ex[i, j, k] + updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1])
+                    Ex[i, j, k] = (updatecoeffsE[materialEx, 0] * Ex[i, j, k] + 
+                                   updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - 
+                                   updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]))
 
     # 2D - Ey component
     elif ny == 1:
@@ -65,7 +65,9 @@ cpdef void update_electric(
             for j in range(0, ny):
                 for k in range(1, nz):
                     materialEy = ID[1, i, j, k]
-                    Ey[i, j, k] = updatecoeffsE[materialEy, 0] * Ey[i, j, k] + updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k])
+                    Ey[i, j, k] = (updatecoeffsE[materialEy, 0] * Ey[i, j, k] + 
+                                   updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - 
+                                   updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]))
 
     # 2D - Ez component
     elif nz == 1:
@@ -73,7 +75,9 @@ cpdef void update_electric(
             for j in range(1, ny):
                 for k in range(0, nz):
                     materialEz = ID[2, i, j, k]
-                    Ez[i, j, k] = updatecoeffsE[materialEz, 0] * Ez[i, j, k] + updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k])
+                    Ez[i, j, k] = (updatecoeffsE[materialEz, 0] * Ez[i, j, k] + 
+                                   updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - 
+                                   updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]))
 
     # 3D
     else:
@@ -83,52 +87,62 @@ cpdef void update_electric(
                     materialEx = ID[0, i, j, k]
                     materialEy = ID[1, i, j, k]
                     materialEz = ID[2, i, j, k]
-                    Ex[i, j, k] = updatecoeffsE[materialEx, 0] * Ex[i, j, k] + updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1])
-                    Ey[i, j, k] = updatecoeffsE[materialEy, 0] * Ey[i, j, k] + updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k])
-                    Ez[i, j, k] = updatecoeffsE[materialEz, 0] * Ez[i, j, k] + updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k])
+                    Ex[i, j, k] = (updatecoeffsE[materialEx, 0] * Ex[i, j, k] + 
+                                   updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - 
+                                   updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]))
+                    Ey[i, j, k] = (updatecoeffsE[materialEy, 0] * Ey[i, j, k] + 
+                                   updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - 
+                                   updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]))
+                    Ez[i, j, k] = (updatecoeffsE[materialEz, 0] * Ez[i, j, k] + 
+                                   updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - 
+                                   updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]))
 
         # Ex components at i = 0
         for j in prange(1, ny, nogil=True, schedule='static', num_threads=nthreads):
             for k in range(1, nz):
                 materialEx = ID[0, 0, j, k]
-                Ex[0, j, k] = updatecoeffsE[materialEx, 0] * Ex[0, j, k] + updatecoeffsE[materialEx, 2] * (Hz[0, j, k] - Hz[0, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[0, j, k] - Hy[0, j, k - 1])
+                Ex[0, j, k] = (updatecoeffsE[materialEx, 0] * Ex[0, j, k] + 
+                               updatecoeffsE[materialEx, 2] * (Hz[0, j, k] - Hz[0, j - 1, k]) - 
+                               updatecoeffsE[materialEx, 3] * (Hy[0, j, k] - Hy[0, j, k - 1]))
 
         # Ey components at j = 0
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for k in range(1, nz):
                 materialEy = ID[1, i, 0, k]
-                Ey[i, 0, k] = updatecoeffsE[materialEy, 0] * Ey[i, 0, k] + updatecoeffsE[materialEy, 3] * (Hx[i, 0, k] - Hx[i, 0, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, 0, k] - Hz[i - 1, 0, k])
+                Ey[i, 0, k] = (updatecoeffsE[materialEy, 0] * Ey[i, 0, k] + 
+                               updatecoeffsE[materialEy, 3] * (Hx[i, 0, k] - Hx[i, 0, k - 1]) - 
+                               updatecoeffsE[materialEy, 1] * (Hz[i, 0, k] - Hz[i - 1, 0, k]))
 
         # Ez components at k = 0
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 materialEz = ID[2, i, j, 0]
-                Ez[i, j, 0] = updatecoeffsE[materialEz, 0] * Ez[i, j, 0] + updatecoeffsE[materialEz, 1] * (Hy[i, j, 0] - Hy[i - 1, j, 0]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, 0] - Hx[i, j - 1, 0])
+                Ez[i, j, 0] = (updatecoeffsE[materialEz, 0] * Ez[i, j, 0] + 
+                               updatecoeffsE[materialEz, 1] * (Hy[i, j, 0] - Hy[i - 1, j, 0]) - 
+                               updatecoeffsE[materialEz, 2] * (Hx[i, j, 0] - Hx[i, j - 1, 0]))
 
 
-##########################
-# Magnetic field updates #
-##########################
 cpdef void update_magnetic(
-                    int nx,
-                    int ny,
-                    int nz,
-                    int nthreads,
-                    float_or_double[:, ::1] updatecoeffsH,
-                    np.uint32_t[:, :, :, ::1] ID,
-                    float_or_double[:, :, ::1] Ex,
-                    float_or_double[:, :, ::1] Ey,
-                    float_or_double[:, :, ::1] Ez,
-                    float_or_double[:, :, ::1] Hx,
-                    float_or_double[:, :, ::1] Hy,
-                    float_or_double[:, :, ::1] Hz
-            ):
-    """This function updates the magnetic field components.
+    int nx,
+    int ny,
+    int nz,
+    int nthreads,
+    float_or_double[:, ::1] updatecoeffsH,
+    np.uint32_t[:, :, :, ::1] ID,
+    float_or_double[:, :, ::1] Ex,
+    float_or_double[:, :, ::1] Ey,
+    float_or_double[:, :, ::1] Ez,
+    float_or_double[:, :, ::1] Hx,
+    float_or_double[:, :, ::1] Hy,
+    float_or_double[:, :, ::1] Hz
+):
+    """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
+        nx, ny, nz: ints for grid size in cells.
+        nthreads: int for number of threads to use.
+        updatecoeffs, ID, E, H: memoryviews to access update coefficients, 
+                                ID and field component arrays
     """
 
     cdef Py_ssize_t i, j, k
@@ -142,7 +156,9 @@ cpdef void update_magnetic(
                 for j in range(0, ny):
                     for k in range(0, nz):
                         materialHx = ID[3, i, j, k]
-                        Hx[i, j, k] = updatecoeffsH[materialHx, 0] * Hx[i, j, k] - updatecoeffsH[materialHx, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) + updatecoeffsH[materialHx, 3] * (Ey[i, j, k + 1] - Ey[i, j, k])
+                        Hx[i, j, k] = (updatecoeffsH[materialHx, 0] * Hx[i, j, k] - 
+                                       updatecoeffsH[materialHx, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) + 
+                                       updatecoeffsH[materialHx, 3] * (Ey[i, j, k + 1] - Ey[i, j, k]))
 
         # Hy component
         if nx == 1 or nz == 1:
@@ -150,7 +166,9 @@ cpdef void update_magnetic(
                 for j in range(1, ny):
                     for k in range(0, nz):
                         materialHy = ID[4, i, j, k]
-                        Hy[i, j, k] = updatecoeffsH[materialHy, 0] * Hy[i, j, k] - updatecoeffsH[materialHy, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) + updatecoeffsH[materialHy, 1] * (Ez[i + 1, j, k] - Ez[i, j, k])
+                        Hy[i, j, k] = (updatecoeffsH[materialHy, 0] * Hy[i, j, k] - 
+                                       updatecoeffsH[materialHy, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) + 
+                                       updatecoeffsH[materialHy, 1] * (Ez[i + 1, j, k] - Ez[i, j, k]))
 
         # Hz component
         if nx == 1 or ny == 1:
@@ -158,7 +176,9 @@ cpdef void update_magnetic(
                 for j in range(0, ny):
                     for k in range(1, nz):
                         materialHz = ID[5, i, j, k]
-                        Hz[i, j, k] = updatecoeffsH[materialHz, 0] * Hz[i, j, k] - updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k] - Ey[i, j, k]) + updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k] - Ex[i, j, k])
+                        Hz[i, j, k] = (updatecoeffsH[materialHz, 0] * Hz[i, j, k] - 
+                                       updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k] - Ey[i, j, k]) + 
+                                       updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k] - Ex[i, j, k]))
     # 3D
     else:
         for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
@@ -167,6 +187,12 @@ cpdef void update_magnetic(
                     materialHx = ID[3, i + 1, j, k]
                     materialHy = ID[4, i, j + 1, k]
                     materialHz = ID[5, i, j, k + 1]
-                    Hx[i + 1, j, k] = updatecoeffsH[materialHx, 0] * Hx[i + 1, j, k] - updatecoeffsH[materialHx, 2] * (Ez[i + 1, j + 1, k] - Ez[i + 1, j, k]) + updatecoeffsH[materialHx, 3] * (Ey[i + 1, j, k + 1] - Ey[i + 1, j, k])
-                    Hy[i, j + 1, k] = updatecoeffsH[materialHy, 0] * Hy[i, j + 1, k] - updatecoeffsH[materialHy, 3] * (Ex[i, j + 1, k + 1] - Ex[i, j + 1, k]) + updatecoeffsH[materialHy, 1] * (Ez[i + 1, j + 1, k] - Ez[i, j + 1, k])
-                    Hz[i, j, k + 1] = updatecoeffsH[materialHz, 0] * Hz[i, j, k + 1] - updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k + 1] - Ey[i, j, k + 1]) + updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k + 1] - Ex[i, j, k + 1])
+                    Hx[i + 1, j, k] = (updatecoeffsH[materialHx, 0] * Hx[i + 1, j, k] - 
+                                       updatecoeffsH[materialHx, 2] * (Ez[i + 1, j + 1, k] - Ez[i + 1, j, k]) + 
+                                       updatecoeffsH[materialHx, 3] * (Ey[i + 1, j, k + 1] - Ey[i + 1, j, k]))
+                    Hy[i, j + 1, k] = (updatecoeffsH[materialHy, 0] * Hy[i, j + 1, k] - 
+                                       updatecoeffsH[materialHy, 3] * (Ex[i, j + 1, k + 1] - Ex[i, j + 1, k]) + 
+                                       updatecoeffsH[materialHy, 1] * (Ez[i + 1, j + 1, k] - Ez[i, j + 1, k]))
+                    Hz[i, j, k + 1] = (updatecoeffsH[materialHz, 0] * Hz[i, j, k + 1] - 
+                                       updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k + 1] - Ey[i, j, k + 1]) + 
+                                       updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k + 1] - Ex[i, j, k + 1]))
diff --git a/gprMax/cython/fractals_generate.pyx b/gprMax/cython/fractals_generate.pyx
index 0c59622a..a9fa610f 100644
--- a/gprMax/cython/fractals_generate.pyx
+++ b/gprMax/cython/fractals_generate.pyx
@@ -23,17 +23,29 @@ from cython.parallel import prange
 from gprMax.config cimport float_or_double_complex
 
 
-cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, ::1] A, float_or_double_complex[:, ::1] fractalsurface):
-    """This function generates a fractal surface for a 2D array.
+cpdef void generate_fractal2D(
+    int nx, 
+    int ny, 
+    int nthreads, 
+    int b, 
+    np.float64_t[:] weighting, 
+    np.float64_t[:] v1, 
+    np.complex128_t[:, ::1] A, 
+    float_or_double_complex[:, ::1] fractalsurface
+):
+    """Generates a fractal surface for a 2D array.
 
     Args:
-        nx, ny (int): Fractal surface size in cells
-        nthreads (int): Number of threads to use
-        b (int): Constant related to fractal dimension
-        weighting (memoryview): Access to weighting vector
-        v1 (memoryview): Access to positional vector at centre of array, scaled by weighting
-        A (memoryview): Access to array containing random numbers (to be convolved with fractal function)
-        fractalsurface (memoryview): Access to array containing fractal surface data
+        nx, ny: int for fractal surface size in cells.
+        nthreads: int for number of threads to use
+        b: int for constant related to fractal dimension.
+        weighting: memoryview for access to weighting vector.
+        v1: memoryview for access to positional vector at centre of array, 
+            scaled by weighting.
+        A: memoryview for access to array containing random numbers 
+            (to be convolved with fractal function).
+        fractalsurface: memoryview for access to array containing fractal 
+                        surface data.
     """
 
     cdef Py_ssize_t i, j
@@ -54,17 +66,30 @@ cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[
                 fractalsurface[i, j] = A[i, j] / B
 
 
-cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, :, ::1] A, float_or_double_complex[:, :, ::1] fractalvolume):
-    """This function generates a fractal volume for a 3D array.
+cpdef void generate_fractal3D(
+    int nx, 
+    int ny, 
+    int nz, 
+    int nthreads, 
+    int b, 
+    np.float64_t[:] weighting, 
+    np.float64_t[:] v1, 
+    np.complex128_t[:, :, ::1] A, 
+    float_or_double_complex[:, :, ::1] fractalvolume
+):
+    """Generates a fractal volume for a 3D array.
 
     Args:
-        nx, ny, nz (int): Fractal volume size in cells
-        nthreads (int): Number of threads to use
-        b (int): Constant related to fractal dimension
-        weighting (memoryview): Access to weighting vector
-        v1 (memoryview): Access to positional vector at centre of array, scaled by weighting
-        A (memoryview): Access to array containing random numbers (to be convolved with fractal function)
-        fractalvolume (memoryview): Access to array containing fractal volume data
+        nx, ny, nz: int for fractal volume size in cells.
+        nthreads: int for number of threads to use
+        b: int for constant related to fractal dimension.
+        weighting: memoryview for access to weighting vector.
+        v1: memoryview for access to positional vector at centre of array, 
+            scaled by weighting.
+        A: memoryview for access to array containing random numbers 
+            (to be convolved with fractal function).
+        fractalsurface: memoryview for access to array containing fractal 
+                        volume data.
     """
 
     cdef Py_ssize_t i, j, k
diff --git a/gprMax/cython/geometry_outputs.pyx b/gprMax/cython/geometry_outputs.pyx
index b1c9894a..5b45f227 100644
--- a/gprMax/cython/geometry_outputs.pyx
+++ b/gprMax/cython/geometry_outputs.pyx
@@ -20,21 +20,30 @@ import numpy as np
 cimport numpy as np
 
 
-cpdef write_lines(float xs, float ys, float zs, int nx, int ny, int nz, 
-                  float dx, float dy, float dz, np.uint32_t[:, :, :, :] ID):
-    """This function generates arrays with to be written as lines (cell edges) 
-        to a VTK file.
+cpdef write_lines(
+    float xs, 
+    float ys, 
+    float zs, 
+    int nx, 
+    int ny, 
+    int nz, 
+    float dx, 
+    float dy, 
+    float dz, 
+    np.uint32_t[:, :, :, :] ID
+):
+    """Generates arrays with to be written as lines (cell edges) to a VTK file.
 
     Args:
-        xs, ys, zs (float): Starting coordinates of geometry view in metres
-        nx, ny, nz (int): Size of the volume in cells
-        dx, dy, dz (float): Spatial discretisation of geometry view in metres
-        ID (nparray): Sampled ID array according to geometry view spatial 
-                            discretisation
+        xs, ys, zs: float for starting coordinates of geometry view in metres.
+        nx, ny, nz: int for size of the volume in cells.
+        dx, dy, dz: float for spatial discretisation of geometry view in metres.
+        ID: memoryview of sampled ID array according to geometry view spatial 
+            discretisation.
 
     Returns:
-        x, y, z (nparray): 1D arrays with coordinates of the vertex of the lines
-        lines (nparray): array of material IDs for each line (cell edge) required
+        x, y, z: 1D nparrays with coordinates of the vertex of the lines.
+        lines: nparray of material IDs for each line (cell edge) required.
     """
 
     cdef Py_ssize_t i, j, k
@@ -99,4 +108,4 @@ cpdef write_lines(float xs, float ys, float zs, int nx, int ny, int nz,
     y += ys
     z += zs
 
-    return x, y, z, lines 
+    return x, y, z, lines
\ No newline at end of file
diff --git a/gprMax/cython/geometry_primitives.pyx b/gprMax/cython/geometry_primitives.pyx
index 5fbb7aa8..66f3530a 100644
--- a/gprMax/cython/geometry_primitives.pyx
+++ b/gprMax/cython/geometry_primitives.pyx
@@ -34,15 +34,15 @@ from gprMax.utilities.utilities import round_value
 
 
 cpdef bint are_clockwise(
-                    float v1x,
-                    float v1y,
-                    float v2x,
-                    float v2y
-            ):
+    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.
+        v1x, v1y, v2x, v2y: floats of coordinates of vectors.
 
     Returns:
         (boolean)
@@ -52,15 +52,15 @@ cpdef bint are_clockwise(
 
 
 cpdef bint is_within_radius(
-                    float vx,
-                    float vy,
-                    float 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.
+        vx, vy: floats of coordinates of vector.
+        radius: float for radius.
 
     Returns:
         (boolean)
@@ -70,14 +70,14 @@ cpdef bint is_within_radius(
 
 
 cpdef bint is_inside_sector(
-                    float px,
-                    float py,
-                    float ctrx,
-                    float ctry,
-                    float sectorstartangle,
-                    float sectorangle,
-                    float radius
-            ):
+    float px,
+    float py,
+    float ctrx,
+    float ctry,
+    float sectorstartangle,
+    float sectorangle,
+    float radius
+):
     """For a point to be inside a circular sector, it has to meet the following tests:
         It has to be positioned anti-clockwise from the start "arm" of the sector
         It has to be positioned clockwise from the end arm of the sector
@@ -86,11 +86,11 @@ cpdef bint is_inside_sector(
         i.e. sector defined in an anti-clockwise direction
 
     Args:
-        px, py (float): Coordinates of point.
-        ctrx, ctry (float): Coordinates of centre of circle.
-        sectorstartangle (float): Angle (in radians) of start of sector.
-        sectorangle (float): Angle (in radians) that sector makes.
-        radius (float): Radius.
+        px, py: floats for coordinates of point.
+        ctrx, ctry: floats for coordinates of centre of circle.
+        sectorstartangle: float for angle (in radians) of start of sector.
+        sectorangle: float for angle (in radians) that sector makes.
+        radius: float for radius.
 
     Returns:
         (boolean)
@@ -111,18 +111,18 @@ cpdef bint is_inside_sector(
 
 
 cpdef bint point_in_polygon(
-                    float px,
-                    float py,
-                    list polycoords
-            ):
+    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.
+        px, py: float for coordinates of point to test.
+        polycoords: list of x, y tuples of coordinates that define the polygon.
 
     Returns:
-        inside (boolean)
+        inside: boolean
     """
 
     cdef int i
@@ -164,20 +164,20 @@ cpdef bint point_in_polygon(
 
 
 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
-            ):
+    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.
+        i, j, k: ints for cell coordinates of edge.
+        numIDx: int for numeric ID of material.
+        rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
     """
 
     set_rigid_Ex(i, j, k, rigidE)
@@ -185,20 +185,20 @@ cpdef void build_edge_x(
 
 
 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
-            ):
+    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.
+        i, j, k: ints for cell coordinates of edge.
+        numIDy: int for numeric ID of material.
+        rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
     """
 
     set_rigid_Ey(i, j, k, rigidE)
@@ -206,20 +206,20 @@ cpdef void build_edge_y(
 
 
 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
-            ):
+    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.
+        i, j, k: ints for cell coordinates of edge.
+        numIDz: int for numeric ID of material.
+        rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
     """
 
     set_rigid_Ez(i, j, k, rigidE)
@@ -227,21 +227,21 @@ cpdef void build_edge_z(
 
 
 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
-            ):
+    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.
+        i, j, k: ints for cell coordinates of the face.
+        numIDy, numIDz: int for numeric ID of material.
+        rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
     """
 
     set_rigid_Ey(i, j, k, rigidE)
@@ -255,21 +255,21 @@ cpdef void build_face_yz(
 
 
 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
-            ):
+    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.
+        i, j, k: ints for cell coordinates of the face.
+        numIDx, numIDz: int for numeric ID of material.
+        rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
     """
 
     set_rigid_Ex(i, j, k, rigidE)
@@ -283,21 +283,21 @@ cpdef void build_face_xz(
 
 
 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
-            ):
+    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.
+        i, j, k: ints for cell coordinates of the face.
+        numIDx, numIDy: int for numeric ID of material.
+        rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
     """
 
     set_rigid_Ex(i, j, k, rigidE)
@@ -311,26 +311,27 @@ cpdef void build_face_xy(
 
 
 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
-            ):
+    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.
+        i, j, k: ints for cell coordinates of voxel.
+        numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+        averaging: bint for whether material property averaging will occur for 
+                    the object.
+        solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
     """
 
     if averaging:
@@ -375,43 +376,44 @@ cpdef void build_voxel(
 
 
 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
-            ):
+    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.
+    Builds triangles and triangular prisms 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.
-        thickness (float): Thickness of the triangular prism.
-        dx, dy, dz (float): Spatial discretisation.
-        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.
+        x1, y1, z1, x2, y2, z2, x3, y3, z3: floats of coordinates of the vertices
+                                                of the triangular prism.
+        normal: string for normal direction to the plane of the triangular prism.
+        thickness: float for thickness of the triangular prism.
+        dx, dy, dz: floats for spatial discretisation.
+        numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+        averaging: bint for whether material property averaging will occur for 
+                    the object.
+        solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
     """
 
     cdef Py_ssize_t i, j, k
@@ -449,7 +451,8 @@ cpdef void build_triangle(
     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
+            # 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
                 jr = (j + 0.5) * dz
@@ -470,70 +473,78 @@ cpdef void build_triangle(
             if s > 0 and t > 0 and (s + t) < 2 * area * sign:
                 if thicknesscells == 0:
                     if normal == 'x':
-                        build_face_yz(level, i, j, numIDy, numIDz, rigidE, rigidH, ID)
+                        build_face_yz(level, i, j, numIDy, numIDz, 
+                                      rigidE, rigidH, ID)
                     elif normal == 'y':
-                        build_face_xz(i, level, j, numIDx, numIDz, rigidE, rigidH, ID)
+                        build_face_xz(i, level, j, numIDx, numIDz, 
+                                      rigidE, rigidH, ID)
                     elif normal == 'z':
-                        build_face_xy(i, j, level, numIDx, numIDy, rigidE, rigidH, ID)
+                        build_face_xy(i, j, level, numIDx, numIDy, 
+                                      rigidE, rigidH, ID)
                 else:
                     for k in range(level, level + thicknesscells):
                         if normal == 'x':
-                            build_voxel(k, i, j, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                            build_voxel(k, i, j, numID, numIDx, numIDy, numIDz, 
+                                        averaging, solid, rigidE, rigidH, ID)
                         elif normal == 'y':
-                            build_voxel(i, k, j, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                            build_voxel(i, k, j, numID, numIDx, numIDy, numIDz, 
+                                        averaging, solid, rigidE, rigidH, ID)
                         elif normal == 'z':
-                            build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                            build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                        averaging, solid, rigidE, rigidH, ID)
 
 
 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
-            ):
+    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.
+    Builds cylindrical sectors 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.
-        sectorstartangle (float): Angle (in radians) of start of sector.
-        sectorangle (float): Angle (in radians) that sector makes.
-        radius (float): Radius of the cylindrical sector.
-        normal (char): Normal direction to the plane of the cylindrical sector.
-        thickness (float): Thickness of the cylindrical sector.
-        dx, dy, dz (float): Spatial discretisation.
-        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.
+        ctr1, ctr2: floats for coordinates of centre of circle.
+        level: int for the third dimensional coordinate.
+        sectorstartangle: float for angle (in radians) of start of sector.
+        sectorangle: float for angle (in radians) that sector makes.
+        radius: float for radius of the cylindrical sector.
+        normal: string for normal direction to the plane of the triangular prism.
+        thickness: float for thickness of the triangular prism.
+        dx, dy, dz: floats for spatial discretisation.
+        numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+        averaging: bint for whether material property averaging will occur for 
+                    the object.
+        solid, rigidE, rigidH, ID: memoryviews to access 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
+        # Angles are defined from zero degrees on the positive y-axis going 
+        # towards positive z-axis.
         y1 = round_value((ctr1 - radius)/dy)
         y2 = round_value((ctr1 + radius)/dy)
         z1 = round_value((ctr2 - radius)/dz)
@@ -552,15 +563,19 @@ cpdef void build_cylindrical_sector(
 
         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):
+                if is_inside_sector(y * dy + 0.5 * dy, z * dz + 0.5 * dz, ctr1, 
+                                    ctr2, sectorstartangle, sectorangle, radius):
                     if thicknesscells == 0:
-                        build_face_yz(level, y, z, numIDy, numIDz, rigidE, rigidH, ID)
+                        build_face_yz(level, y, z, numIDy, numIDz, 
+                                      rigidE, rigidH, ID)
                     else:
                         for x in range(level, level + thicknesscells):
-                            build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                            build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, 
+                                        averaging, solid, rigidE, rigidH, ID)
 
     elif normal == 'y':
-        # Angles are defined from zero degrees on the positive x-axis going towards positive z-axis
+        # Angles are defined from zero degrees on the positive x-axis going 
+        # towards positive z-axis.
         x1 = round_value((ctr1 - radius)/dx)
         x2 = round_value((ctr1 + radius)/dx)
         z1 = round_value((ctr2 - radius)/dz)
@@ -579,15 +594,19 @@ cpdef void build_cylindrical_sector(
 
         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):
+                if is_inside_sector(x * dx + 0.5 * dx, z * dz + 0.5 * dz, ctr1, 
+                                    ctr2, sectorstartangle, sectorangle, radius):
                     if thicknesscells == 0:
-                        build_face_xz(x, level, z, numIDx, numIDz, rigidE, rigidH, ID)
+                        build_face_xz(x, level, z, numIDx, numIDz, 
+                                      rigidE, rigidH, ID)
                     else:
                         for y in range(level, level + thicknesscells):
-                            build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                            build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, 
+                                        averaging, solid, rigidE, rigidH, ID)
 
     elif normal == 'z':
-        # Angles are defined from zero degrees on the positive x-axis going towards positive y-axis
+        # Angles are defined from zero degrees on the positive x-axis going 
+        # towards positive y-axis.
         x1 = round_value((ctr1 - radius)/dx)
         x2 = round_value((ctr1 + radius)/dx)
         y1 = round_value((ctr2 - radius)/dy)
@@ -606,38 +625,42 @@ cpdef void build_cylindrical_sector(
 
         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):
+                if is_inside_sector(x * dx + 0.5 * dx, y * dy + 0.5 * dy, ctr1, 
+                                    ctr2, sectorstartangle, sectorangle, radius):
                     if thicknesscells == 0:
-                        build_face_xy(x, y, level, numIDx, numIDy, rigidE, rigidH, ID)
+                        build_face_xy(x, y, level, numIDx, numIDy, 
+                                      rigidE, rigidH, ID)
                     else:
                         for z in range(level, level + thicknesscells):
-                            build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                            build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, 
+                                        averaging, solid, rigidE, rigidH, ID)
 
 
 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.
+    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 boxes 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.
+        xs, xf, ys, yf, zs, zf: ints for cell coordinates of entire box.
+        numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+        averaging: bint for whether material property averaging will occur for 
+                    the object.
+        solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
     """
 
     cdef Py_ssize_t i, j, k
@@ -689,35 +712,38 @@ cpdef void build_box(
 
 
 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.
+    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 cylinders 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.
-        dx, dy, dz (float): Spatial discretisation.
-        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.
+        x1, y1, z1, x2, y2, z2: floats for coordinates of the centres of cylinder 
+                                faces.
+        r: float for radius of the cylinder.
+        dx, dy, dz: floats for spatial discretisation.
+        numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+        averaging: bint for whether material property averaging will occur for 
+                    the object.
+        solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
     """
 
     cdef Py_ssize_t i, j, k
@@ -729,15 +755,18 @@ cpdef void build_cylinder(
     # Check if cylinder is aligned with an axis
     x_align = y_align = z_align = 0
     # x-aligned
-    if round_value(y1 / dy) == round_value(y2 / dy) and round_value(z1 / dz) == round_value(z2 / dz):
+    if (round_value(y1 / dy) == round_value(y2 / dy) and 
+        round_value(z1 / dz) == round_value(z2 / dz)):
         x_align = 1
 
     # y-aligned
-    elif round_value(x1 / dx) == round_value(x2 / dx) and round_value(z1 / dz) == round_value(z2 / dz):
+    elif (round_value(x1 / dx) == round_value(x2 / dx) and 
+          round_value(z1 / dz) == round_value(z2 / dz)):
         y_align = 1
 
     # z-aligned
-    elif round_value(x1 / dx) == round_value(x2 / dx) and round_value(y1 / dy) == round_value(y2 / dy):
+    elif (round_value(x1 / dx) == round_value(x2 / dx) and 
+          round_value(y1 / dy) == round_value(y2 / dy):)
         z_align = 1
 
     # Calculate a bounding box for the cylinder
@@ -804,21 +833,24 @@ cpdef void build_cylinder(
             for k in range(zs, zf):
                 if np.sqrt((j * dy + 0.5 * dy - y1)**2 + (k * dz + 0.5 * dz - z1)**2) <= r:
                     for i in range(xs, xf):
-                        build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                        build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                    averaging, solid, rigidE, rigidH, ID)
     # y-aligned cylinder
     elif y_align:
         for i in range(xs, xf):
             for k in range(zs, zf):
                 if np.sqrt((i * dx + 0.5 * dx - x1)**2 + (k * dz + 0.5 * dz - z1)**2) <= r:
                     for j in range(ys, yf):
-                        build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                        build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                    averaging, solid, rigidE, rigidH, ID)
     # z-aligned cylinder
     elif z_align:
         for i in range(xs, xf):
             for j in range(ys, yf):
                 if np.sqrt((i * dx + 0.5 * dx - x1)**2 + (j * dy + 0.5 * dy - y1)**2) <= r:
                     for k in range(zs, zf):
-                        build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                        build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                    averaging, solid, rigidE, rigidH, ID)
 
     # Not aligned with any axis
     else:
@@ -836,9 +868,13 @@ cpdef void build_cylinder(
                     # Build flag - default false, set to True if point is in cylinder
                     build = 0
                     # Vector from centre of first cylinder face to test point
-                    f1pt = np.array([i * dx + 0.5 * dx - x1, j * dy + 0.5 * dy - y1, k * dz + 0.5 * dz - z1], dtype=np.float32)
+                    f1pt = np.array([i * dx + 0.5 * dx - x1, 
+                                     j * dy + 0.5 * dy - y1, 
+                                     k * dz + 0.5 * dz - z1], dtype=np.float32)
                     # Vector from centre of second cylinder face to test point
-                    f2pt = np.array([i * dx + 0.5 * dx - x2, j * dy + 0.5 * dy - y2, k * dz + 0.5 * dz - z2], dtype=np.float32)
+                    f2pt = np.array([i * dx + 0.5 * dx - x2, 
+                                     j * dy + 0.5 * dy - y2, 
+                                     k * dz + 0.5 * dz - z2], dtype=np.float32)
                     # Magnitudes
                     f1ptmag = np.sqrt((f1pt*f1pt).sum(axis=0))
                     f2ptmag = np.sqrt((f2pt*f2pt).sum(axis=0))
@@ -862,40 +898,44 @@ cpdef void build_cylinder(
                             theta2 = 0
                         distance1 = f1ptmag * np.sin(theta1)
                         distance2 = f2ptmag * np.sin(theta2)
-                        if (distance1 <= r or distance2 <= r) and theta1 <= np.pi/2 and theta2 <= np.pi/2:
+                        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)
+                        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.
+    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 spheres 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.
-        dx, dy, dz (float): Spatial discretisation.
-        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.
+        xc, yc, zc: ints for cell coordinates of the centre of the sphere.
+        r: float for radius of the sphere.
+        dx, dy, dz: floats for spatial discretisation.
+        numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+        averaging: bint for whether material property averaging will occur for 
+                    the object.
+        solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
     """
 
     cdef Py_ssize_t i, j, k
@@ -926,30 +966,36 @@ cpdef void build_sphere(
     for i in range(xs, xf):
         for j in range(ys, yf):
             for k in range(zs, zf):
-                if np.sqrt((i + 0.5 - xc)**2 * dx**2 + (j + 0.5 - yc)**2 * dy**2 + (k + 0.5 - zc)**2 * dz**2) <= r:
-                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                if (np.sqrt((i + 0.5 - xc)**2 * dx**2 + 
+                            (j + 0.5 - yc)**2 * dy**2 + 
+                            (k + 0.5 - zc)**2 * dz**2) <= r):
+                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                averaging, solid, rigidE, rigidH, ID)
 
 
 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
-            ):
+    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.
-        averaging (bint): Whether material property averaging will occur for the object.
-        data (memoryview): Access to array containing numeric IDs of voxels to create.
-        solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
+        xs, ys, zs: ints for cell coordinates of position of start of array in 
+                    domain.
+        numexistmaterials: int for number of existing materials in model prior 
+                            to building voxels.
+        averaging: bint for whether material property averaging will occur for 
+                    the object.
+        data: memoryview to access array containing numeric IDs of voxels to create.
+        solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
     """
 
     cdef Py_ssize_t i, j, k
@@ -987,28 +1033,29 @@ cpdef void build_voxels_from_array(
 
 
 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
-            ):
+    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.
-        averaging (bint): Whether material property averaging will occur for the object.
-        data (memoryview): Access to array containing numeric IDs of voxels to create.
-        mask (memoryview): Access to array containing a mask of voxels to create.
-        solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
+        xs, ys, zs: ints for cell coordinates of position of start of array in domain.
+        waternumID, grassnumID: ints for numeric ID of water and grass materials.
+        averaging: bint for whether material property averaging will occur for 
+                the object.
+        data: memoryview to access array containing numeric IDs of voxels to create.
+        mask: memoryview to access to array containing a mask of voxels to create.
+        solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
     """
 
     cdef Py_ssize_t i, j, k
@@ -1024,10 +1071,13 @@ cpdef void build_voxels_from_array_mask(
             for k in range(zs, zf):
                 if mask[i - xs, j - ys, k - zs] == 1:
                     numID = numIDx = numIDy = numIDz = data[i - xs, j - ys, k - zs]
-                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                averaging, solid, rigidE, rigidH, ID)
                 elif mask[i - xs, j - ys, k - zs] == 2:
                     numID = numIDx = numIDy = numIDz = waternumID
-                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                averaging, solid, rigidE, rigidH, ID)
                 elif mask[i - xs, j - ys, k - zs] == 3:
                     numID = numIDx = numIDy = numIDz = grassnumID
-                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
+                    build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, 
+                                averaging, solid, rigidE, rigidH, ID)
diff --git a/gprMax/cython/snapshots.pyx b/gprMax/cython/snapshots.pyx
index e621ada3..e14c41a8 100644
--- a/gprMax/cython/snapshots.pyx
+++ b/gprMax/cython/snapshots.pyx
@@ -20,31 +20,29 @@ from gprMax.config cimport float_or_double
 
 
 cpdef void calculate_snapshot_fields(
-                    int nx,
-                    int ny,
-                    int nz,
-                    float_or_double[:, :, ::1] sliceEx,
-                    float_or_double[:, :, ::1] sliceEy,
-                    float_or_double[:, :, ::1] sliceEz,
-                    float_or_double[:, :, ::1] sliceHx,
-                    float_or_double[:, :, ::1] sliceHy,
-                    float_or_double[:, :, ::1] sliceHz,
-                    float_or_double[:, :, ::1] snapEx,
-                    float_or_double[:, :, ::1] snapEy,
-                    float_or_double[:, :, ::1] snapEz,
-                    float_or_double[:, :, ::1] snapHx,
-                    float_or_double[:, :, ::1] snapHy,
-                    float_or_double[:, :, ::1] snapHz
-            ):
-    """This function calculates electric and magnetic values at points from
-        averaging values in cells.
+    int nx,
+    int ny,
+    int nz,
+    float_or_double[:, :, ::1] sliceEx,
+    float_or_double[:, :, ::1] sliceEy,
+    float_or_double[:, :, ::1] sliceEz,
+    float_or_double[:, :, ::1] sliceHx,
+    float_or_double[:, :, ::1] sliceHy,
+    float_or_double[:, :, ::1] sliceHz,
+    float_or_double[:, :, ::1] snapEx,
+    float_or_double[:, :, ::1] snapEy,
+    float_or_double[:, :, ::1] snapEz,
+    float_or_double[:, :, ::1] snapHx,
+    float_or_double[:, :, ::1] snapHy,
+    float_or_double[:, :, ::1] snapHz
+):
+    """Calculates electric and magnetic values at points from averaging values 
+        in cells.
 
     Args:
-        nx, ny, nz (int): Size of snapshot array
-        sliceEx, sliceEy, sliceEz,
-            sliceHx, sliceHy, sliceHz (memoryview): Access to slices of field arrays
-        snapEx, snapEy, snapEz,
-            snapHx, snapHy, snapHz (memoryview): Access to snapshot arrays
+        nx, ny, nz: ints for size of snapshot array.
+        slice: memoryviews to access slices of field arrays.
+        snap: memoryviews to access snapshot arrays.
     """
 
     cdef Py_ssize_t i, j, k
diff --git a/gprMax/cython/yee_cell_build.pyx b/gprMax/cython/yee_cell_build.pyx
index f597d822..194d45d5 100644
--- a/gprMax/cython/yee_cell_build.pyx
+++ b/gprMax/cython/yee_cell_build.pyx
@@ -28,18 +28,30 @@ from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Hz
 from gprMax.materials import Material
 
 
-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.
+cpdef void create_electric_average(
+    int i, 
+    int j, 
+    int k, 
+    int numID1, 
+    int numID2, 
+    int numID3, 
+    int numID4, 
+    int componentID, 
+    G
+):
+    """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.
-        componentID (int): Numeric ID for electric field component.
-        G (class): Grid class instance - holds essential parameters describing the model.
+        i, j, k: ints for cell coordinates.
+        numID: ints for numeric IDs for materials in surrounding cells.
+        componentID: int for numeric ID for electric field component.
+        G: FDTDGrid class describing a grid in a model.
     """
 
     # 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
+    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('+')
@@ -57,10 +69,14 @@ cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2,
         m = Material(newNumID, requiredID)
         m.type = 'dielectric-smoothed'
         # Create averaged constituents for material
-        m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er, G.materials[numID3].er, G.materials[numID4].er), axis=0)
-        m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se, G.materials[numID3].se, G.materials[numID4].se), axis=0)
-        m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr, G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
-        m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm, G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
+        m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er, 
+                        G.materials[numID3].er, G.materials[numID4].er), axis=0)
+        m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se, 
+                        G.materials[numID3].se, G.materials[numID4].se), axis=0)
+        m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr, 
+                        G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
+        m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm, 
+                        G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
 
         # Append the new material object to the materials list
         G.materials.append(m)
@@ -68,14 +84,23 @@ cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2,
         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.
+cpdef void create_magnetic_average(
+    int i, 
+    int j, 
+    int k, 
+    int numID1, 
+    int numID2, 
+    int componentID, 
+    G
+):
+    """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.
-        componentID (int): Numeric ID for electric field component.
-        G (class): Grid class instance - holds essential parameters describing the model.
+        i, j, k: ints for cell coordinates.
+        numID: ints for numeric IDs for materials in surrounding cells.
+        componentID: int for numeric ID for magnetic field component.
+        G: FDTDGrid class describing a grid in a model.
     """
 
     # Make an ID composed of the names of the two materials that will be averaged
@@ -107,12 +132,17 @@ cpdef void create_magnetic_average(int i, int j, int k, int numID1, int numID2,
         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.
+cpdef void build_electric_components(
+    np.uint32_t[:, :, ::1] solid, 
+    np.int8_t[:, :, :, ::1] rigidE, 
+    np.uint32_t[:, :, :, ::1] ID, 
+    G
+):
+    """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.
+        solid, rigid, ID: memoryviews to access solid, rigid and ID arrays.
+        G: FDTDGrid class describing a grid in a model.
     """
 
     cdef Py_ssize_t i, j, k
@@ -138,7 +168,8 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
+                        create_electric_average(i, j, k, numID1, numID2, 
+                                                numID3, numID4, componentID, G)
 
     # Ey component
     componentID = G.IDlookup['Ey']
@@ -160,7 +191,8 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
+                        create_electric_average(i, j, k, numID1, numID2, 
+                                                numID3, numID4, componentID, G)
 
     # Ez component
     componentID = G.IDlookup['Ez']
@@ -182,15 +214,21 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
+                        create_electric_average(i, j, k, numID1, numID2, 
+                                                numID3, numID4, componentID, G)
 
 
-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.
+cpdef void build_magnetic_components(
+    np.uint32_t[:, :, ::1] solid, 
+    np.int8_t[:, :, :, ::1] rigidH, 
+    np.uint32_t[:, :, :, ::1] ID, 
+    G
+):
+    """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.
+        solid, rigid, ID: memoryviews to access solid, rigid and ID arrays.
+        G: FDTDGrid class describing a grid in a model.
     """
 
     cdef Py_ssize_t i, j, k
@@ -214,7 +252,8 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
+                        create_magnetic_average(i, j, k, numID1, numID2, 
+                                                componentID, G)
 
     # Hy component
     componentID = G.IDlookup['Hy']
@@ -234,7 +273,8 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[4, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
+                        create_magnetic_average(i, j, k, numID1, numID2, 
+                                                componentID, G)
 
     # Hz component
     componentID = G.IDlookup['Hz']
@@ -254,4 +294,5 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[5, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
+                        create_magnetic_average(i, j, k, numID1, numID2, 
+                                                componentID, G)
diff --git a/gprMax/cython/yee_cell_setget_rigid.pyx b/gprMax/cython/yee_cell_setget_rigid.pyx
index 8eb0f21d..28fb888f 100644
--- a/gprMax/cython/yee_cell_setget_rigid.pyx
+++ b/gprMax/cython/yee_cell_setget_rigid.pyx
@@ -23,7 +23,13 @@ cimport numpy as np
 # Get and set functions for the rigid electric component array. The rigid array
 # is 4D with the 1st dimension holding the 12 electric edge components of a
 # cell - Ex1, Ex2, Ex3, Ex4, Ey1, Ey2, Ey3, Ey4, Ez1, Ez2, Ez3, Ez4
-cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+cdef bint get_rigid_Ex(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
+
     cdef bint result
     result = False
     if rigidE[0, i, j, k]:
@@ -39,7 +45,13 @@ cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
             result = True
     return result
 
-cdef bint get_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+
+cdef bint get_rigid_Ey(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
     cdef bint result
     result = False
     if rigidE[4, i, j, k]:
@@ -55,7 +67,13 @@ cdef bint get_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
             result = True
     return result
 
-cdef bint get_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+
+cdef bint get_rigid_Ez(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
     cdef bint result
     result = False
     if rigidE[8, i, j, k]:
@@ -71,7 +89,13 @@ cdef bint get_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
             result = True
     return result
 
-cdef void set_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+
+cdef void set_rigid_Ex(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
     rigidE[0, i, j, k] = True
     if j != 0:
         rigidE[1, i, j - 1, k] = True
@@ -80,7 +104,13 @@ cdef void set_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
     if j != 0 and k != 0:
         rigidE[2, i, j - 1, k - 1] = True
 
-cdef void set_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+
+cdef void set_rigid_Ey(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
     rigidE[4, i, j, k] = True
     if i != 0:
         rigidE[7, i - 1, j, k] = True
@@ -89,7 +119,13 @@ cdef void set_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
     if i != 0 and k != 0:
         rigidE[6, i - 1, j, k - 1] = True
 
-cdef void set_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+
+cdef void set_rigid_Ez(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
     rigidE[8, i, j, k] = True
     if i != 0:
         rigidE[9, i - 1, j, k] = True
@@ -98,16 +134,33 @@ cdef void set_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
     if i != 0 and j != 0:
         rigidE[10, i - 1, j - 1, k] = True
 
-cdef void set_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+
+cdef void set_rigid_E(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
     rigidE[:, i, j, k] = True
 
-cdef void unset_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
+
+cdef void unset_rigid_E(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidE
+):
     rigidE[:, i, j, k] = False
 
 # Get and set functions for the rigid magnetic component array. The rigid array
 # is 4D with the 1st dimension holding the 6 magnetic edge components - Hx1,
 # Hx2, Hy1, Hy2, Hz1, Hz2
-cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+cdef bint get_rigid_Hx(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     cdef bint result
     result = False
     if rigidH[0, i, j, k]:
@@ -117,7 +170,13 @@ cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
             result = True
     return result
 
-cdef bint get_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+
+cdef bint get_rigid_Hy(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     cdef bint result
     result = False
     if rigidH[2, i, j, k]:
@@ -127,7 +186,13 @@ cdef bint get_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
             result = True
     return result
 
-cdef bint get_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+
+cdef bint get_rigid_Hz(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     cdef bint result
     result = False
     if rigidH[4, i, j, k]:
@@ -137,23 +202,53 @@ cdef bint get_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
             result = True
     return result
 
-cdef void set_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+
+cdef void set_rigid_Hx(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     rigidH[0, i, j, k] = True
     if i != 0:
         rigidH[1, i - 1, j, k] = True
 
-cdef void set_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+
+cdef void set_rigid_Hy(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     rigidH[2, i, j, k] = True
     if j != 0:
         rigidH[3, i, j - 1, k] = True
 
-cdef void set_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+
+cdef void set_rigid_Hz(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     rigidH[4, i, j, k] = True
     if k != 0:
         rigidH[5, i, j, k - 1] = True
 
-cdef void set_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+
+cdef void set_rigid_H(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     rigidH[:, i, j, k] = True
 
-cdef void unset_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
+
+cdef void unset_rigid_H(
+    int i, 
+    int j, 
+    int k, 
+    np.int8_t[:, :, :, ::1] rigidH
+):
     rigidH[:, i, j, k] = False
diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py
index a163b2eb..f4c79fc5 100644
--- a/gprMax/subgrids/base.py
+++ b/gprMax/subgrids/base.py
@@ -17,7 +17,6 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import logging
-import numpy as np
 
 from ..grid import FDTDGrid
 
@@ -57,7 +56,7 @@ class SubGridBase(FDTDGrid):
 
         # Number of sub cells to extend the sub grid beyond the IS boundary
         d_to_pml = self.s_is_os_sep + self.pml_separation
-        # index of the is
+        # Index of the is
         self.n_boundary_cells = d_to_pml + self.pmlthickness['x0']
         self.n_boundary_cells_x = d_to_pml + self.pmlthickness['x0']
         self.n_boundary_cells_y = d_to_pml + self.pmlthickness['y0']
diff --git a/gprMax/subgrids/multi.py b/gprMax/subgrids/multi.py
index ad908041..598baf14 100644
--- a/gprMax/subgrids/multi.py
+++ b/gprMax/subgrids/multi.py
@@ -34,14 +34,14 @@ class ReferenceRx(Rx):
     logger.debug('ReferenceRx has no offset member.')
 
     def get_field(self, str_id, field):
-        """Return the field value at the equivalent coarse yee cell.
+        """Field value at the equivalent coarse yee cell.
 
         Args:
-            str_id (str): 'Ex' etc...
-            field (nparray): Numpy array of grid.Ez
+            str_id: string of 'Ex' etc...
+            field: nparray of grid.Ez
 
         Returns:
-            e (float): Field value.
+            e: float of field value.
         """
 
         d = {
@@ -58,13 +58,13 @@ class ReferenceRx(Rx):
         return e
 
     def get_Ex_from_field(self, Ex):
-        """Return the Ex field value from the equivalent coarse yee cell.
+        """Ex field value from the equivalent coarse yee cell.
 
         Args:
-            Ex (nparray): Numpy array of Ex field.
+            Ex: nparray of Ex field.
 
         Returns:
-            e (float): Ex field value.
+            e: float of Ex field value.
         """
 
         # offset = ratio // 2
@@ -72,62 +72,62 @@ class ReferenceRx(Rx):
         return e
 
     def get_Ey_from_field(self, Ey):
-        """Return the Ey field value from the equivalent coarse yee cell.
+        """Ey field value from the equivalent coarse yee cell.
 
         Args:
-            Ey (nparray): Numpy array of Ey field.
+            Ey: nparray of Ey field.
 
         Returns:
-            e (float): Ey field value.
+            e: float of Ey field value.
         """
         e = Ey[self.xcoord, self.ycoord + self.offset, self.zcoord]
         return e
 
     def get_Ez_from_field(self, Ez):
-        """Return the Ez field value from the equivalent coarse yee cell.
+        """Ez field value from the equivalent coarse yee cell.
 
         Args:
-            Ez (nparray): Numpy array of Ez field.
+            Ez: nparray of Ez field.
 
         Returns:
-            e (float): Ez field value.
+            e: float of Ez field value.
         """
 
         e = Ez[self.xcoord, self.ycoord, self.zcoord + self.offset]
         return e
 
     def get_Hx_from_field(self, Hx):
-        """Return the Hx field value from the equivalent coarse yee cell.
+        """Hx field value from the equivalent coarse yee cell.
 
         Args:
-            Hx (nparray): Numpy array of Hx field.
+            Hx: nparray of Hx field.
 
         Returns:
-            e (float): Hx field value.
+            e: float of Hx field value.
         """
         e = Hx[self.xcoord, self.ycoord + self.offset, self.zcoord + self.offset]
         return e
 
     def get_Hy_from_field(self, Hy):
-        """Return the Hy field value from the equivalent coarse yee cell.
+        """Hy field value from the equivalent coarse yee cell.
 
         Args:
-            Hy (nparray): Numpy array of Hy field.
+            Hy: nparray of Hy field.
 
         Returns:
-            e (float): Hy field value.
+            e: float of Hy field value.
         """
         e = Hy[self.xcoord + self.offset, self.ycoord, self.zcoord + self.offset]
         return e
 
     def get_Hz_from_field(self, Hz):
-        """Return the Hz field value from the equivalent coarse yee cell.
+        """Hz field value from the equivalent coarse yee cell.
 
         Args:
-            Hz (nparray): Numpy array of Hz field.
+            Hz: nparray of Hz field.
 
         Returns:
-            e (float): Hz field value.
+            e: float of Hz field value.
         """
         e = Hz[self.xcoord + self.offset, self.ycoord + self.offset, self.zcoord]
-        return e
+        return e
\ No newline at end of file
diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py
index 597fa338..844c3ac1 100644
--- a/gprMax/subgrids/precursor_nodes.py
+++ b/gprMax/subgrids/precursor_nodes.py
@@ -180,9 +180,8 @@ class PrecursorNodesBase:
             setattr(self, f, val)
 
     def calc_exact_field(self, field_names):
-        """Function to set the fields used in update calculations to the
-            values at the current main time step.
-            i.e. ey_left = copy.ey_left_1
+        """Sets the fields used in update calculations to the values at the 
+        current main time step, i.e. ey_left = copy.ey_left_1
         """
         for f in field_names:
             val = np.copy(getattr(self, f + '_1'))
diff --git a/gprMax/subgrids/subgrid_hsg.py b/gprMax/subgrids/subgrid_hsg.py
index 71c19137..2f6e34a7 100644
--- a/gprMax/subgrids/subgrid_hsg.py
+++ b/gprMax/subgrids/subgrid_hsg.py
@@ -36,7 +36,7 @@ class SubGridHSG(SubGridBase):
         self.gridtype = SubGridHSG.gridtype
 
     def update_magnetic_is(self, precursors):
-        """Update the subgrid nodes at the IS with the currents derived
+        """Updates the subgrid nodes at the IS with the currents derived
             from the main grid.
 
         Args:
@@ -59,7 +59,7 @@ class SubGridHSG(SubGridBase):
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwz, self.nwy, 3, self.Hx, precursors.ez_front, precursors.ez_back, self.IDlookup['Hx'], 1, -1, 2, config.get_model_config().ompthreads)
 
     def update_electric_is(self, precursors):
-        """Update the subgrid nodes at the IS with the currents derived
+        """Updates the subgrid nodes at the IS with the currents derived
             from the main grid.
 
         Args:
@@ -142,7 +142,7 @@ class SubGridHSG(SubGridBase):
         cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
 
     def print_info(self):
-        """Print information about the subgrid."""
+        """Prints information about the subgrid."""
         
         logger.info('')
         logger.info(f'[{self.name}] Type: {self.gridtype}')
diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py
index 5c80fb8c..f5a991cc 100644
--- a/gprMax/subgrids/updates.py
+++ b/gprMax/subgrids/updates.py
@@ -47,37 +47,37 @@ def create_updates(G):
 
 
 class SubgridUpdates(CPUUpdates):
-    """Provides update functions for the Sub gridding simulation."""
+    """Update functions for the Sub gridding simulation."""
 
     def __init__(self, G, updaters):
         super().__init__(G)
         self.updaters = updaters
 
     def hsg_1(self):
-        """Method to update the subgrids over the first phase."""
-        # updaters update each subgrid
+        """Updates the subgrids over the first phase."""
         for sg_updater in self.updaters:
             sg_updater.hsg_1()
 
     def hsg_2(self):
-        """Method to update the subgrids over the second phase."""
+        """Updates the subgrids over the second phase."""
         for sg_updater in self.updaters:
             sg_updater.hsg_2()
 
 
 class SubgridUpdater(CPUUpdates):
-    """Class to handle updating the electric and magnetic fields of an HSG
-        subgrid. The IS, OS, subgrid region and the electric/magnetic sources
-        are updated using the precursor regions.
+    """Handles updating the electric and magnetic fields of an HSG subgrid. 
+        The IS, OS, subgrid region and the electric/magnetic sources are updated 
+        using the precursor regions.
     """
 
     def __init__(self, subgrid, precursors, G):
         """
         Args:
-            subgrid (SubGrid3d): Subgrid to be updated.
-            precursors (PrecursorNodes): Precursor nodes associated with
-            the subgrid - contain interpolated fields.
-            G (FDTDGrid): Parameters describing a grid in a model.
+            subgrid: SubGrid3d instance to be updated.
+            precursors (PrecursorNodes): PrecursorNodes instance nodes associated 
+                                            with the subgrid - contain interpolated 
+                                            fields.
+            G: FDTDGrid class describing a grid in a model.
         """
         super().__init__(subgrid)
         self.precursors = precursors
@@ -85,8 +85,8 @@ class SubgridUpdater(CPUUpdates):
         self.source_iteration = 0
 
     def hsg_1(self):
-        """This is the first half of the subgrid update. Takes the time step
-            up to the main grid magnetic update.
+        """First half of the subgrid update. Takes the time step up to the main 
+            grid magnetic update.
         """
 
         G = self.G
@@ -126,8 +126,8 @@ class SubgridUpdater(CPUUpdates):
         sub_grid.update_electric_os(G)
 
     def hsg_2(self):
-        """This is the first half of the subgrid update. Takes the time step
-            up to the main grid electric update.
+        """Second half of the subgrid update. Takes the time step up to the main 
+            grid electric update.
         """
 
         G = self.G
diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py
index e91c85b7..65f13c92 100644
--- a/gprMax/subgrids/user_objects.py
+++ b/gprMax/subgrids/user_objects.py
@@ -30,18 +30,17 @@ logger = logging.getLogger(__name__)
 
 
 class SubGridBase(UserObjectMulti):
-    """Class to allow UserObjectMulti and UserObjectGeometry to be nested
-        in SubGrid type user objects.
+    """Allows UserObjectMulti and UserObjectGeometry to be nested in SubGrid 
+        type user objects.
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.children_multiple = []
         self.children_geometry = []
 
     def add(self, node):
-        """Function to add other user objects. Geometry and multi only."""
+        """Adds other user objects. Geometry and multi only."""
         if isinstance(node, UserObjectMulti):
             self.children_multiple.append(node)
         elif isinstance(node, UserObjectGeometry):
@@ -51,14 +50,13 @@ class SubGridBase(UserObjectMulti):
             raise ValueError
 
     def set_discretisation(self, sg, grid):
-        """Set the spatial discretisation."""
         sg.dx = grid.dx / sg.ratio
         sg.dy = grid.dy / sg.ratio
         sg.dz = grid.dz / sg.ratio
         sg.dl = np.array([sg.dx, sg.dy, sg.dz])
 
     def set_main_grid_indices(self, sg, grid, uip, p1, p2):
-        """Set subgrid indices related to main grid placement."""
+        """Sets subgrid indices related to main grid placement."""
         # location of the IS
         sg.i0, sg.j0, sg.k0 = p1
         sg.i1, sg.j1, sg.k1 = p2
@@ -82,7 +80,7 @@ class SubGridBase(UserObjectMulti):
         sg.nz = 2 * sg.n_boundary_cells_z + sg.nwz
 
     def set_iterations(self, sg, main):
-        """Set number of iterations that will take place in the subgrid."""
+        """Sets number of iterations that will take place in the subgrid."""
         sg.iterations = main.iterations * sg.ratio
 
     def setup(self, sg, grid, uip):
@@ -97,7 +95,7 @@ class SubGridBase(UserObjectMulti):
         # Set the temporal discretisation
         sg.calculate_dt()
 
-        # set the indices related to the subgrids main grid placement
+        # Set the indices related to the subgrids main grid placement
         self.set_main_grid_indices(sg, grid, uip, p1, p2)
 
         """
@@ -121,8 +119,8 @@ class SubGridBase(UserObjectMulti):
 
         sg.timewindow = grid.timewindow
 
-        # Copy a subgrid reference to self so that children.create(grid, uip) can access
-        # the correct grid
+        # Copy a subgrid reference to self so that children.create(grid, uip) 
+        # can access the correct grid.
         self.subgrid = sg
 
         # Copy over built in materials
@@ -141,25 +139,25 @@ class SubGridBase(UserObjectMulti):
 class SubGridHSG(SubGridBase):
     """Huygens Surface subgridding (HSG) user object.
 
-    :param p1: Position of the lower left corner of the Inner Surface (x, y, z) in the main grid.
-    :type p1: list, non-optional
-    :param p2: Position of the upper right corner of the Inner Surface (x, y, z) in the main grid.
-    :type p2: list, non-optional
-    :param ratio: Ratio of the main grid spatial step to the sub-grid spatial step. Must be an odd integer.
-    :type ratio: int, non-optional
-    :param id: Identifier for the sub-grid.
-    :type id: str, non-optional
-    :param is_os_sep: Number of main grid cells between the Inner Surface and the Outer Surface. Defaults to 3.
-    :type is_os_sep: str, optional
-    :param pml_separation: Number of sub-grid cells between the Outer Surface and the PML. Defaults to ratio // 2 + 2
-    :type pml_separation: int, optional
-    :param subgrid_pml_thickness: Thickness of the PML on each of the 6 sides of the sub-grid. Defaults to 6.
-    :type subgrid_pml_thickness: int, optional
-    :param interpolation: Degree of the interpolation scheme used for spatial interpolation of the fields at the Inner Surface. Defaults to Linear
-    :type interpolation: str, optional
-    :param filter: Turn on the 3-pole filter. Increases numerical stability. Defaults to True
-    :type filter: bool, optional
-
+    Attributes:
+        p1: list of the position of the lower left corner of the Inner Surface 
+            (x, y, z) in the main grid.
+        p2: list of the position of the upper right corner of the Inner Surface 
+            (x, y, z) in the main grid.
+        ratio: int of the ratio of the main grid spatial step to the sub-grid 
+                spatial step. Must be an odd integer.
+        id: string identifier for the sub-grid.
+        is_os_sep: int for the number of main grid cells between the Inner 
+                    Surface and the Outer Surface. Defaults to 3.
+        pml_separation: int for the number of sub-grid cells between the Outer 
+                        Surface and the PML. Defaults to ratio // 2 + 2
+        subgrid_pml_thickness: int for the thickness of the PML on each of the 
+                                6 sides of the sub-grid. Defaults to 6.
+        interpolation: string for the degree of the interpolation scheme used 
+                        for spatial interpolation of the fields at the Inner 
+                        Surface. Defaults to Linear.
+        filter: boolean to turn on the 3-pole filter. Increases numerical 
+                stability. Defaults to True.
     """
     def __init__(self,
                  p1=None,
@@ -175,7 +173,7 @@ class SubGridHSG(SubGridBase):
 
         pml_separation = ratio // 2 + 2
 
-        # copy over the optional parameters
+        # Copy over the optional parameters
         kwargs['p1'] = p1
         kwargs['p2'] = p2
         kwargs['ratio'] = ratio