diff --git a/.github/FUNDING.yml b/.github/FUNDING.yml
new file mode 100644
index 00000000..409e993f
--- /dev/null
+++ b/.github/FUNDING.yml
@@ -0,0 +1,13 @@
+# These are supported funding model platforms
+
+github: # Replace with up to 4 GitHub Sponsors-enabled usernames e.g., [user1, user2]
+patreon: # Replace with a single Patreon username
+open_collective: gprMax
+ko_fi: # Replace with a single Ko-fi username
+tidelift: # Replace with a single Tidelift platform-name/package-name e.g., npm/babel
+community_bridge: # Replace with a single Community Bridge project-name e.g., cloud-foundry
+liberapay: # Replace with a single Liberapay username
+issuehunt: # Replace with a single IssueHunt username
+otechie: # Replace with a single Otechie username
+lfx_crowdfunding: # Replace with a single LFX Crowdfunding project-name e.g., cloud-foundry
+custom: # Replace with up to 4 custom sponsorship URLs e.g., ['link1', 'link2']
\ No newline at end of file
diff --git a/examples/antenna_like_GSSI_1500.vti b/examples/antenna_like_GSSI_1500.vti
deleted file mode 100644
index 56e7a075..00000000
Binary files a/examples/antenna_like_GSSI_1500.vti and /dev/null differ
diff --git a/examples/antenna_like_GSSI_1500_fs.h5 b/examples/antenna_like_GSSI_1500_fs.h5
deleted file mode 100644
index 0012b8b2..00000000
Binary files a/examples/antenna_like_GSSI_1500_fs.h5 and /dev/null differ
diff --git a/examples/antenna_like_GSSI_1500_fs.py b/examples/antenna_like_GSSI_1500_fs.py
index 0ee0453b..7ae28dc2 100644
--- a/examples/antenna_like_GSSI_1500_fs.py
+++ b/examples/antenna_like_GSSI_1500_fs.py
@@ -41,7 +41,7 @@ gv2 = gprMax.GeometryView(p1=(ant_pos[0] - 0.170/2, ant_pos[1] - 0.108/2, ant_po
                           p2=(ant_pos[0] + 0.170/2, ant_pos[1] + 0.108/2, ant_pos[2] + 0.010),
                           dl=(dl, dl, dl), filename='antenna_like_GSSI_1500_pcb',
                           output_type='f')
-scene.add(gv1)
+#scene.add(gv1)
 #scene.add(gv2)
 
 # Run model
diff --git a/examples/heterogeneous_soil_v2.h5 b/examples/heterogeneous_soil_v2.h5
deleted file mode 100644
index c0f6519a..00000000
Binary files a/examples/heterogeneous_soil_v2.h5 and /dev/null differ
diff --git a/examples/subgrids/cylinder_fs.py b/examples/subgrids/cylinder_fs.py
index 82d02b88..4800eb60 100644
--- a/examples/subgrids/cylinder_fs.py
+++ b/examples/subgrids/cylinder_fs.py
@@ -89,17 +89,10 @@ gv1 = gprMax.GeometryView(p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl),
 scene.add(gv1)
 
 # Create some snapshots of entire domain
-# for i in range(5):
-#     s = gprMax.Snapshot(p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl),
-#                         time=(i + 0.5) * 1e-9,
-#                         filename=fn.with_suffix('').parts[-1] + '_' + str(i + 1))
-#     scene.add(s)
-
 for i in range(5):
-    s = gprMax.Snapshot(p1=sg1, p2=sg2, dl=(dl_sg, dl_sg, dl_sg),
+    s = gprMax.Snapshot(p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl),
                         time=(i + 0.5) * 1e-9,
-                        filename=fn.with_suffix('').parts[-1] + '_' + str(i + 1),
-                        outputs=['Ez'])
-    subgrid.add(s)
+                        filename=fn.with_suffix('').parts[-1] + '_' + str(i + 1))
+    scene.add(s)
 
 gprMax.run(scenes=[scene], n=1, geometry_only=False, outputfile=fn, subgrid=True, autotranslate=True)
diff --git a/examples/subgrids/gssi_400_over_fractal_subsurface.py b/examples/subgrids/gssi_400_over_fractal_subsurface.py
index 0ed7ae67..9a4546a8 100644
--- a/examples/subgrids/gssi_400_over_fractal_subsurface.py
+++ b/examples/subgrids/gssi_400_over_fractal_subsurface.py
@@ -102,9 +102,11 @@ halfspace = gprMax.Box(p1=(0, 0, 0), p2=(x, y, antenna_p[2]), material_id='soil'
 scene.add(halfspace)
 
 for i in range(1, 51):
-    snap = gprMax.Snapshot(p1=(0, y / 2, 0), p2=(x, y / 2 + dl, z), dl=(dl, dl, dl),
-                           filename=Path(*parts[:-1], parts[-1] + '_' + str(i)).name,
-                           time=i * tw / 50)
+    snap = gprMax.Snapshot(p1=(0, y / 2, 0),
+                           p2=(x, y / 2 + dl, z),
+                           dl=(dl, dl, dl),
+                           filename=Path(*parts[:-1], f'{parts[-1]}_{str(i)}').name,
+                           time=i * tw / 50,)
     scene.add(snap)
 
 # create a geometry view of the main grid and the sub grid stitched together
diff --git a/gprMax/__init__.py b/gprMax/__init__.py
index b5ccc50b..55463eeb 100644
--- a/gprMax/__init__.py
+++ b/gprMax/__init__.py
@@ -25,8 +25,9 @@ from .cmds_geometry.triangle import Triangle
 from .cmds_multiuse import (PMLCFS, AddDebyeDispersion, AddDrudeDispersion,
                             AddLorentzDispersion, GeometryObjectsWrite,
                             GeometryView, HertzianDipole, MagneticDipole,
-                            Material, Rx, RxArray, Snapshot, SoilPeplinski,
-                            TransmissionLine, VoltageSource, Waveform, MaterialRange, MaterialList)
+                            Material, MaterialList, MaterialRange, Rx, RxArray, 
+                            Snapshot, SoilPeplinski, TransmissionLine, 
+                            VoltageSource, Waveform)
 from .cmds_singleuse import (Discretisation, Domain, ExcitationFile,
                              OMPThreads, PMLProps, RxSteps, SrcSteps,
                              TimeStepStabilityFactor, TimeWindow, Title)
diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py
index d932f680..43b80f25 100644
--- a/gprMax/cmds_geometry/add_grass.py
+++ b/gprMax/cmds_geometry/add_grass.py
@@ -74,7 +74,7 @@ class AddGrass(UserObjectGeometry):
             limits = self.kwargs['limits']
             n_blades = self.kwargs['n_blades']
         except KeyError:
-            logger.exception(self.__str__() + ' requires at least eleven parameters')
+            logger.exception(f'{self.__str__()} requires at least eleven parameters')
             raise
 
         try:
@@ -90,7 +90,7 @@ class AddGrass(UserObjectGeometry):
         try:
             volume = volumes[0]
         except NameError:
-            logger.exception(self.__str__() + f' cannot find FractalBox {fractal_box_id}')
+            logger.exception(f'{self.__str__()} cannot find FractalBox {fractal_box_id}')
             raise
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
@@ -98,32 +98,32 @@ class AddGrass(UserObjectGeometry):
         xf, yf, zf = p2
 
         if frac_dim < 0:
-            logger.exception(self.__str__() + ' requires a positive value for ' +
+            logger.exception(f'{self.__str__()} requires a positive value for ' +
                              'the fractal dimension')
             raise ValueError
         if limits[0] < 0 or limits[1] < 0:
-            logger.exception(self.__str__() + ' requires a positive value for ' +
+            logger.exception(f'{self.__str__()} requires a positive value for ' +
                              'the minimum and maximum heights for grass blades')
             raise ValueError
 
         # Check for valid orientations
         if xs == xf:
             if ys == yf or zs == zf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
             if xs != volume.xs and xs != volume.xf:
-                logger.exception(self.__str__() + ' must specify external surfaces on a fractal box')
+                logger.exception(f'{self.__str__()} must specify external surfaces on a fractal box')
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
             # xminus surface
             if xs == volume.xs:
-                logger.exception(self.__str__() + ' grass can only be specified ' +
+                logger.exception(f'{self.__str__()} grass can only be specified ' +
                                  'on surfaces in the positive axis direction')
                 raise ValueError
             # xplus surface
             elif xf == volume.xf:
                 if fractalrange[1] > grid.nx:
-                    logger.exception(self.__str__() + ' cannot apply grass to ' +
+                    logger.exception(f'{self.__str__()} cannot apply grass to ' +
                                      'fractal box as it would exceed the domain ' +
                                      'size in the x direction')
                     raise ValueError
@@ -131,21 +131,21 @@ class AddGrass(UserObjectGeometry):
 
         elif ys == yf:
             if xs == xf or zs == zf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
             if ys != volume.ys and ys != volume.yf:
-                logger.exception(self.__str__() + ' must specify external surfaces on a fractal box')
+                logger.exception(f'{self.__str__()} must specify external surfaces on a fractal box')
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
             # yminus surface
             if ys == volume.ys:
-                logger.exception(self.__str__() + ' grass can only be specified ' +
+                logger.exception(f'{self.__str__()} grass can only be specified ' +
                                  'on surfaces in the positive axis direction')
                 raise ValueError
             # yplus surface
             elif yf == volume.yf:
                 if fractalrange[1] > grid.ny:
-                    logger.exception(self.__str__() + ' cannot apply grass to ' + 
+                    logger.exception(f'{self.__str__()} cannot apply grass to ' + 
                                      'fractal box as it would exceed the domain ' +
                                      'size in the y direction')
                     raise ValueError
@@ -153,28 +153,28 @@ class AddGrass(UserObjectGeometry):
 
         elif zs == zf:
             if xs == xf or ys == yf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
             if zs != volume.zs and zs != volume.zf:
-                logger.exception(self.__str__() + ' must specify external surfaces on a fractal box')
+                logger.exception(f'{self.__str__()} must specify external surfaces on a fractal box')
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
             # zminus surface
             if zs == volume.zs:
-                logger.exception(self.__str__() + ' grass can only be specified ' +
+                logger.exception(f'{self.__str__()} grass can only be specified ' +
                                  'on surfaces in the positive axis direction')
                 raise ValueError
             # zplus surface
             elif zf == volume.zf:
                 if fractalrange[1] > grid.nz:
-                    logger.exception(self.__str__() + ' cannot apply grass to ' +
+                    logger.exception(f'{self.__str__()} cannot apply grass to ' +
                                      'fractal box as it would exceed the domain ' +
                                      'size in the z direction')
                     raise ValueError
                 requestedsurface = 'zplus'
 
         else:
-            logger.exception(self.__str__() + ' dimensions are not specified correctly')
+            logger.exception(f'{self.__str__()} dimensions are not specified correctly')
             raise ValueError
 
         surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim)
@@ -187,7 +187,7 @@ class AddGrass(UserObjectGeometry):
         surface.operatingonID = volume.ID
         surface.generate_fractal_surface()
         if n_blades > surface.fractalsurface.shape[0] * surface.fractalsurface.shape[1]:
-            logger.exception(self.__str__() + ' the specified surface is not large ' +
+            logger.exception(f'{self.__str__()} the specified surface is not large ' +
                              'enough for the number of grass blades/roots specified')
             raise ValueError
 
@@ -234,13 +234,13 @@ class AddGrass(UserObjectGeometry):
         grass = next((x for x in grid.materials if x.ID == 'grass'))
         testgrass = next((x for x in grass.tau if x < grid.dt), None)
         if testgrass:
-            logger.exception(self.__str__() + ' requires the time step for the ' +
+            logger.exception(f'{self.__str__()} requires the time step for the ' +
                              'model to be less than the relaxation time required to model grass.')
             raise ValueError
 
         volume.fractalsurfaces.append(surface)
 
-        logger.info(self.grid_name(grid) + f'{n_blades} blades of grass on surface from ' + 
+        logger.info(f'{self.grid_name(grid)}{n_blades} blades of grass on surface from ' + 
                     f'{xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, ' +
                     f'to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m ' +
                     f'with fractal dimension {surface.dimension:g}, fractal seeding ' +
diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py
index 0a9f0443..df4b02da 100644
--- a/gprMax/cmds_geometry/add_surface_roughness.py
+++ b/gprMax/cmds_geometry/add_surface_roughness.py
@@ -74,13 +74,13 @@ class AddSurfaceRoughness(UserObjectGeometry):
             limits = np.array(self.kwargs['limits'])
             fractal_box_id = self.kwargs['fractal_box_id']
         except KeyError:
-            logger.exception(self.__str__() + ' incorrect parameters')
+            logger.exception(f'{self.__str__()} incorrect parameters')
             raise
 
         try:
             seed = self.kwargs['seed']
         except KeyError:
-            logger.warning(self.__str__() + ' no value for seed detected. This ' +
+            logger.warning(f'{self.__str__()} no value for seed detected. This ' +
                            'means you will get a different fractal distribution ' +
                            'every time the model runs.')
             seed = None
@@ -93,7 +93,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
         if volumes:
             volume = volumes[0]
         else:
-            logger.exception(self.__str__() + f' cannot find FractalBox {fractal_box_id}')
+            logger.exception(f'{self.__str__()} cannot find FractalBox {fractal_box_id}')
             raise ValueError
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
@@ -101,25 +101,25 @@ class AddSurfaceRoughness(UserObjectGeometry):
         xf, yf, zf = p2
 
         if frac_dim < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'fractal dimension')
             raise ValueError
         if weighting[0] < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'fractal weighting in the first direction of the surface')
             raise ValueError
         if weighting[1] < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'fractal weighting in the second direction of the surface')
             raise ValueError
 
         # Check for valid orientations
         if xs == xf:
             if ys == yf or zs == zf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
             if xs != volume.xs and xs != volume.xf:
-                logger.exception(self.__str__() + ' can only be used on the external ' + 
+                logger.exception(f'{self.__str__()} can only be used on the external ' + 
                                  'surfaces of a fractal box')
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dx), 
@@ -127,7 +127,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             # xminus surface
             if xs == volume.xs:
                 if fractalrange[0] < 0 or fractalrange[1] > volume.xf:
-                    logger.exception(self.__str__() + ' cannot apply fractal surface ' +
+                    logger.exception(f'{self.__str__()} cannot apply fractal surface ' +
                                      'to fractal box as it would exceed either the ' +
                                      'upper coordinates of the fractal box or the ' +
                                      'domain in the x direction')
@@ -136,7 +136,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             # xplus surface
             elif xf == volume.xf:
                 if fractalrange[0] < volume.xs or fractalrange[1] > grid.nx:
-                    logger.exception(self.__str__() + ' cannot apply fractal surface ' +
+                    logger.exception(f'{self.__str__()} cannot apply fractal surface ' +
                                      'to fractal box as it would exceed either the ' +
                                      'lower coordinates of the fractal box or the ' +
                                      'domain in the x direction')
@@ -145,10 +145,10 @@ class AddSurfaceRoughness(UserObjectGeometry):
 
         elif ys == yf:
             if xs == xf or zs == zf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
             if ys != volume.ys and ys != volume.yf:
-                logger.exception(self.__str__() + ' can only be used on the external ' + 
+                logger.exception(f'{self.__str__()} can only be used on the external ' + 
                                  'surfaces of a fractal box')
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dy), 
@@ -156,7 +156,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             # yminus surface
             if ys == volume.ys:
                 if fractalrange[0] < 0 or fractalrange[1] > volume.yf:
-                    logger.exception(self.__str__() + ' cannot apply fractal surface ' +
+                    logger.exception(f'{self.__str__()} cannot apply fractal surface ' +
                                      'to fractal box as it would exceed either the ' +
                                      'upper coordinates of the fractal box or the ' +
                                      'domain in the y direction')
@@ -165,7 +165,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             # yplus surface
             elif yf == volume.yf:
                 if fractalrange[0] < volume.ys or fractalrange[1] > grid.ny:
-                    logger.exception(self.__str__() + ' cannot apply fractal surface ' +
+                    logger.exception(f'{self.__str__()} cannot apply fractal surface ' +
                                      'to fractal box as it would exceed either the ' +
                                      'lower coordinates of the fractal box or the ' +
                                      'domain in the y direction')
@@ -174,10 +174,10 @@ class AddSurfaceRoughness(UserObjectGeometry):
 
         elif zs == zf:
             if xs == xf or ys == yf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
             if zs != volume.zs and zs != volume.zf:
-                logger.exception(self.__str__() + ' can only be used on the external ' +
+                logger.exception(f'{self.__str__()} can only be used on the external ' +
                                  'surfaces of a fractal box')
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dz), 
@@ -185,7 +185,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             # zminus surface
             if zs == volume.zs:
                 if fractalrange[0] < 0 or fractalrange[1] > volume.zf:
-                    logger.exception(self.__str__() + ' cannot apply fractal surface ' +
+                    logger.exception(f'{self.__str__()} cannot apply fractal surface ' +
                                      'to fractal box as it would exceed either the ' +
                                      'upper coordinates of the fractal box or the ' +
                                      'domain in the x direction')
@@ -194,7 +194,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             # zplus surface
             elif zf == volume.zf:
                 if fractalrange[0] < volume.zs or fractalrange[1] > grid.nz:
-                    logger.exception(self.__str__() + ' cannot apply fractal surface ' +
+                    logger.exception(f'{self.__str__()} cannot apply fractal surface ' +
                                      'to fractal box as it would exceed either the ' +
                                      'lower coordinates of the fractal box or the ' +
                                      'domain in the z direction')
@@ -202,7 +202,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
                 requestedsurface = 'zplus'
 
         else:
-            logger.exception(self.__str__() + ' dimensions are not specified correctly')
+            logger.exception(f'{self.__str__()} dimensions are not specified correctly')
             raise ValueError
 
         surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim)
@@ -218,14 +218,14 @@ 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 ' +
+            logger.exception(f'{self.__str__()} has already been used on the ' +
                              f'{surface.surfaceID} surface')
             raise ValueError
 
         surface.generate_fractal_surface()
         volume.fractalsurfaces.append(surface)
 
-        logger.info(self.grid_name(grid) + f'Fractal surface from {xs * grid.dx:g}m, ' +
+        logger.info(f'{self.grid_name(grid)}Fractal surface from {xs * grid.dx:g}m, ' +
                     f'{ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, ' +
                     f'{yf * grid.dy:g}m, {zf * grid.dz:g}m with fractal dimension ' +
                     f'{surface.dimension:g}, fractal weightings {surface.weighting[0]:g}, ' +
diff --git a/gprMax/cmds_geometry/add_surface_water.py b/gprMax/cmds_geometry/add_surface_water.py
index 4e8960ea..00d33a48 100644
--- a/gprMax/cmds_geometry/add_surface_water.py
+++ b/gprMax/cmds_geometry/add_surface_water.py
@@ -68,18 +68,16 @@ class AddSurfaceWater(UserObjectGeometry):
             fractal_box_id = self.kwargs['fractal_box_id']
             depth = self.kwargs['depth']
         except KeyError:
-            logger.exception(self.__str__() + ' requires exactly eight parameters')
+            logger.exception(f'{self.__str__()} requires exactly eight parameters')
             raise
 
         if self.do_rotate:
             self._do_rotate()
 
-        # Get the correct fractal volume
-        volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
-        if volumes:
+        if volumes := [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]:
             volume = volumes[0]
         else:
-            logger.exception(self.__str__() + f' cannot find FractalBox {fractal_box_id}')
+            logger.exception(f'{self.__str__()} cannot find FractalBox {fractal_box_id}')
             raise ValueError
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
@@ -87,18 +85,18 @@ 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(f'{self.__str__()} requires a positive value for the ' + 
+                             f'depth of water')
             raise ValueError
 
         # Check for valid orientations
         if xs == xf:
             if ys == yf or zs == zf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
-            if xs != volume.xs and xs != volume.xf:
-                logger.exception(self.__str__() + ' can only be used on the ' +
-                                 'external surfaces of a fractal box')
+            if xs not in [volume.xs, volume.xf]:
+                logger.exception(f'{self.__str__()} can only be used on the external surfaces '
+                                 f'of a fractal box')
                 raise ValueError
             # xminus surface
             if xs == volume.xs:
@@ -110,12 +108,12 @@ class AddSurfaceWater(UserObjectGeometry):
             filldepth = filldepthcells * grid.dx
 
         elif ys == yf:
-            if xs == xf or zs == zf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
+            if zs == zf:
+                logger.exception(f'{self.__str__()} dimensions are not specified correctly')
                 raise ValueError
-            if ys != volume.ys and ys != volume.yf:
-                logger.exception(self.__str__() + ' can only be used on the ' +
-                                 'external surfaces of a fractal box')
+            if ys not in [volume.ys, volume.yf]:
+                logger.exception(f'{self.__str__()} can only be used on the external surfaces ' +
+                                 f'of a fractal box')
                 raise ValueError
             # yminus surface
             if ys == volume.ys:
@@ -127,12 +125,9 @@ class AddSurfaceWater(UserObjectGeometry):
             filldepth = filldepthcells * grid.dy
 
         elif zs == zf:
-            if xs == xf or ys == yf:
-                logger.exception(self.__str__() + ' dimensions are not specified correctly')
-                raise ValueError
-            if zs != volume.zs and zs != volume.zf:
-                logger.exception(self.__str__() + ' can only be used on the ' +
-                                 'external surfaces of a fractal box')
+            if zs not in [volume.zs, volume.zf]:
+                logger.exception(f'{self.__str__()} can only be used on the external surfaces '
+                                f'of a fractal box')
                 raise ValueError
             # zminus surface
             if zs == volume.zs:
@@ -144,38 +139,35 @@ class AddSurfaceWater(UserObjectGeometry):
             filldepth = filldepthcells * grid.dz
 
         else:
-            logger.exception(self.__str__() + ' dimensions are not specified correctly')
+            logger.exception(f'{self.__str__()} dimensions are not specified correctly')
             raise ValueError
 
         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(f'{self.__str__()} specified surface {requestedsurface} ' +
+                             f'does not have a rough surface applied')
             raise ValueError
 
         surface.filldepth = filldepthcells
 
         # Check that requested fill depth falls within range of surface roughness
         if surface.filldepth < surface.fractalrange[0] or surface.filldepth > surface.fractalrange[1]:
-            logger.exception(self.__str__() + ' requires a value for the depth ' +
-                             'of water that lies with the range of the requested ' +
-                             'surface roughness')
+            logger.exception(f'{self.__str__()} requires a value for the depth of water that lies with the ' +
+                             f'range of the requested surface roughness')
             raise ValueError
 
         # Check to see if water has been already defined as a material
-        if not any(x.ID == 'water' for x in grid.materials):
+        if all(x.ID != 'water' for x in grid.materials):
             create_water(grid)
 
         # Check if time step for model is suitable for using water
         water = next((x for x in grid.materials if x.ID == 'water'))
-        testwater = next((x for x in water.tau if x < grid.dt), None)
-        if testwater:
-            logger.exception(self.__str__() + ' requires the time step for the ' +
-                             'model to be less than the relaxation time required ' +
-                             'to model water.')
+        if testwater := next((x for x in water.tau if x < grid.dt), None):
+            logger.exception(f'{self.__str__()} requires the time step for the model '
+                             f'to be less than the relaxation time required to model water.')
             raise ValueError
 
-        logger.info(self.grid_name(grid) + f'Water on surface from {xs * grid.dx:g}m, ' +
+        logger.info(f'{self.grid_name(grid)}Water on surface from {xs * grid.dx:g}m, ' +
                     f'{ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, ' +
                     f'{yf * grid.dy:g}m, {zf * grid.dz:g}m with depth {filldepth:g}m, ' +
                     f'added to {surface.operatingonID}.')
diff --git a/gprMax/cmds_geometry/box.py b/gprMax/cmds_geometry/box.py
index aa509bf7..de07b862 100644
--- a/gprMax/cmds_geometry/box.py
+++ b/gprMax/cmds_geometry/box.py
@@ -65,7 +65,7 @@ class Box(UserObjectGeometry):
             p1 = self.kwargs['p1']
             p2 = self.kwargs['p2']
         except KeyError:
-            logger.exception(self.__str__() + ' Please specify two points.')
+            logger.exception(f'{self.__str__()} Please specify two points.')
             raise
 
         if self.do_rotate:
@@ -80,7 +80,7 @@ class Box(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' No materials have been specified')
+                logger.exception(f'{self.__str__()} No materials have been specified')
                 raise
 
         # Check averaging
@@ -103,7 +103,7 @@ class Box(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         # Isotropic case
@@ -144,7 +144,7 @@ class Box(UserObjectGeometry):
 
         dielectricsmoothing = 'on' if averaging else 'off'
 
-        logger.info(self.grid_name(grid) + f"Box from {p5[0]:g}m, {p5[1]:g}m, " +
+        logger.info(f"{self.grid_name(grid)}Box from {p5[0]:g}m, {p5[1]:g}m, " +
                     f"{p5[2]:g}m, to {p6[0]:g}m, {p6[1]:g}m, {p6[2]:g}m of " +
                     f"material(s) {', '.join(materialsrequested)} created, " +
                     f"dielectric smoothing is {dielectricsmoothing}.")
diff --git a/gprMax/cmds_geometry/build_templates.py b/gprMax/cmds_geometry/build_templates.py
index b78543e6..7f0cbd2b 100644
--- a/gprMax/cmds_geometry/build_templates.py
+++ b/gprMax/cmds_geometry/build_templates.py
@@ -66,6 +66,5 @@ r = template.render(
         }]
 )
 
-f = open('cython/dispersive_updates_test.pyx', 'w')
-f.write(r)
-f.close()
+with open('cython/dispersive_updates_test.pyx', 'w') as f:
+    f.write(r)
diff --git a/gprMax/cmds_geometry/cmds_geometry.py b/gprMax/cmds_geometry/cmds_geometry.py
index 0ad5bc3d..63980390 100644
--- a/gprMax/cmds_geometry/cmds_geometry.py
+++ b/gprMax/cmds_geometry/cmds_geometry.py
@@ -38,9 +38,9 @@ class UserObjectGeometry:
         """Readable string of parameters given to object."""
         s = ''
         for _, v in self.kwargs.items():
-            if isinstance(v, tuple) or isinstance(v, list):
+            if isinstance(v, (tuple, list)):
                 v = ' '.join([str(el) for el in v])
-            s += str(v) + ' '
+            s += f'{str(v)} '
 
         return f'{self.hash}: {s[:-1]}'
 
@@ -109,7 +109,7 @@ def rotate_2point_object(pts, axis, angle, origin=None):
     # Use origin at centre of object if not given
     if not origin:
         origin = pts[0,:] + (pts[1,:] - pts[0,:]) / 2
-    
+
     # Check angle value is suitable
     angle = int(angle)
     if angle < 0 or angle > 360:
@@ -120,7 +120,7 @@ def rotate_2point_object(pts, axis, angle, origin=None):
         raise ValueError
 
     # Check axis is valid
-    if axis != 'x' and axis != 'y' and axis != 'z':
+    if axis not in ['x', 'y', 'z']:
         logger.exception('Axis of rotation must be x, y, or z')
         raise ValueError
 
diff --git a/gprMax/cmds_geometry/cone.py b/gprMax/cmds_geometry/cone.py
index 34e5ec13..609505a1 100644
--- a/gprMax/cmds_geometry/cone.py
+++ b/gprMax/cmds_geometry/cone.py
@@ -55,7 +55,7 @@ class Cone(UserObjectGeometry):
             r1 = self.kwargs['r1']
             r2 = self.kwargs['r2']
         except KeyError:
-            logger.exception(self.__str__() + ' please specify 2 points and two radii')
+            logger.exception(f'{self.__str__()} please specify 2 points and two radii')
             raise
 
         # Check averaging
@@ -75,7 +75,7 @@ class Cone(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' no materials have been specified')
+                logger.exception(f'{self.__str__()} no materials have been specified')
                 raise
 
         p3 = uip.round_to_grid_static_point(p1)
@@ -85,15 +85,15 @@ class Cone(UserObjectGeometry):
         x2, y2, z2 = uip.round_to_grid(p2)
 
         if r1 < 0:
-            logger.exception(self.__str__() + f' the radius of the first face {r1:g} should be a positive value.')
+            logger.exception(f'{self.__str__()} the radius of the first face {r1:g} should be a positive value.')
             raise ValueError
         
         if r2 < 0:
-            logger.exception(self.__str__() + f' the radius of the second face {r2:g} should be a positive value.')
+            logger.exception(f'{self.__str__()} the radius of the second face {r2:g} should be a positive value.')
             raise ValueError
         
         if r1 == 0 and r2 == 0:
-            logger.exception(self.__str__() + f' not both radii can be zero.')
+            logger.exception(f'{self.__str__()} not both radii can be zero.')
             raise ValueError
 
         # Look up requested materials in existing list of material instances
@@ -101,7 +101,7 @@ class Cone(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         # Isotropic case
@@ -141,7 +141,7 @@ class Cone(UserObjectGeometry):
                        grid.rigidE, grid.rigidH, grid.ID)
 
         dielectricsmoothing = 'on' if averaging else 'off'
-        logger.info(self.grid_name(grid) + f"Cone with face centres {p3[0]:g}m, " +
+        logger.info(f"{self.grid_name(grid)}Cone with face centres {p3[0]:g}m, " +
                     f"{p3[1]:g}m, {p3[2]:g}m and {p4[0]:g}m, {p4[1]:g}m, {p4[2]:g}m, " +
                     f"with radii {r1:g}m and {r2:g}, of material(s) {', '.join(materialsrequested)} " +
                     f"created, dielectric smoothing is {dielectricsmoothing}.")
diff --git a/gprMax/cmds_geometry/cylinder.py b/gprMax/cmds_geometry/cylinder.py
index d70b824a..006b46a4 100644
--- a/gprMax/cmds_geometry/cylinder.py
+++ b/gprMax/cmds_geometry/cylinder.py
@@ -52,7 +52,7 @@ class Cylinder(UserObjectGeometry):
             p2 = self.kwargs['p2']
             r = self.kwargs['r']
         except KeyError:
-            logger.exception(self.__str__() + ' please specify 2 points and a radius')
+            logger.exception(f'{self.__str__()} please specify 2 points and a radius')
             raise
 
         # Check averaging
@@ -72,7 +72,7 @@ class Cylinder(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' no materials have been specified')
+                logger.exception(f'{self.__str__()} no materials have been specified')
                 raise
 
         p3 = uip.round_to_grid_static_point(p1)
@@ -82,7 +82,7 @@ class Cylinder(UserObjectGeometry):
         x2, y2, z2 = uip.round_to_grid(p2)
 
         if r <= 0:
-            logger.exception(self.__str__() + f' the radius {r:g} should be a positive value.')
+            logger.exception(f'{self.__str__()} the radius {r:g} should be a positive value.')
             raise ValueError
 
         # Look up requested materials in existing list of material instances
@@ -90,7 +90,7 @@ class Cylinder(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         # Isotropic case
@@ -130,7 +130,7 @@ class Cylinder(UserObjectGeometry):
                        grid.rigidE, grid.rigidH, grid.ID)
 
         dielectricsmoothing = 'on' if averaging else 'off'
-        logger.info(self.grid_name(grid) + f"Cylinder with face centres {p3[0]:g}m, " +
+        logger.info(f"{self.grid_name(grid)}Cylinder with face centres {p3[0]:g}m, " +
                     f"{p3[1]:g}m, {p3[2]:g}m and {p4[0]:g}m, {p4[1]:g}m, {p4[2]:g}m, " +
                     f"with radius {r:g}m, of material(s) {', '.join(materialsrequested)} " +
                     f"created, dielectric smoothing is {dielectricsmoothing}.")
diff --git a/gprMax/cmds_geometry/cylindrical_sector.py b/gprMax/cmds_geometry/cylindrical_sector.py
index 42005250..8fa52c6c 100644
--- a/gprMax/cmds_geometry/cylindrical_sector.py
+++ b/gprMax/cmds_geometry/cylindrical_sector.py
@@ -88,24 +88,25 @@ class CylindricalSector(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' No materials have been specified')
+                logger.exception(f'{self.__str__()} No materials have been specified')
                 raise
 
         sectorstartangle = 2 * np.pi * (start / 360)
         sectorangle = 2 * np.pi * (end / 360)
 
-        if normal != 'x' and normal != 'y' and normal != 'z':
-            logger.exception(self.__str__() + ' the normal direction must be either x, y or z.')
+        if normal not in ['x', 'y', 'z']:
+            logger.exception(f'{self.__str__()} the normal direction must be either ' +
+                             f'x, y or z.')
             raise ValueError
         if r <= 0:
-            logger.exception(self.__str__() + f' the radius {r:g} should be a positive value.')
+            logger.exception(f'{self.__str__()} the radius {r:g} should be a positive value.')
         if sectorstartangle < 0 or sectorangle <= 0:
-            logger.exception(self.__str__() + ' the starting angle and sector ' +
-                             'angle should be a positive values.')
+            logger.exception(f'{self.__str__()} the starting angle and sector angle should be ' +
+                             f'a positive values.')
             raise ValueError
         if sectorstartangle >= 2 * np.pi or sectorangle >= 2 * np.pi:
-            logger.exception(self.__str__() + ' the starting angle and sector ' +
-                             'angle must be less than 360 degrees.')
+            logger.exception(f'{self.__str__()} the starting angle and sector angle must be ' +
+                             f'less than 360 degrees.')
             raise ValueError
 
         # Look up requested materials in existing list of material instances
@@ -113,7 +114,7 @@ class CylindricalSector(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         if thickness > 0:
@@ -122,13 +123,12 @@ class CylindricalSector(UserObjectGeometry):
                 averaging = materials[0].averagable and averagecylindricalsector
                 numID = numIDx = numIDy = numIDz = materials[0].numID
 
-            # Uniaxial anisotropic case
             elif len(materials) == 3:
                 averaging = False
                 numIDx = materials[0].numID
                 numIDy = materials[1].numID
                 numIDz = materials[2].numID
-                requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+                requiredID = f'{materials[0].ID}+{materials[1].ID}+{materials[2].ID}'
                 averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
                 if averagedmaterial:
                     numID = averagedmaterial.numID
@@ -181,14 +181,14 @@ class CylindricalSector(UserObjectGeometry):
 
         if thickness > 0:
             dielectricsmoothing = 'on' if averaging else 'off'
-            logger.info(self.grid_name(grid) + f"Cylindrical sector with centre " +
+            logger.info(f"{self.grid_name(grid)}Cylindrical sector with centre " +
                         f"{ctr1:g}m, {ctr2:g}m, radius {r:g}m, starting angle " +
                         f"{(sectorstartangle / (2 * np.pi)) * 360:.1f} degrees, " +
                         f"sector angle {(sectorangle / (2 * np.pi)) * 360:.1f} degrees, " +
                         f"thickness {thickness:g}m, of material(s) {', '.join(materialsrequested)} " +
                         f"created, dielectric smoothing is {dielectricsmoothing}.")
         else:
-            logger.info(self.grid_name(grid) + f"Cylindrical sector with centre " +
+            logger.info(f"{self.grid_name(grid)}Cylindrical sector with centre " +
                         f"{ctr1:g}m, {ctr2:g}m, radius {r:g}m, starting angle " +
                         f"{(sectorstartangle / (2 * np.pi)) * 360:.1f} degrees, " +
                         f"sector angle {(sectorangle / (2 * np.pi)) * 360:.1f} " +
diff --git a/gprMax/cmds_geometry/edge.py b/gprMax/cmds_geometry/edge.py
index acba61db..e1af273e 100644
--- a/gprMax/cmds_geometry/edge.py
+++ b/gprMax/cmds_geometry/edge.py
@@ -62,15 +62,15 @@ class Edge(UserObjectGeometry):
             p2 = self.kwargs['p2']
             material_id = self.kwargs['material_id']
         except KeyError:
-            logger.exception(self.__str__() + ' requires exactly 3 parameters')
+            logger.exception(f'{self.__str__()} requires exactly 3 parameters')
             raise
 
         if self.do_rotate:
             self._do_rotate()
-        
+
         p3 = uip.round_to_grid_static_point(p1)
         p4 = uip.round_to_grid_static_point(p2)
-        
+
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
         xs, ys, zs = p1
         xf, yf, zf = p2
@@ -83,32 +83,23 @@ class Edge(UserObjectGeometry):
 
         # Check for valid orientations
         # x-orientated edge
-        if xs != xf:
-            if ys != yf or zs != zf:
-                logger.exception(self.__str__() + ' the edge is not specified correctly')
-                raise ValueError
-            else:
-                for i in range(xs, xf):
-                    build_edge_x(i, ys, zs, material.numID, grid.rigidE, grid.rigidH, grid.ID)
+        if ((xs != xf and (ys != yf or zs != zf))
+            or (ys != yf and (xs != xf or zs != zf))
+            or (zs != zf and (xs != xf or ys != yf))):
+            logger.exception(f'{self.__str__()} the edge is not specified correctly')
+            raise ValueError
+        elif xs != xf:
+            for i in range(xs, xf):
+                build_edge_x(i, ys, zs, material.numID, grid.rigidE, grid.rigidH, grid.ID)
 
-        # y-orientated edge
         elif ys != yf:
-            if xs != xf or zs != zf:
-                logger.exception(self.__str__() + ' the edge is not specified correctly')
-                raise ValueError
-            else:
-                for j in range(ys, yf):
-                    build_edge_y(xs, j, zs, material.numID, grid.rigidE, grid.rigidH, grid.ID)
+            for j in range(ys, yf):
+                build_edge_y(xs, j, zs, material.numID, grid.rigidE, grid.rigidH, grid.ID)
 
-        # z-orientated edge
         elif zs != zf:
-            if xs != xf or ys != yf:
-                logger.exception(self.__str__() + ' the edge is not specified correctly')
-                raise ValueError
-            else:
-                for k in range(zs, zf):
-                    build_edge_z(xs, ys, k, material.numID, grid.rigidE, grid.rigidH, grid.ID)
+            for k in range(zs, zf):
+                build_edge_z(xs, ys, k, material.numID, grid.rigidE, grid.rigidH, grid.ID)
 
-        logger.info(self.grid_name(grid) + f'Edge from {p3[0]:g}m, {p3[1]:g}m, ' +
+        logger.info(f'{self.grid_name(grid)}Edge from {p3[0]:g}m, {p3[1]:g}m, ' +
                     f'{p3[2]:g}m, to {p4[0]:g}m, {p4[1]:g}m, {p4[2]:g}m of ' +
                     f'material {material_id} created.')
diff --git a/gprMax/cmds_geometry/ellipsoid.py b/gprMax/cmds_geometry/ellipsoid.py
index 115ee10c..ba5e7264 100644
--- a/gprMax/cmds_geometry/ellipsoid.py
+++ b/gprMax/cmds_geometry/ellipsoid.py
@@ -53,7 +53,7 @@ class Ellipsoid(UserObjectGeometry):
             zr = self.kwargs['zr']
 
         except KeyError:
-            logger.exception(self.__str__() + ' please specify a point and the three semiaxes.')
+            logger.exception(f'{self.__str__()} please specify a point and the three semiaxes.')
             raise
 
         # Check averaging
@@ -73,7 +73,7 @@ class Ellipsoid(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' no materials have been specified')
+                logger.exception(f'{self.__str__()} no materials have been specified')
                 raise
 
         # Centre of sphere
@@ -86,7 +86,7 @@ class Ellipsoid(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         # Isotropic case
@@ -126,7 +126,7 @@ class Ellipsoid(UserObjectGeometry):
                      grid.rigidE, grid.rigidH, grid.ID)
 
         dielectricsmoothing = 'on' if averaging else 'off'
-        logger.info(self.grid_name(grid) + f"Ellipsoid with centre {p2[0]:g}m, " +
+        logger.info(f"{self.grid_name(grid)}Ellipsoid with centre {p2[0]:g}m, " +
                     f"{p2[1]:g}m, {p2[2]:g}m, x-semiaxis {xr:g}m, y-semiaxis {yr:g}m and z-semiaxis {zr:g}m of material(s) " +
                     f"{', '.join(materialsrequested)} created, dielectric " +
                     f"smoothing is {dielectricsmoothing}.")
diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py
index 55d18277..37306ec7 100644
--- a/gprMax/cmds_geometry/fractal_box.py
+++ b/gprMax/cmds_geometry/fractal_box.py
@@ -78,13 +78,13 @@ class FractalBox(UserObjectGeometry):
             mixing_model_id = self.kwargs['mixing_model_id']
             ID = self.kwargs['id']
         except KeyError:
-            logger.exception(self.__str__() + ' Incorrect parameters')
+            logger.exception(f'{self.__str__()} Incorrect parameters')
             raise
 
         try:
             seed = self.kwargs['seed']
         except KeyError:
-            logger.warning(self.__str__() + ' no value for seed detected. This ' +
+            logger.warning(f'{self.__str__()} no value for seed detected. This ' +
                            'means you will get a different fractal distribution ' +
                            'every time the model runs.')
             seed = None
@@ -107,24 +107,24 @@ class FractalBox(UserObjectGeometry):
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
         xs, ys, zs = p1
         xf, yf, zf = p2
-        
+
         if frac_dim < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'fractal dimension')
             raise ValueError
         if weighting[0] < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'fractal weighting in the x direction')
             raise ValueError
         if weighting[1] < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'fractal weighting in the y direction')
             raise ValueError
         if weighting[2] < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'fractal weighting in the z direction')
         if n_materials < 0:
-            logger.exception(self.__str__() + ' requires a positive value for the ' +
+            logger.exception(f'{self.__str__()} requires a positive value for the ' +
                              'number of bins')
             raise ValueError
 
@@ -136,14 +136,14 @@ class FractalBox(UserObjectGeometry):
 
         if mixingmodel:
             if nbins == 1:
-                logger.exception(self.__str__() + ' must be used with more than ' +
+                logger.exception(f'{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.
             mixingmodel.calculate_properties(nbins, grid)
         elif not material:
-            logger.exception(self.__str__() + f' mixing model or material with ' +
+            logger.exception(f'{self.__str__()} mixing model or material with ' +
                              'ID {mixing_model_id} does not exist')
             raise ValueError
 
@@ -160,7 +160,7 @@ class FractalBox(UserObjectGeometry):
         volume.mixingmodel = mixingmodel
 
         dielectricsmoothing = 'on' if volume.averaging else 'off'
-        logger.info(self.grid_name(grid) + f'Fractal box {volume.ID} from ' +
+        logger.info(f'{self.grid_name(grid)}Fractal box {volume.ID} from ' +
                     f'{p3[0]:g}m, {p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, ' +
                     f'{p4[1]:g}m, {p4[2]:g}m with {volume.operatingonID}, ' +
                     f'fractal dimension {volume.dimension:g}, fractal weightings ' +
diff --git a/gprMax/cmds_geometry/fractal_box_builder.py b/gprMax/cmds_geometry/fractal_box_builder.py
index 5c804f71..f2243fa5 100644
--- a/gprMax/cmds_geometry/fractal_box_builder.py
+++ b/gprMax/cmds_geometry/fractal_box_builder.py
@@ -319,7 +319,7 @@ class FractalBoxBuilder(UserObjectGeometry):
 
             else:
                 if volume.nbins == 1:
-                    logger.exception(self.__str__() + ' is being used with a ' +
+                    logger.exception(f'{self.__str__()} is being used with a ' +
                                      'single material and no modifications, ' +
                                      'therefore please use a #box command instead.')
                     raise ValueError
diff --git a/gprMax/cmds_geometry/geometry_objects_read.py b/gprMax/cmds_geometry/geometry_objects_read.py
index 936d7349..a53f84af 100644
--- a/gprMax/cmds_geometry/geometry_objects_read.py
+++ b/gprMax/cmds_geometry/geometry_objects_read.py
@@ -46,7 +46,7 @@ class GeometryObjectsRead(UserObjectGeometry):
             geofile = self.kwargs['geofile']
             matfile = self.kwargs['matfile']
         except KeyError:
-            logger.exception(self.__str__() + 'requires exactly five parameters')
+            logger.exception(f'{self.__str__()} requires exactly five parameters')
             raise
 
         # Discretise the point using uip object. This has different behaviour
@@ -100,7 +100,7 @@ class GeometryObjectsRead(UserObjectGeometry):
         if round_value((dx_dy_dz[0] / grid.dx) != 1 or 
                        round_value(dx_dy_dz[1] / grid.dy) != 1 or 
                        round_value(dx_dy_dz[2] / grid.dz) != 1):
-            logger.exception(self.__str__() + ' requires the spatial resolution ' +
+            logger.exception(f'{self.__str__()} requires the spatial resolution ' +
                              'of the geometry objects file to match the spatial ' +
                              'resolution of the model')
             raise ValueError
@@ -122,7 +122,7 @@ class GeometryObjectsRead(UserObjectGeometry):
             grid.rigidE[:, xs:xs + rigidE.shape[1], ys:ys + rigidE.shape[2], zs:zs + rigidE.shape[3]] = rigidE
             grid.rigidH[:, xs:xs + rigidH.shape[1], ys:ys + rigidH.shape[2], zs:zs + rigidH.shape[3]] = rigidH
             grid.ID[:, xs:xs + ID.shape[1], ys:ys + ID.shape[2], zs:zs + ID.shape[3]] = ID + numexistmaterials
-            logger.info(self.grid_name(grid) + f'Geometry objects from file {geofile} ' +
+            logger.info(f'{self.grid_name(grid)}Geometry objects from file {geofile} ' +
                         f'inserted at {xs * grid.dx:g}m, {ys * grid.dy:g}m, ' +
                         f'{zs * grid.dz:g}m, with corresponding materials file ' +
                         f'{matfile}.')
@@ -131,7 +131,7 @@ class GeometryObjectsRead(UserObjectGeometry):
             build_voxels_from_array(xs, ys, zs, config.get_model_config().ompthreads,
                                     numexistmaterials, averaging, data, 
                                     grid.solid, grid.rigidE, grid.rigidH, grid.ID)
-            logger.info(self.grid_name(grid) + f'Geometry objects from file ' +
+            logger.info(f'{self.grid_name(grid)}Geometry objects from file ' +
                         f'(voxels only){geofile} inserted at {xs * grid.dx:g}m, ' +
                         f'{ys * grid.dy:g}m, {zs * grid.dz:g}m, with corresponding ' +
                         f'materials file {matfile}.')
diff --git a/gprMax/cmds_geometry/plate.py b/gprMax/cmds_geometry/plate.py
index cd157604..9a4b9313 100644
--- a/gprMax/cmds_geometry/plate.py
+++ b/gprMax/cmds_geometry/plate.py
@@ -61,7 +61,7 @@ class Plate(UserObjectGeometry):
             p1 = self.kwargs['p1']
             p2 = self.kwargs['p2']
         except KeyError:
-            logger.exception(self.__str__() + ' 2 points must be specified')
+            logger.exception(f'{self.__str__()} 2 points must be specified')
             raise
 
         # isotropic
@@ -72,7 +72,7 @@ class Plate(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' No materials have been specified')
+                logger.exception(f'{self.__str__()} No materials have been specified')
                 raise
 
         if self.do_rotate:
@@ -86,23 +86,10 @@ class Plate(UserObjectGeometry):
         xf, yf, zf = p2
 
         # Check for valid orientations
-        if xs == xf:
-            if ys == yf or zs == zf:
-                logger.exception(self.__str__() + ' the plate is not specified correctly')
-                raise ValueError
-
-        elif ys == yf:
-            if xs == xf or zs == zf:
-                logger.exception(self.__str__() + ' the plate is not specified correctly')
-                raise ValueError
-
-        elif zs == zf:
-            if xs == xf or ys == yf:
-                logger.exception(self.__str__() + ' the plate is not specified correctly')
-                raise ValueError
-
-        else:
-            logger.exception(self.__str__() + ' the plate is not specified correctly')
+        if ((xs == xf and (ys == yf or zs == zf))
+            or (ys == yf and (xs == xf or zs == zf))
+            or (zs == zf and (ys != yf and xs != xf))):
+            logger.exception(f'{self.__str__()} the plate is not specified correctly')
             raise ValueError
 
         # Look up requested materials in existing list of material instances
@@ -110,7 +97,7 @@ class Plate(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         # yz-plane plate
@@ -161,6 +148,6 @@ class Plate(UserObjectGeometry):
                     build_face_xy(i, j, zs, numIDx, numIDy, grid.rigidE, 
                                   grid.rigidH, grid.ID)
 
-        logger.info(self.grid_name(grid) + f"Plate from {p3[0]:g}m, {p3[1]:g}m, " +
+        logger.info(f"{self.grid_name(grid)}Plate from {p3[0]:g}m, {p3[1]:g}m, " +
                     f"{p3[2]:g}m, to {p4[0]:g}m, {p4[1]:g}m, {p4[2]:g}m of " +
                     f"material(s) {', '.join(materialsrequested)} created.")
diff --git a/gprMax/cmds_geometry/sphere.py b/gprMax/cmds_geometry/sphere.py
index 014ec3b0..31e7c465 100644
--- a/gprMax/cmds_geometry/sphere.py
+++ b/gprMax/cmds_geometry/sphere.py
@@ -48,7 +48,7 @@ class Sphere(UserObjectGeometry):
             p1 = self.kwargs['p1']
             r = self.kwargs['r']
         except KeyError:
-            logger.exception(self.__str__() + ' please specify a point and a radius.')
+            logger.exception(f'{self.__str__()} please specify a point and a radius.')
             raise
 
         # Check averaging
@@ -68,7 +68,7 @@ class Sphere(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' no materials have been specified')
+                logger.exception(f'{self.__str__()} no materials have been specified')
                 raise
 
         # Centre of sphere
@@ -85,7 +85,7 @@ class Sphere(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         # Isotropic case
@@ -125,7 +125,7 @@ class Sphere(UserObjectGeometry):
                      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, " +
+        logger.info(f"{self.grid_name(grid)}Sphere with centre {p2[0]:g}m, " +
                     f"{p2[1]:g}m, {p2[2]:g}m, radius {r:g}m, of material(s) " +
                     f"{', '.join(materialsrequested)} created, dielectric " +
                     f"smoothing is {dielectricsmoothing}.")
diff --git a/gprMax/cmds_geometry/triangle.py b/gprMax/cmds_geometry/triangle.py
index 17ab6077..77a334db 100644
--- a/gprMax/cmds_geometry/triangle.py
+++ b/gprMax/cmds_geometry/triangle.py
@@ -70,7 +70,7 @@ class Triangle(UserObjectGeometry):
             up3 = self.kwargs['p3']
             thickness = self.kwargs['thickness']
         except KeyError:
-            logger.exception(self.__str__() + ' specify 3 points and a thickness')
+            logger.exception(f'{self.__str__()} specify 3 points and a thickness')
             raise
         
         if self.do_rotate:
@@ -93,7 +93,7 @@ class Triangle(UserObjectGeometry):
             try:
                 materialsrequested = self.kwargs['material_ids']
             except KeyError:
-                logger.exception(self.__str__() + ' no materials have been specified')
+                logger.exception(f'{self.__str__()} no materials have been specified')
                 raise
 
         p4 = uip.round_to_grid_static_point(up1)
@@ -108,21 +108,21 @@ class Triangle(UserObjectGeometry):
         x3, y3, z3 = uip.round_to_grid(up3)
         
         if thickness < 0:
-            logger.exception(self.__str__() + ' requires a positive value for thickness')
+            logger.exception(f'{self.__str__()} requires a positive value for thickness')
             raise ValueError
 
         # Check for valid orientations
         # yz-plane triangle
-        if x1 == x2 and x2 == x3:
+        if x1 == x2 == x3:
             normal = 'x'
         # xz-plane triangle
-        elif y1 == y2 and y2 == y3:
+        elif y1 == y2 == y3:
             normal = 'y'
         # xy-plane triangle
-        elif z1 == z2 and z2 == z3:
+        elif z1 == z2 == z3:
             normal = 'z'
         else:
-            logger.exception(self.__str__() + ' the triangle is not specified correctly')
+            logger.exception(f'{self.__str__()} the triangle is not specified correctly')
             raise ValueError
 
         # Look up requested materials in existing list of material instances
@@ -130,7 +130,7 @@ class Triangle(UserObjectGeometry):
 
         if len(materials) != len(materialsrequested):
             notfound = [x for x in materialsrequested if x not in materials]
-            logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+            logger.exception(f'{self.__str__()} material(s) {notfound} do not exist')
             raise ValueError
 
         if thickness > 0:
@@ -185,13 +185,13 @@ class Triangle(UserObjectGeometry):
 
         if thickness > 0:
             dielectricsmoothing = 'on' if averaging else 'off'
-            logger.info(self.grid_name(grid) + f"Triangle with coordinates " +
+            logger.info(f"{self.grid_name(grid)}Triangle with coordinates " +
                         f"{p4[0]:g}m {p4[1]:g}m {p4[2]:g}m, {p5[0]:g}m {p5[1]:g}m " +
                         f"{p5[2]:g}m, {p6[0]:g}m {p6[1]:g}m {p6[2]:g}m and thickness " +
                         f"{thickness:g}m of material(s) {', '.join(materialsrequested)} " +
                         f"created, dielectric smoothing is {dielectricsmoothing}.")
         else:
-            logger.info(self.grid_name(grid) + f"Triangle with coordinates " +
+            logger.info(f"{self.grid_name(grid)}Triangle with coordinates " +
                         f"{p4[0]:g}m {p4[1]:g}m {p4[2]:g}m, {p5[0]:g}m {p5[1]:g}m " +
                         f"{p5[2]:g}m, {p6[0]:g}m {p6[1]:g}m {p6[2]:g}m of material(s) " +
                         f"{', '.join(materialsrequested)} created.")
diff --git a/gprMax/cmds_multiuse.py b/gprMax/cmds_multiuse.py
index 31be50ee..b9f26adc 100644
--- a/gprMax/cmds_multiuse.py
+++ b/gprMax/cmds_multiuse.py
@@ -61,9 +61,9 @@ class UserObjectMulti:
         """Readable user string as per hash commands."""
         s = ''
         for _, v in self.kwargs.items():
-            if isinstance(v, tuple) or isinstance(v, list):
+            if isinstance(v, (tuple, list)):
                 v = ' '.join([str(el) for el in v])
-            s += str(v) + ' '
+            s += f'{str(v)} '
 
         return f'{self.hash}: {s[:-1]}'
 
@@ -77,7 +77,7 @@ class UserObjectMulti:
 
     def params_str(self):
         """Readable string of parameters given to object."""
-        return self.hash + ': ' + str(self.kwargs)
+        return f'{self.hash}: {str(self.kwargs)}'
 
     def grid_name(self, grid):
         """Returns subgrid name for use with logging info. Returns an empty
@@ -114,12 +114,12 @@ class Waveform(UserObjectMulti):
         try:
             wavetype = self.kwargs['wave_type'].lower()
         except KeyError:
-            logger.exception(self.params_str() + (f" must have one of the " 
-                             f"following types {','.join(WaveformUser.types)}."))
+            logger.exception(f"{self.params_str()} must have one of the " + 
+                             f"following types {','.join(WaveformUser.types)}.")
             raise
         if wavetype not in WaveformUser.types:
-            logger.exception(self.params_str() + (f" must have one of the "
-                             f"following types {','.join(WaveformUser.types)}."))
+            logger.exception(f"{self.params_str()} must have one of the " +
+                             f"following types {','.join(WaveformUser.types)}.")
             raise ValueError
 
         if wavetype != 'user':
@@ -243,18 +243,15 @@ class VoltageSource(UserObjectMulti):
             logger.exception(self.params_str() + (' polarisation must be '
                                                   'x, y, or z.'))
             raise ValueError
-        if ('2D TMx' in config.get_model_config().mode and 
-                        (polarisation == 'y' or polarisation == 'z')):
+        if '2D TMx' in config.get_model_config().mode and polarisation in ['y','z',]:
             logger.exception(self.params_str() + (' polarisation must be x in '
                                                   '2D TMx mode.'))
             raise ValueError
-        elif ('2D TMy' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'z')):
+        elif '2D TMy' in config.get_model_config().mode and polarisation in ['x','z',]:
             logger.exception(self.params_str() + (' polarisation must be y in '
                                                   '2D TMy mode.'))
             raise ValueError
-        elif ('2D TMz' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'y')):
+        elif '2D TMz' in config.get_model_config().mode and polarisation in ['x','y',]:
             logger.exception(self.params_str() + (' polarisation must be z in '
                                                   '2D TMz mode.'))
             raise ValueError
@@ -305,10 +302,7 @@ class VoltageSource(UserObjectMulti):
                                                       'less.'))
                 raise ValueError
             v.start = start
-            if stop > grid.timewindow:
-                v.stop = grid.timewindow
-            else:
-                v.stop = stop
+            v.stop = min(stop, grid.timewindow)
             startstop = (f' start time {v.start:g} secs, finish time '
                          f'{v.stop:g} secs ')
         except KeyError:
@@ -318,7 +312,7 @@ class VoltageSource(UserObjectMulti):
 
         v.calculate_waveform_values(grid)
 
-        logger.info(self.grid_name(grid) + f'Voltage source with polarity '
+        logger.info(f'{self.grid_name(grid)}Voltage source with polarity '
                     f'{v.polarisation} at {p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m, '
                     f'resistance {v.resistance:.1f} Ohms,' + startstop + 
                     f'using waveform {v.waveformID} created.')
@@ -367,8 +361,7 @@ class HertzianDipole(UserObjectMulti):
             p1 = self.kwargs['p1']
             waveform_id = self.kwargs['waveform_id']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at least 3 ' + 
-                             'parameters.')
+            logger.exception(f'{self.params_str()} requires at least 3 parameters.')
             raise
 
         if self.do_rotate:
@@ -379,18 +372,15 @@ class HertzianDipole(UserObjectMulti):
             logger.exception(self.params_str() + (' polarisation must be '
                                                   'x, y, or z.'))
             raise ValueError
-        if ('2D TMx' in config.get_model_config().mode and 
-                        (polarisation == 'y' or polarisation == 'z')):
+        if '2D TMx' in config.get_model_config().mode and polarisation in ['y','z',]:
             logger.exception(self.params_str() + (' polarisation must be x in '
                                                   '2D TMx mode.'))
             raise ValueError
-        elif ('2D TMy' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'z')):
+        elif '2D TMy' in config.get_model_config().mode and polarisation in ['x','z',]:
             logger.exception(self.params_str() + (' polarisation must be y in '
                                                   '2D TMy mode.'))
             raise ValueError
-        elif ('2D TMz' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'y')):
+        elif '2D TMz' in config.get_model_config().mode and polarisation in ['x','y',]:
             logger.exception(self.params_str() + (' polarisation must be z in '
                                                   '2D TMz mode.'))
             raise ValueError
@@ -401,7 +391,7 @@ class HertzianDipole(UserObjectMulti):
 
         # Check if there is a waveformID in the waveforms list
         if not any(x.ID == waveform_id for x in grid.waveforms):
-            logger.exception(self.params_str() + ' there is no waveform ' + 
+            logger.exception(f'{self.params_str()} there is no waveform ' + 
                              f'with the identifier {waveform_id}.')
             raise ValueError
 
@@ -422,8 +412,7 @@ class HertzianDipole(UserObjectMulti):
         h.xcoordorigin = xcoord
         h.ycoordorigin = ycoord
         h.zcoordorigin = zcoord
-        h.ID = (h.__class__.__name__ + '(' + str(h.xcoord) + ',' + 
-                str(h.ycoord) + ',' + str(h.zcoord) + ')')
+        h.ID = f'{h.__class__.__name__}({str(h.xcoord)},{str(h.ycoord)},{str(h.zcoord)})'
         h.waveformID = waveform_id
 
         try:
@@ -431,25 +420,19 @@ class HertzianDipole(UserObjectMulti):
             start = self.kwargs['start']
             stop = self.kwargs['stop']
             if start < 0:
-                logger.exception(self.params_str() + (' delay of the initiation '
-                                                      'of the source should not '
-                                                      'be less than zero.'))
+                logger.exception(f'{self.params_str()} delay of the initiation of the ' +
+                                 f'source should not be less than zero.')
                 raise ValueError
             if stop < 0:
-                logger.exception(self.params_str() + (' time to remove the '
-                                                      'source should not be '
-                                                      'less than zero.'))
+                logger.exception(f'{self.params_str()} time to remove the source ' +
+                                 f'should not be less than zero.')
                 raise ValueError
             if stop - start <= 0:
-                logger.exception(self.params_str() + (' duration of the source '
-                                                      'should not be zero or '
-                                                      'less.'))
+                logger.exception(
+                    f'{self.params_str()} duration of the source should not be zero or less.')
                 raise ValueError
             h.start = start
-            if stop > grid.timewindow:
-                h.stop = grid.timewindow
-            else:
-                h.stop = stop
+            h.stop = min(stop, grid.timewindow)
             startstop = (f' start time {h.start:g} secs, finish time '
                          f'{h.stop:g} secs ')
         except KeyError:
@@ -460,12 +443,12 @@ class HertzianDipole(UserObjectMulti):
         h.calculate_waveform_values(grid)
 
         if config.get_model_config().mode == '2D':
-            logger.info(self.grid_name(grid) + f'Hertzian dipole is a line ' +
+            logger.info(f'{self.grid_name(grid)}Hertzian dipole is a line ' +
                         f'source in 2D with polarity {h.polarisation} at ' + 
                         f'{p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m,' + startstop + 
                         f'using waveform {h.waveformID} created.')
         else:
-            logger.info(self.grid_name(grid) + f'Hertzian dipole with ' + 
+            logger.info(f'{self.grid_name(grid)}Hertzian dipole with ' + 
                         f'polarity {h.polarisation} at {p2[0]:g}m, ' + 
                         f'{p2[1]:g}m, {p2[2]:g}m,' + startstop + 
                         f'using waveform {h.waveformID} created.')
@@ -514,8 +497,7 @@ class MagneticDipole(UserObjectMulti):
             p1 = self.kwargs['p1']
             waveform_id = self.kwargs['waveform_id']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at least five ' 
-                             'parameters.')
+            logger.exception(f'{self.params_str()} requires at least five parameters.')
             raise
 
         if self.do_rotate:
@@ -526,18 +508,15 @@ class MagneticDipole(UserObjectMulti):
             logger.exception(self.params_str() + (' polarisation must be '
                                                   'x, y, or z.'))
             raise ValueError
-        if ('2D TMx' in config.get_model_config().mode and 
-                        (polarisation == 'y' or polarisation == 'z')):
+        if '2D TMx' in config.get_model_config().mode and polarisation in ['y','z',]:
             logger.exception(self.params_str() + (' polarisation must be x in '
                                                   '2D TMx mode.'))
             raise ValueError
-        elif ('2D TMy' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'z')):
+        elif '2D TMy' in config.get_model_config().mode and polarisation in ['x','z',]:
             logger.exception(self.params_str() + (' polarisation must be y in '
                                                   '2D TMy mode.'))
             raise ValueError
-        elif ('2D TMz' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'y')):
+        elif '2D TMz' in config.get_model_config().mode and polarisation in ['x','y',]:
             logger.exception(self.params_str() + (' polarisation must be z in '
                                                   '2D TMz mode.'))
             raise ValueError
@@ -548,7 +527,7 @@ class MagneticDipole(UserObjectMulti):
 
         # Check if there is a waveformID in the waveforms list
         if not any(x.ID == waveform_id for x in grid.waveforms):
-            logger.exception(self.params_str() + ' there is no waveform ' + 
+            logger.exception(f'{self.params_str()} there is no waveform ' + 
                              f'with the identifier {waveform_id}.')
             raise ValueError
 
@@ -584,10 +563,7 @@ class MagneticDipole(UserObjectMulti):
                                                       'less.'))
                 raise ValueError
             m.start = start
-            if stop > grid.timewindow:
-                m.stop = grid.timewindow
-            else:
-                m.stop = stop
+            m.stop = min(stop, grid.timewindow)
             startstop = (f' start time {m.start:g} secs, '
                          f'finish time {m.stop:g} secs ')
         except KeyError:
@@ -597,7 +573,7 @@ class MagneticDipole(UserObjectMulti):
 
         m.calculate_waveform_values(grid)
 
-        logger.info(self.grid_name(grid) + f'Magnetic dipole with polarity ' +
+        logger.info(f'{self.grid_name(grid)}Magnetic dipole with polarity ' +
                     f'{m.polarisation} at {p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m,' + 
                     startstop + f'using waveform {m.waveformID} created.')
 
@@ -647,17 +623,15 @@ class TransmissionLine(UserObjectMulti):
             waveform_id = self.kwargs['waveform_id']
             resistance = self.kwargs['resistance']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at least six '
-                             'parameters.')
+            logger.exception(f'{self.params_str()} requires at least six parameters.')
             raise
 
         if self.do_rotate:
             self._do_rotate(grid)
 
         # Warn about using a transmission line on GPU
-        if (config.sim_config.general['solver'] == 'cuda' or 
-            config.sim_config.general['solver'] == 'opencl'):
-            logger.exception(self.params_str() + ' cannot currently be used ' +
+        if (config.sim_config.general['solver'] in ['cuda', 'opencl']):
+            logger.exception(f'{self.params_str()} cannot currently be used ' +
                              'with the CUDA or OpenCL-based solver. Consider ' +
                              'using a #voltage_source instead.')
             raise ValueError
@@ -667,18 +641,15 @@ class TransmissionLine(UserObjectMulti):
             logger.exception(self.params_str() + (' polarisation must be '
                                                   'x, y, or z.'))
             raise ValueError
-        if ('2D TMx' in config.get_model_config().mode and 
-                        (polarisation == 'y' or polarisation == 'z')):
+        if '2D TMx' in config.get_model_config().mode and polarisation in ['y','z',]:
             logger.exception(self.params_str() + (' polarisation must be x in '
                                                   '2D TMx mode.'))
             raise ValueError
-        elif ('2D TMy' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'z')):
+        elif '2D TMy' in config.get_model_config().mode and polarisation in ['x','z',]:
             logger.exception(self.params_str() + (' polarisation must be y in '
                                                   '2D TMy mode.'))
             raise ValueError
-        elif ('2D TMz' in config.get_model_config().mode and 
-                          (polarisation == 'x' or polarisation == 'y')):
+        elif '2D TMz' in config.get_model_config().mode and polarisation in ['x','y',]:
             logger.exception(self.params_str() + (' polarisation must be z in '
                                                   '2D TMz mode.'))
             raise ValueError
@@ -688,14 +659,14 @@ class TransmissionLine(UserObjectMulti):
 
 
         if resistance <= 0 or resistance >= config.sim_config.em_consts['z0']:
-            logger.exception(self.params_str() + ' requires a resistance ' + 
+            logger.exception(f'{self.params_str()} requires a resistance ' + 
                              'greater than zero and less than the impedance ' +
                              'of free space, i.e. 376.73 Ohms.')
             raise ValueError
 
         # Check if there is a waveformID in the waveforms list
         if not any(x.ID == waveform_id for x in grid.waveforms):
-            logger.exception(self.params_str() + ' there is no waveform ' + 
+            logger.exception(f'{self.params_str()} there is no waveform ' + 
                              f'with the identifier {waveform_id}.')
             raise ValueError
 
@@ -729,11 +700,8 @@ class TransmissionLine(UserObjectMulti):
                                                       'less.'))
                 raise ValueError
             t.start = start
-            if stop > grid.timewindow:
-                t.stop = grid.timewindow
-            else:
-                t.stop = stop
-            startstop = (f' start time {t.start:g} secs, finish time '
+            t.stop = min(stop, grid.timewindow)
+            startstop = (f' start time {t.start:g} secs, finish time ' +
                          f'{t.stop:g} secs ')
         except KeyError:
             t.start = 0
@@ -743,7 +711,7 @@ class TransmissionLine(UserObjectMulti):
         t.calculate_waveform_values(grid)
         t.calculate_incident_V_I(grid)
 
-        logger.info(self.grid_name(grid) + f'Transmission line with polarity ' +
+        logger.info(f'{self.grid_name(grid)}Transmission line with polarity ' +
                     f'{t.polarisation} at {p2[0]:g}m, {p2[1]:g}m, ' + 
                     f'{p2[2]:g}m, resistance {t.resistance:.1f} Ohms,' + 
                     startstop + f'using waveform {t.waveformID} created.')
@@ -817,16 +785,14 @@ class Rx(UserObjectMulti):
             outputs = [self.kwargs['outputs']]
         except KeyError:
             # If no ID or outputs are specified, use default
-            r.ID = (r.__class__.__name__ + '(' + str(r.xcoord) + ',' + 
-                    str(r.ycoord) + ',' + str(r.zcoord) + ')')
+            r.ID = f'{r.__class__.__name__}({str(r.xcoord)},{str(r.ycoord)},{str(r.zcoord)})'
             for key in RxUser.defaultoutputs:
                 r.outputs[key] = np.zeros(grid.iterations, 
                                           dtype=config.sim_config.dtypes['float_or_double'])
         else:
             outputs.sort()
             # Get allowable outputs
-            if (config.sim_config.general['solver'] =='cuda' or 
-                config.sim_config.general['solver'] =='opencl'):
+            if config.sim_config.general['solver'] in ['cuda', 'opencl']:
                 allowableoutputs = RxUser.allowableoutputs_dev
             else:
                 allowableoutputs = RxUser.allowableoutputs
@@ -836,14 +802,14 @@ class Rx(UserObjectMulti):
                     r.outputs[field] = np.zeros(grid.iterations, 
                                                 dtype=config.sim_config.dtypes['float_or_double'])
                 else:
-                    logger.exception(self.params_str() + ' contains an '
-                                     'output type that is not allowable. '
-                                     'Allowable outputs in current context are '
+                    logger.exception(f'{self.params_str()} contains an output '
+                                     f'type that is not allowable. Allowable '
+                                     f'outputs in current context are '
                                      f'{allowableoutputs}.')
                     raise ValueError
 
-        logger.info(self.grid_name(grid) + f"Receiver at {p2[0]:g}m, " 
-                    f"{p2[1]:g}m, {p2[2]:g}m with output component(s) "
+        logger.info(f"{self.grid_name(grid)}Receiver at {p2[0]:g}m, {p2[1]:g}m, "
+                    f"{p2[2]:g}m with output component(s) "
                     f"{', '.join(r.outputs)} created.")
 
         grid.rxs.append(r)
@@ -871,8 +837,7 @@ class RxArray(UserObjectMulti):
             p2 = self.kwargs['p2']
             dl = self.kwargs['dl']
         except KeyError:
-            logger.exception(self.params_str() + ' requires exactly 9 '
-                             'parameters')
+            logger.exception(f'{self.params_str()} requires exactly 9 parameters')
             raise
 
         xs, ys, zs = uip.check_src_rx_point(p1, self.params_str(), 'lower')
@@ -882,39 +847,39 @@ class RxArray(UserObjectMulti):
         dx, dy, dz = uip.discretise_point(dl)
 
         if xs > xf or ys > yf or zs > zf:
-            logger.exception(self.params_str() + ' the lower coordinates ' 
+            logger.exception(f'{self.params_str()} the lower coordinates ' + 
                              'should be less than the upper coordinates.')
             raise ValueError
         if dx < 0 or dy < 0 or dz < 0:
-            logger.exception(self.params_str() + ' the step size should not ' 
+            logger.exception(f'{self.params_str()} the step size should not ' + 
                              'be less than zero.')
             raise ValueError
         if dx < 1:
             if dx == 0:
                 dx = 1
             else:
-                logger.exception(self.params_str() + ' the step size should ' 
+                logger.exception(f'{self.params_str()} the step size should ' +
                                  'not be less than the spatial discretisation.')
                 raise ValueError
         if dy < 1:
             if dy == 0:
                 dy = 1
             else:
-                logger.exception(self.params_str() + ' the step size should ' 
+                logger.exception(f'{self.params_str()} the step size should ' +
                                  'not be less than the spatial discretisation.')
                 raise ValueError
         if dz < 1:
             if dz == 0:
                 dz = 1
             else:
-                logger.exception(self.params_str() + ' the step size should '
+                logger.exception(f'{self.params_str()} the step size should ' +
                                  'not be less than the spatial discretisation.')
                 raise ValueError
 
-        logger.info(self.grid_name(grid) + f'Receiver array {p3[0]:g}m, ' 
-                    f'{p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, {p4[1]:g}m, '
-                    f'{p4[2]:g}m with steps {dx * grid.dx:g}m, '
-                    f'{dy * grid.dy:g}m, {dz * grid.dz:g}m')
+        logger.info(f'{self.grid_name(grid)}Receiver array '
+                    f'{p3[0]:g}m, {p3[1]:g}m, {p3[2]:g}m, to '
+                    f'{p4[0]:g}m, {p4[1]:g}m, {p4[2]:g}m with steps '
+                    f'{dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m')
 
         for x in range(xs, xf + 1, dx):
             for y in range(ys, yf + 1, dy):
@@ -930,8 +895,7 @@ class RxArray(UserObjectMulti):
                     p5 = np.array([x, y, z])
                     p5 = uip.descretised_to_continuous(p5)
                     p5 = uip.round_to_grid_static_point(p5)
-                    r.ID = (r.__class__.__name__ + '(' + str(x) + ',' + 
-                            str(y) + ',' + str(z) + ')')
+                    r.ID = f'{r.__class__.__name__}({str(x)},{str(y)},{str(z)})'
                     for key in RxUser.defaultoutputs:
                         r.outputs[key] = np.zeros(grid.iterations, dtype=config.sim_config.dtypes['float_or_double'])
                     logger.info(f"  Receiver at {p5[0]:g}m, {p5[1]:g}m, "
@@ -968,8 +932,7 @@ class Snapshot(UserObjectMulti):
 
     def create(self, grid, uip):
         if isinstance(grid, SubGridBaseGrid):
-            logger.exception(self.params_str() + ' do not add snapshots to '
-                             'subgrids.')
+            logger.exception(f'{self.params_str()} do not add snapshots to subgrids.')
             raise ValueError
         try:
             p1 = self.kwargs['p1']
@@ -977,8 +940,7 @@ class Snapshot(UserObjectMulti):
             dl = self.kwargs['dl']
             filename = self.kwargs['filename']
         except KeyError:
-            logger.exception(self.params_str() + ' requires exactly 11 '
-                             'parameters.')
+            logger.exception(f'{self.params_str()} requires exactly 11 parameters.')
             raise
 
         try:
@@ -986,7 +948,7 @@ class Snapshot(UserObjectMulti):
             p4 = uip.round_to_grid_static_point(p2)
             p1, p2 = uip.check_box_points(p1, p2, self.params_str())
         except ValueError:
-            logger.exception(self.params_str() + ' point is outside the domain.')
+            logger.exception(f'{self.params_str()} point is outside the domain.')
             raise
         xs, ys, zs = p1
         xf, yf, zf = p2
@@ -1001,13 +963,12 @@ class Snapshot(UserObjectMulti):
             try:
                 time = self.kwargs['time']
             except KeyError:
-                logger.exception(self.params_str() + ' requires exactly 5 '
-                                 'parameters.')
+                logger.exception(f'{self.params_str()} requires exactly 5 parameters.')
                 raise
             if time > 0:
                 iterations = round_value((time / grid.dt)) + 1
             else:
-                logger.exception(self.params_str() + ' time value must be '
+                logger.exception(f'{self.params_str()} time value must be ' + 
                                  'greater than zero.')
                 raise ValueError
 
@@ -1027,10 +988,10 @@ class Snapshot(UserObjectMulti):
             # Check and set output names
             for output in tmp:
                 if output not in SnapshotUser.allowableoutputs.keys():
-                    logger.exception(self.params_str() + " contains an "
-                                        "output type that is not allowable. "
-                                        "Allowable outputs in current context are "
-                                        f"{', '.join(SnapshotUser.allowableoutputs.keys())}.")
+                    logger.exception(f"{self.params_str()} contains an output "
+                                     f"type that is not allowable. Allowable "
+                                     f"outputs in current context are "
+                                     f"{', '.join(SnapshotUser.allowableoutputs.keys())}.")
                     raise ValueError
                 else:
                     outputs[output] = True
@@ -1039,15 +1000,15 @@ class Snapshot(UserObjectMulti):
             outputs = dict.fromkeys(SnapshotUser.allowableoutputs, True)
 
         if dx < 0 or dy < 0 or dz < 0:
-            logger.exception(self.params_str() + ' the step size should not '
+            logger.exception(f'{self.params_str()} the step size should not ' +
                              'be less than zero.')
             raise ValueError
         if dx < 1 or dy < 1 or dz < 1:
-            logger.exception(self.params_str() + ' the step size should not '
+            logger.exception(f'{self.params_str()} the step size should not ' +
                              'be less than the spatial discretisation.')
             raise ValueError
         if iterations <= 0 or iterations > grid.iterations:
-            logger.exception(self.params_str() + ' time value is not valid.')
+            logger.exception(f'{self.params_str()} time value is not valid.')
             raise ValueError
 
         s = SnapshotUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, iterations, 
@@ -1088,39 +1049,37 @@ class Material(UserObjectMulti):
             sm = self.kwargs['sm']
             material_id = self.kwargs['id']
         except KeyError:
-            logger.exception(self.params_str() + ' requires exactly five '
-                             'parameters.')
+            logger.exception(f'{self.params_str()} requires exactly five parameters.')
             raise
 
         if er < 1:
-            logger.exception(self.params_str() + ' requires a positive value '
-                             'of one or greater for static (DC) permittivity.')
+            logger.exception(f'{self.params_str()} requires a positive value of ' +
+                             f'one or greater for static (DC) permittivity.')
             raise ValueError
         if se != 'inf':
             se = float(se)
             if se < 0:
-                logger.exception(self.params_str() + ' requires a positive '
-                                 'value for electric conductivity.')
+                logger.exception(f'{self.params_str()} requires a positive ' +
+                                 f'value for electric conductivity.')
                 raise ValueError
         else:
             se = float('inf')
         if mr < 1:
-            logger.exception(self.params_str() + ' requires a positive value '
-                             'of one or greater for magnetic permeability.')
+            logger.exception(f'{self.params_str()} requires a positive value of ' +
+                             f'one or greater for magnetic permeability.')
             raise ValueError
         if sm < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
-                             'for magnetic loss.')
+            logger.exception(f'{self.params_str()} requires a positive value ' +
+                             f'for magnetic loss.')
             raise ValueError
         if any(x.ID == material_id for x in grid.materials):
-            logger.exception(self.params_str() + f' with ID {material_id} '
-                             'already exists')
+            logger.exception(f'{self.params_str()} with ID {material_id} ' +
+                             f'already exists')
             raise ValueError
 
         # Create a new instance of the Material class material 
         # (start index after pec & free_space)
         m = MaterialUser(len(grid.materials), material_id)
-        m.er = er
         m.se = se
         m.mr = mr
         m.sm = sm
@@ -1129,9 +1088,10 @@ class Material(UserObjectMulti):
         if m.se == float('inf'):
             m.averagable = False
 
-        logger.info(self.grid_name(grid) + f'Material {m.ID} with '
-                    f'eps_r={m.er:g}, sigma={m.se:g} S/m; mu_r={m.mr:g}, '
-                    f'sigma*={m.sm:g} Ohm/m created.')
+        m.er = er
+        logger.info(f'{self.grid_name(grid)}Material {m.ID} with eps_r={m.er:g}, '
+                    f'sigma={m.se:g} S/m; mu_r={m.mr:g}, sigma*={m.sm:g} Ohm/m '
+                    f'created.')
 
         grid.materials.append(m)
 
@@ -1162,12 +1122,11 @@ class AddDebyeDispersion(UserObjectMulti):
             tau = self.kwargs['tau']
             material_ids = self.kwargs['material_ids']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at least four '
-                             'parameters.')
+            logger.exception(f'{self.params_str()} requires at least four parameters.')
             raise
 
         if poles < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for number of poles.')
             raise ValueError
 
@@ -1176,8 +1135,7 @@ class AddDebyeDispersion(UserObjectMulti):
 
         if len(materials) != len(material_ids):
             notfound = [x for x in material_ids if x not in materials]
-            logger.exception(self.params_str() + f' material(s) {notfound} do '
-                             'not exist')
+            logger.exception(f'{self.params_str()} material(s) {notfound} do not exist')
             raise ValueError
 
         for material in materials:
@@ -1189,14 +1147,14 @@ class AddDebyeDispersion(UserObjectMulti):
             disp_material.type = 'debye'
             disp_material.poles = poles
             disp_material.averagable = False
-            for i in range(0, poles):
+            for i in range(poles):
                 if tau[i] > 0:
                     logger.debug('Not checking if relaxation times are '
                                  'greater than time-step.')
                     disp_material.deltaer.append(er_delta[i])
                     disp_material.tau.append(tau[i])
                 else:
-                    logger.exception(self.params_str() + ' requires positive '
+                    logger.exception(f'{self.params_str()} requires positive ' +
                                      'values for the permittivity difference.')
                     raise ValueError
             if disp_material.poles > config.get_model_config().materials['maxpoles']:
@@ -1205,8 +1163,8 @@ class AddDebyeDispersion(UserObjectMulti):
             # Replace original material with newly created DispersiveMaterial
             grid.materials = [disp_material if mat.numID==material.numID else mat for mat in grid.materials]
 
-            logger.info(self.grid_name(grid) + f"Debye disperion added to {disp_material.ID} "
-                        f"with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in disp_material.deltaer)}, " 
+            logger.info(f"{self.grid_name(grid)}Debye disperion added to {disp_material.ID} "
+                        f"with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in disp_material.deltaer)}, "
                         f"and tau={', '.join('%4.3e' % tau for tau in disp_material.tau)} secs created.")
 
 
@@ -1238,12 +1196,11 @@ class AddLorentzDispersion(UserObjectMulti):
             delta = self.kwargs['delta']
             material_ids = self.kwargs['material_ids']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at least five '
-                             'parameters.')
+            logger.exception(f'{self.params_str()} requires at least five parameters.')
             raise
 
         if poles < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for number of poles.')
             raise ValueError
 
@@ -1252,8 +1209,7 @@ class AddLorentzDispersion(UserObjectMulti):
 
         if len(materials) != len(material_ids):
             notfound = [x for x in material_ids if x not in materials]
-            logger.exception(self.params_str() + f' material(s) {notfound} do '
-                             'not exist')
+            logger.exception(f'{self.params_str()} material(s) {notfound} do not exist')
             raise ValueError
 
         for material in materials:
@@ -1265,13 +1221,13 @@ class AddLorentzDispersion(UserObjectMulti):
             disp_material.type = 'lorentz'
             disp_material.poles = poles
             disp_material.averagable = False
-            for i in range(0, poles):
+            for i in range(poles):
                 if er_delta[i] > 0 and omega[i] > grid.dt and delta[i] > grid.dt:
                     disp_material.deltaer.append(er_delta[i])
                     disp_material.tau.append(omega[i])
                     disp_material.alpha.append(delta[i])
                 else:
-                    logger.exception(self.params_str() + ' requires positive '
+                    logger.exception(f'{self.params_str()} requires positive '
                                      'values for the permittivity difference '
                                      'and frequencies, and associated times '
                                      'that are greater than the time step for '
@@ -1283,10 +1239,11 @@ class AddLorentzDispersion(UserObjectMulti):
             # Replace original material with newly created DispersiveMaterial
             grid.materials = [disp_material if mat.numID==material.numID else mat for mat in grid.materials]
 
-            logger.info(self.grid_name(grid) + f"Lorentz disperion added to {disp_material.ID} "
-                        f"with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in disp_material.deltaer)}, "
-                        f"omega={', '.join('%4.3e' % omega for omega in disp_material.tau)} secs, "
-                        f"and gamma={', '.join('%4.3e' % delta for delta in disp_material.alpha)} created.")
+            logger.info(
+                f"{self.grid_name(grid)}Lorentz disperion added to {disp_material.ID} " +
+                f"with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in disp_material.deltaer)}, " +
+                f"omega={', '.join('%4.3e' % omega for omega in disp_material.tau)} secs, " +
+                f"and gamma={', '.join('%4.3e' % delta for delta in disp_material.alpha)} created.")
 
 
 class AddDrudeDispersion(UserObjectMulti):
@@ -1313,12 +1270,12 @@ class AddDrudeDispersion(UserObjectMulti):
             alpha = self.kwargs['alpha']
             material_ids = self.kwargs['material_ids']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at least four '
+            logger.exception(f'{self.params_str()} requires at least four ' +
                              'parameters.')
             raise
 
         if poles < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for number of poles.')
             raise ValueError
 
@@ -1327,8 +1284,7 @@ class AddDrudeDispersion(UserObjectMulti):
 
         if len(materials) != len(material_ids):
             notfound = [x for x in material_ids if x not in materials]
-            logger.exception(self.params_str() + f' material(s) {notfound} do '
-                             'not exist.')
+            logger.exception(f'{self.params_str()} material(s) {notfound} do not exist.')
             raise ValueError
 
         for material in materials:
@@ -1340,14 +1296,14 @@ class AddDrudeDispersion(UserObjectMulti):
             disp_material.type = 'drude'
             disp_material.poles = poles
             disp_material.averagable = False
-            for i in range(0, poles):
+            for i in range(poles):
                 if omega[i] > 0 and alpha[i] > grid.dt:
                     disp_material.tau.append(omega[i])
                     disp_material.alpha.append(alpha[i])
                 else:
-                    logger.exception(self.params_str() + ' requires positive '
-                                     'values for the frequencies, and '
-                                     'associated times that are greater than '
+                    logger.exception(f'{self.params_str()} requires positive ' +
+                                     'values for the frequencies, and ' +
+                                     'associated times that are greater than ' +
                                      'the time step for the model.')
                     raise ValueError
             if disp_material.poles > config.get_model_config().materials['maxpoles']:
@@ -1356,9 +1312,10 @@ class AddDrudeDispersion(UserObjectMulti):
             # Replace original material with newly created DispersiveMaterial
             grid.materials = [disp_material if mat.numID==material.numID else mat for mat in grid.materials]
 
-            logger.info(self.grid_name(grid) + f"Drude disperion added to {disp_material.ID} "
-                        f"with omega={', '.join('%4.3e' % omega for omega in disp_material.tau)} secs, "
-                        f"and gamma={', '.join('%4.3e' % alpha for alpha in disp_material.alpha)} secs created.")
+            logger.info(
+                f"{self.grid_name(grid)}Drude disperion added to {disp_material.ID} "
+                f"with omega={', '.join('%4.3e' % omega for omega in disp_material.tau)} secs, "
+                f"and gamma={', '.join('%4.3e' % alpha for alpha in disp_material.alpha)} secs created.")
 
 
 class SoilPeplinski(UserObjectMulti):
@@ -1392,38 +1349,38 @@ class SoilPeplinski(UserObjectMulti):
             water_fraction_upper = self.kwargs['water_fraction_upper']
             ID = self.kwargs['id']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at exactly seven '
+            logger.exception(f'{self.params_str()} requires at exactly seven ' +
                              'parameters.')
             raise
 
         if sand_fraction < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the sand fraction.')
             raise ValueError
         if clay_fraction < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the clay fraction.')
             raise ValueError
         if bulk_density < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the bulk density.')
             raise ValueError
         if sand_density < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the sand particle density.')
             raise ValueError
         if water_fraction_lower < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the lower limit of the water volumetric '
                              'fraction.')
             raise ValueError
         if water_fraction_upper < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the upper limit of the water volumetric '
                              'fraction.')
             raise ValueError
         if any(x.ID == ID for x in grid.mixingmodels):
-            logger.exception(self.params_str() + f' with ID {ID} already exists')
+            logger.exception(f'{self.params_str()} with ID {ID} already exists')
             raise ValueError
 
         # Create a new instance of the Material class material 
@@ -1431,10 +1388,10 @@ class SoilPeplinski(UserObjectMulti):
         s = PeplinskiSoilUser(ID, sand_fraction, clay_fraction, bulk_density, 
                               sand_density, (water_fraction_lower, water_fraction_upper))
 
-        logger.info(self.grid_name(grid) + 'Mixing model (Peplinski) used to '
-                    f'create {s.ID} with sand fraction {s.S:g}, clay fraction '
-                    f'{s.C:g}, bulk density {s.rb:g}g/cm3, sand particle '
-                    f'density {s.rs:g}g/cm3, and water volumetric fraction '
+        logger.info(f'{self.grid_name(grid)}Mixing model (Peplinski) used to ' +
+                    f'create {s.ID} with sand fraction {s.S:g}, clay fraction ' +
+                    f'{s.C:g}, bulk density {s.rb:g}g/cm3, sand particle ' +
+                    f'density {s.rs:g}g/cm3, and water volumetric fraction ' +
                     f'{s.mu[0]:g} to {s.mu[1]:g} created.')
 
         grid.mixingmodels.append(s)
@@ -1474,58 +1431,59 @@ class MaterialRange(UserObjectMulti):
             ro_upper = self.kwargs['ro_upper']
             ID = self.kwargs['id']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at exactly nine '
+            logger.exception(f'{self.params_str()} requires at exactly nine ' +
                              'parameters.')
             raise
 
         if er_lower < 1:
-            logger.exception(self.params_str() + ' requires a value greater or equal to 1 '
+            logger.exception(f'{self.params_str()} requires a value greater or equal to 1 ' +
                              'for the lower range of relative permittivity.')
             raise ValueError
         if mr_lower < 1:
-            logger.exception(self.params_str() + ' requires a value greater or equal to 1 '
+            logger.exception(f'{self.params_str()} requires a value greater or equal to 1 ' +
                              'for the lower range of relative magnetic permeability.')
             raise ValueError
         if sigma_lower < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the lower limit of conductivity.')
             raise ValueError
         if ro_lower < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the lower range magnetic loss.')
             raise ValueError
         if er_upper < 1:
-            logger.exception(self.params_str() + ' requires a value greater or equal to 1'
+            logger.exception(f'{self.params_str()} requires a value greater or equal to 1' +
                             'for the upper range of relative permittivity.')
             raise ValueError
-        
+
         if mr_upper < 1:
-            logger.exception(self.params_str() + ' requires a value greater or equal to 1'
+            logger.exception(f'{self.params_str()} requires a value greater or equal to 1' +
                              'for the upper range of relative magnetic permeability')
             raise ValueError
-        
+
         if sigma_upper < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the upper range of conductivity.')
             raise ValueError
         if ro_upper < 0:
-            logger.exception(self.params_str() + ' requires a positive value '
+            logger.exception(f'{self.params_str()} requires a positive value ' +
                              'for the upper range of magnetic loss.')
 
-        
+
         if any(x.ID == ID for x in grid.mixingmodels):
-            logger.exception(self.params_str() + f' with ID {ID} already exists')
+            logger.exception(f'{self.params_str()} with ID {ID} already exists')
             raise ValueError
 
         # Create a new instance of the Material class material 
         # (start index after pec & free_space)
-        s = RangeMaterialUser(ID, (er_lower, er_upper), (sigma_lower, sigma_upper), (mr_lower, mr_upper), (ro_lower, ro_upper))
+        s = RangeMaterialUser(ID, (er_lower, er_upper), (sigma_lower, sigma_upper), 
+                              (mr_lower, mr_upper), (ro_lower, ro_upper))
 
-        logger.info(self.grid_name(grid) + 'Material properties used to '
-                    f'create {s.ID} with range(s) {s.er[0]:g} to {s.er[1]:g}, relative permittivity '
-                    f'{s.sig[0]:g} to {s.sig[1]:g}, S/m conductivity, {s.mu[0]:g} to {s.mu[1]:g} relative magnetic permeability '
+        logger.info(f'{self.grid_name(grid)}Material properties used to ' +
+                    f'create {s.ID} with range(s) {s.er[0]:g} to {s.er[1]:g}, relative permittivity ' +
+                    f'{s.sig[0]:g} to {s.sig[1]:g}, S/m conductivity, {s.mu[0]:g} to {s.mu[1]:g} relative magnetic permeability ' +
                     f'{s.ro[0]:g} to {s.ro[1]:g} Ohm/m magnetic loss, created')
-        
+
 
         grid.mixingmodels.append(s)
 
@@ -1550,27 +1508,25 @@ class MaterialList(UserObjectMulti):
             list_of_materials = self.kwargs['list_of_materials']
             ID = self.kwargs['id']
         except KeyError:
-            logger.exception(self.params_str() + ' requires at at least 2 '
+            logger.exception(f'{self.params_str()} requires at at least 2 ' +
                              'parameters.')
             raise
 
-               
         if any(x.ID == ID for x in grid.mixingmodels):
-            logger.exception(self.params_str() + f' with ID {ID} already exists')
+            logger.exception(f'{self.params_str()} with ID {ID} already exists')
             raise ValueError
 
         # Create a new instance of the Material class material 
         # (start index after pec & free_space)
         s = ListMaterialUser(ID, list_of_materials)
 
-        logger.info(self.grid_name(grid) + 'A list of materials used to '
+        logger.info(f'{self.grid_name(grid)}A list of materials used to ' +
                     f'create {s.ID} that includes {s.mat}, created')
-        
+
 
         grid.mixingmodels.append(s)
 
 
-
 class GeometryView(UserObjectMulti):
     """Outputs to file(s) information about the geometry (mesh) of model.
 
@@ -1613,7 +1569,7 @@ class GeometryView(UserObjectMulti):
             output_type = self.kwargs['output_type'].lower()
             filename = self.kwargs['filename']
         except KeyError:
-            logger.exception(self.params_str() + ' requires exactly eleven '
+            logger.exception(f'{self.params_str()} requires exactly eleven ' +
                              'parameters.')
             raise
 
@@ -1624,45 +1580,44 @@ class GeometryView(UserObjectMulti):
             p4 = uip.round_to_grid_static_point(p2)
             p1, p2 = uip.check_box_points(p1, p2, self.params_str())
         except ValueError:
-            logger.exception(self.params_str() + ' point is outside the domain.')
+            logger.exception(f'{self.params_str()} point is outside the domain.')
             raise
         xs, ys, zs = p1
         xf, yf, zf = p2
 
         dx, dy, dz = uip.discretise_static_point(dl)
 
-
         if dx < 0 or dy < 0 or dz < 0:
-            logger.exception(self.params_str() + ' the step size should not be '
+            logger.exception(f'{self.params_str()} the step size should not be ' +
                              'less than zero.')
             raise ValueError
         if dx > grid.nx or dy > grid.ny or dz > grid.nz:
-            logger.exception(self.params_str() + ' the step size should be '
+            logger.exception(f'{self.params_str()} the step size should be ' +
                              'less than the domain size.')
             raise ValueError
         if dx < 1 or dy < 1 or dz < 1:
-            logger.exception(self.params_str() + ' the step size should not '
+            logger.exception(f'{self.params_str()} the step size should not ' +
                              'be less than the spatial discretisation.')
             raise ValueError
-        if output_type != 'n' and output_type != 'f':
-            logger.exception(self.params_str() + ' requires type to be either '
+        if output_type not in ['n', 'f']:
+            logger.exception(f'{self.params_str()} requires type to be either ' +
                              'n (normal) or f (fine).')
             raise ValueError
         if output_type == 'f' and (dx * grid.dx != grid.dx or 
                                    dy * grid.dy != grid.dy or 
                                    dz * grid.dz != grid.dz):
-            logger.exception(self.params_str() + ' requires the spatial '
-                             'discretisation for the geometry view to be the '
-                             'same as the model for geometry view of '
+            logger.exception(f'{self.params_str()} requires the spatial ' +
+                             'discretisation for the geometry view to be the ' +
+                             'same as the model for geometry view of ' +
                              'type f (fine)')
             raise ValueError
 
         g = GeometryViewUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, grid)
 
-        logger.info(self.grid_name(grid) + f'Geometry view from {p3[0]:g}m, '
-                    f'{p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, {p4[1]:g}m, '
-                    f'{p4[2]:g}m, discretisation {dx * grid.dx:g}m, '
-                    f'{dy * grid.dy:g}m, {dz * grid.dz:g}m, with filename '
+        logger.info(f'{self.grid_name(grid)}Geometry view from {p3[0]:g}m, ' +
+                    f'{p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, {p4[1]:g}m, ' +
+                    f'{p4[2]:g}m, discretisation {dx * grid.dx:g}m, ' +
+                    f'{dy * grid.dy:g}m, {dz * grid.dz:g}m, with filename ' +
                     f'base {g.filename} created.')
 
         grid.geometryviews.append(g)
@@ -1692,7 +1647,7 @@ class GeometryObjectsWrite(UserObjectMulti):
             p2 = self.kwargs['p2']
             basefilename = self.kwargs['filename']
         except KeyError:
-            logger.exception(self.params_str() + ' requires exactly seven '
+            logger.exception(f'{self.params_str()} requires exactly seven ' +
                              'parameters.')
             raise
 
@@ -1756,30 +1711,30 @@ class PMLCFS(UserObjectMulti):
             sigmamin = self.kwargs['sigmamin']
             sigmamax = self.kwargs['sigmamax']
         except KeyError:
-            logger.exception(self.params_str() + ' requires exactly twelve '
+            logger.exception(f'{self.params_str()} requires exactly twelve ' +
                              'parameters.')
             raise
 
         if (alphascalingprofile not in CFSParameter.scalingprofiles.keys() or 
             kappascalingprofile not in CFSParameter.scalingprofiles.keys() or 
             sigmascalingprofile not in CFSParameter.scalingprofiles.keys()):
-            logger.exception(self.params_str() + ' must have scaling type '
+            logger.exception(f'{self.params_str()} must have scaling type ' +
                              f"{','.join(CFSParameter.scalingprofiles.keys())}")
             raise ValueError
         if (alphascalingdirection not in CFSParameter.scalingdirections or 
             kappascalingdirection not in CFSParameter.scalingdirections or 
             sigmascalingdirection not in CFSParameter.scalingdirections):
-            logger.exception(self.params_str() + ' must have scaling type '
+            logger.exception(f'{self.params_str()} must have scaling type ' +
                              f"{','.join(CFSParameter.scalingdirections)}")
             raise ValueError
         if (float(alphamin) < 0 or float(alphamax) < 0 or 
             float(kappamin) < 0 or float(kappamax) < 0 or float(sigmamin) < 0):
-            logger.exception(self.params_str() + ' minimum and maximum scaling '
+            logger.exception(f'{self.params_str()} minimum and maximum scaling ' +
                              'values must be greater than zero.')
             raise ValueError
         # TODO: Fix handling of kappa for 2nd order PMLs
         # if float(kappamin) < 1:
-        #     logger.exception(self.params_str() + ' minimum scaling value for '
+        #     logger.exception(f'{self.params_str()} minimum scaling value for '
         #                      'kappa must be greater than or equal to one.')
         #     raise ValueError
 
@@ -1819,7 +1774,7 @@ class PMLCFS(UserObjectMulti):
         grid.pmls['cfs'].append(cfs)
 
         if len(grid.pmls['cfs']) > 2:
-            logger.exception(self.params_str() + ' can only be used up to two '
+            logger.exception(f'{self.params_str()} can only be used up to two ' +
                              'times, for up to a 2nd order PML.')
             raise ValueError
             
diff --git a/gprMax/cmds_singleuse.py b/gprMax/cmds_singleuse.py
index 2a38beb6..c0de561d 100644
--- a/gprMax/cmds_singleuse.py
+++ b/gprMax/cmds_singleuse.py
@@ -93,20 +93,20 @@ class Discretisation(UserObjectSingle):
             G.dl = np.array(self.kwargs['p1'])
             G.dx, G.dy, G.dz = self.kwargs['p1']
         except KeyError:
-            logger.exception(self.__str__() + ' discretisation requires a point')
+            logger.exception(f'{self.__str__()} discretisation requires a point')
             raise
 
         if G.dl[0] <= 0:
-            logger.exception(self.__str__() + ' discretisation requires the ' +
-                             'x-direction spatial step to be greater than zero')
+            logger.exception(f'{self.__str__()} discretisation requires the '
+                             f'x-direction spatial step to be greater than zero')
             raise ValueError
         if G.dl[1] <= 0:
-            logger.exception(self.__str__() + ' discretisation requires the ' +
-                             'y-direction spatial step to be greater than zero')
+            logger.exception(f'{self.__str__()} discretisation requires the '
+                             f'y-direction spatial step to be greater than zero')
             raise ValueError
         if G.dl[2] <= 0:
-            logger.exception(self.__str__() + ' discretisation requires the ' +
-                             'z-direction spatial step to be greater than zero')
+            logger.exception(f'{self.__str__()} discretisation requires the '
+                             f'z-direction spatial step to be greater than zero')
             raise ValueError
 
         logger.info(f'Spatial discretisation: {G.dl[0]:g} x {G.dl[1]:g} x {G.dl[2]:g}m')
@@ -127,11 +127,12 @@ class Domain(UserObjectSingle):
         try:
             G.nx, G.ny, G.nz = uip.discretise_point(self.kwargs['p1'])
         except KeyError:
-            logger.exception(self.__str__() + ' please specify a point')
+            logger.exception(f'{self.__str__()} please specify a point')
             raise
 
         if G.nx == 0 or G.ny == 0 or G.nz == 0:
-            logger.exception(self.__str__() + ' requires at least one cell in every dimension')
+            logger.exception(f'{self.__str__()} requires at least one cell in '
+                             f'every dimension')
             raise ValueError
 
         logger.info(f"Domain size: {self.kwargs['p1'][0]:g} x {self.kwargs['p1'][1]:g} x " +
@@ -181,12 +182,12 @@ class TimeStepStabilityFactor(UserObjectSingle):
         try:
             f = self.kwargs['f']
         except KeyError:
-            logger.exception(self.__str__() + ' requires exactly one parameter')
+            logger.exception(f'{self.__str__()} requires exactly one parameter')
             raise
 
         if f <= 0 or f > 1:
-            logger.exception(self.__str__() + ' requires the value of the time ' +
-                             'step stability factor to be between zero and one')
+            logger.exception(f'{self.__str__()} requires the value of the time '
+                             f'step stability factor to be between zero and one')
             raise ValueError
         G.dt = G.dt * f
 
@@ -250,12 +251,12 @@ class OMPThreads(UserObjectSingle):
         try:
             n = self.kwargs['n']
         except KeyError:
-            logger.exception(self.__str__() + ' requires exactly one parameter ' +
-                             'to specify the number of CPU OpenMP threads to use')
+            logger.exception(f'{self.__str__()} requires exactly one parameter '
+                             f'to specify the number of CPU OpenMP threads to use')
             raise
         if n < 1:
-            logger.exception(self.__str__() + ' requires the value to be an ' +
-                             'integer not less than one')
+            logger.exception(f'{self.__str__()} requires the value to be an '
+                             f'integer not less than one')
             raise ValueError
 
         config.get_model_config().ompthreads = set_omp_threads(n)
@@ -302,7 +303,7 @@ class PMLProps(UserObjectSingle):
                 G.pmls['thickness']['ymax'] = int(self.kwargs['ymax'])
                 G.pmls['thickness']['zmax'] = int(self.kwargs['zmax'])
             except KeyError:
-                logger.exception(self.__str__() + ' requires either one or six parameter(s)')
+                logger.exception(f'{self.__str__()} requires either one or six parameter(s)')
                 raise
 
         if (2 * G.pmls['thickness']['x0'] >= G.nx or
@@ -311,8 +312,8 @@ class PMLProps(UserObjectSingle):
             2 * G.pmls['thickness']['xmax'] >= G.nx or
             2 * G.pmls['thickness']['ymax'] >= G.ny or
             2 * G.pmls['thickness']['zmax'] >= G.nz):
-                logger.exception(self.__str__() + ' has too many cells for the domain size')
-                raise ValueError
+            logger.exception(f'{self.__str__()} has too many cells for the domain size')
+            raise ValueError
 
 
 class SrcSteps(UserObjectSingle):
@@ -330,12 +331,12 @@ class SrcSteps(UserObjectSingle):
         try:
             G.srcsteps = uip.discretise_point(self.kwargs['p1'])
         except KeyError:
-            logger.exception(self.__str__() + ' requires exactly three parameters')
+            logger.exception(f'{self.__str__()} requires exactly three parameters')
             raise
 
-        logger.info(f'Simple sources will step {G.srcsteps[0] * G.dx:g}m, ' +
-                    f'{G.srcsteps[1] * G.dy:g}m, {G.srcsteps[2] * G.dz:g}m ' +
-                    f'for each model run.')
+        logger.info(f'Simple sources will step {G.srcsteps[0] * G.dx:g}m, ' + 
+                    f'{G.srcsteps[1] * G.dy:g}m, {G.srcsteps[2] * G.dz:g}m ' + 
+                    'for each model run.')
 
 
 class RxSteps(UserObjectSingle):
@@ -353,12 +354,12 @@ class RxSteps(UserObjectSingle):
         try:
             G.rxsteps = uip.discretise_point(self.kwargs['p1'])
         except KeyError:
-            logger.exception(self.__str__() + ' requires exactly three parameters')
+            logger.exception(f'{self.__str__()} requires exactly three parameters')
             raise
 
-        logger.info(f'All receivers will step {G.rxsteps[0] * G.dx:g}m, ' +
-                    f'{G.rxsteps[1] * G.dy:g}m, {G.rxsteps[2] * G.dz:g}m ' +
-                    f'for each model run.')
+        logger.info(f'All receivers will step {G.rxsteps[0] * G.dx:g}m, ' + 
+                    f'{G.rxsteps[1] * G.dy:g}m, {G.rxsteps[2] * G.dz:g}m ' + 
+                    'for each model run.')
 
 
 class ExcitationFile(UserObjectSingle):
@@ -378,7 +379,7 @@ class ExcitationFile(UserObjectSingle):
 
     def create(self, G, uip):
         try:
-            kwargs = dict()
+            kwargs = {}
             excitationfile = self.kwargs['filepath']
             kwargs['kind'] = self.kwargs['kind']
             kwargs['fill_value'] = self.kwargs['fill_value']
@@ -389,7 +390,7 @@ class ExcitationFile(UserObjectSingle):
                 args, varargs, keywords, defaults = inspect.getargspec(interpolate.interp1d)
                 kwargs = dict(zip(reversed(args), reversed(defaults)))
             except KeyError:
-                logger.exception(self.__str__() + ' requires either one or three parameter(s)')
+                logger.exception(f'{self.__str__()} requires either one or three parameter(s)')
                 raise
 
         # See if file exists at specified path and if not try input file directory
diff --git a/gprMax/config.py b/gprMax/config.py
index c0544d24..162b06c8 100644
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -224,14 +224,14 @@ class SimulationConfig:
         #   solver: cpu, cuda, opencl.
         #   subgrid: whether the simulation uses sub-grids.
         #   precision: data type for electromagnetic field output (single/double).
+        #   progressbars: progress bars on stdoout or not - switch off 
+        #                   progressbars when logging level is greater than 
+        #                   info (20) 
 
         self.general = {'solver': 'cpu',
                         'subgrid': False,
-                        'precision': 'single'}
-
-        # Progress bars on stdoout or not - switch off progressbars when
-        # logging level is greater than info (20)
-        self.general['progressbars'] = False if args.log_level > 20 else True
+                        'precision': 'single',
+                        'progressbars': args.log_level <= 20}
 
         self.em_consts = {'c': c, # Speed of light in free space (m/s)
                           'e0': e0, # Permittivity of free space (F/m)
@@ -289,7 +289,6 @@ class SimulationConfig:
 
         # Set more complex parameters
         self._set_precision()
-        self._get_byteorder()
         self._set_input_file_path()
         self._set_model_start_end()
 
@@ -348,12 +347,6 @@ class SimulationConfig:
             elif self.general['solver'] == 'opencl':
                 self.dtypes['C_complex'] = 'cdouble'
 
-    def _get_byteorder(self):
-        """Checks the byte order of system to use for VTK files, i.e. geometry
-            views and snapshots.
-        """
-        self.vtk_byteorder = 'LittleEndian' if sys.byteorder == 'little' else 'BigEndian'
-
     def _set_model_start_end(self):
         """Sets range for number of models to run (internally 0 index)."""
         if self.args.i:
diff --git a/gprMax/contexts.py b/gprMax/contexts.py
index f0e3c952..fbde3698 100644
--- a/gprMax/contexts.py
+++ b/gprMax/contexts.py
@@ -92,7 +92,7 @@ class Context:
 
     def print_logo_copyright(self):
         """Prints gprMax logo, version, and copyright/licencing information."""
-        logo_copyright = logo(__version__ + ' (' + codename + ')')
+        logo_copyright = logo(f'{__version__} ({codename})')
         logger.basic(logo_copyright)
 
     def print_sim_time_taken(self):
@@ -164,19 +164,15 @@ class MPIContext(Context):
         executor = self.MPIExecutor(self._run_model, comm=self.comm)
 
         # Check GPU resources versus number of MPI tasks
-        if executor.is_master():
-            if config.sim_config.general['solver'] == 'cuda':
-                if executor.size - 1 > len(config.sim_config.devices['devs']):
-                    logger.exception('Not enough GPU resources for number of '
-                                     'MPI tasks requested. Number of MPI tasks '
-                                     'should be equal to number of GPUs + 1.')
-                    raise ValueError
-
-        # Create job list
-        jobs = []
-        for i in self.model_range:
-            jobs.append({'i': i})
+        if (executor.is_master() and
+            config.sim_config.general['solver'] == 'cuda' and
+            executor.size - 1 > len(config.sim_config.devices['devs'])):
+            logger.exception('Not enough GPU resources for number of '
+                             'MPI tasks requested. Number of MPI tasks '
+                             'should be equal to number of GPUs + 1.')
+            raise ValueError
 
+        jobs = [{'i': i} for i in self.model_range]
         # Send the workers to their work loop
         executor.start()
         if executor.is_master():
diff --git a/gprMax/cython/pml_build.pyx b/gprMax/cython/pml_build.pyx
new file mode 100644
index 00000000..421ca1a0
--- /dev/null
+++ b/gprMax/cython/pml_build.pyx
@@ -0,0 +1,70 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
+
+cpdef pml_average_er_mr(
+    str dir,
+    int s,
+    G
+):
+    """Calculates average permittivity and permeability for building PML 
+        (based on underlying material er and mr from solid array).
+
+    Args:
+        dir: string identifier for direction of PML.
+        s: int for starting cell of PML.
+        G: FDTDGrid class describing a grid in a model.
+    """
+
+    sumer = 0  # Sum of relative permittivities in PML slab
+    summr = 0  # Sum of relative permeabilities in PML slab
+
+    if dir == 'x':
+        for j in range(G.ny):
+            for k in range(G.nz):
+                numID = G.solid[s, j, k]
+                material = [x for x in G.materials if x.numID == numID]
+                material = material[0]
+                sumer += material.er
+                summr += material.mr
+        averageer = sumer / (G.ny * G.nz)
+        averagemr = summr / (G.ny * G.nz)
+
+    elif dir == 'y':
+        for i in range(G.nx):
+            for k in range(G.nz):
+                numID = G.solid[i, s, k]
+                material = [x for x in G.materials if x.numID == numID]
+                material = material[0]
+                sumer += material.er
+                summr += material.mr
+        averageer = sumer / (G.nx * G.nz)
+        averagemr = summr / (G.nx * G.nz)
+
+    elif dir == 'z':
+        for i in range(G.nx):
+            for j in range(G.ny):
+                numID = G.solid[i, j, s]
+                material = [x for x in G.materials if x.numID == numID]
+                material = material[0]
+                sumer += material.er
+                summr += material.mr
+        averageer = sumer / (G.nx * G.ny)
+        averagemr = summr / (G.nx * G.ny)
+
+    return averageer, averagemr
diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py
index 7ddc9a5d..d8050c9a 100644
--- a/gprMax/fields_outputs.py
+++ b/gprMax/fields_outputs.py
@@ -76,7 +76,7 @@ def write_hdf5_outputfile(outputfile, G):
     # Write meta data and data for any subgrids
     if sg_rxs:
         for sg in G.subgrids:
-            grp = f.create_group('/subgrids/' + sg.name)
+            grp = f.create_group(f'/subgrids/{sg.name}')
             write_hd5_data(grp, sg, is_subgrid=True)
 
     if G.rxs or sg_rxs:
@@ -115,7 +115,7 @@ def write_hd5_data(basegrp, G, is_subgrid=False):
     # Create group for sources (except transmission lines); add type and positional data attributes
     srclist = G.voltagesources + G.hertziandipoles + G.magneticdipoles
     for srcindex, src in enumerate(srclist):
-        grp = basegrp.create_group('srcs/src' + str(srcindex + 1))
+        grp = basegrp.create_group(f'srcs/src{str(srcindex + 1)}')
         grp.attrs['Type'] = type(src).__name__
         grp.attrs['Position'] = (src.xcoord * G.dx, src.ycoord * G.dy, src.zcoord * G.dz)
 
diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py
index f93d30a9..dbe16a84 100644
--- a/gprMax/geometry_outputs.py
+++ b/gprMax/geometry_outputs.py
@@ -277,14 +277,10 @@ class Comments():
         """
 
         # Comments for Paraview macro
-        comments = {}
-
-        comments['gprMax_version'] = __version__
-        comments['dx_dy_dz'] = self.dx_dy_dz_comment()
-        comments['nx_ny_nz'] = self.nx_ny_nz_comment()
-
-        # Write the name and numeric ID for each material
-        comments['Materials'] = self.materials_comment()
+        comments = {'gprMax_version': __version__, 
+                    'dx_dy_dz': self.dx_dy_dz_comment(),
+                    'nx_ny_nz': self.nx_ny_nz_comment(),
+                    'Materials': self.materials_comment()} # Write the name and numeric ID for each material
 
         # Write information on PMLs, sources, and receivers
         if not self.materials_only:
@@ -345,7 +341,7 @@ class Comments():
 
     def materials_comment(self):
         if not self.averaged_materials:
-            return [m.ID for m in self.grid.materials if '+' not in m.ID]
+            return [m.ID for m in self.grid.materials if m.type is not 'dielectric-smoothed']
         else:
             return [m.ID for m in self.grid.materials]
 
@@ -375,8 +371,7 @@ class GeometryObjects:
         parts = config.sim_config.input_file_path.with_suffix('').parts
         self.filename_hdf5 = Path(*parts[:-1], self.basefilename)
         self.filename_hdf5 = self.filename_hdf5.with_suffix('.h5')
-        self.filename_materials = Path(
-            *parts[:-1], self.basefilename + '_materials')
+        self.filename_materials = Path(*parts[:-1], f'{self.basefilename}_materials')
         self.filename_materials = self.filename_materials.with_suffix('.txt')
 
         # Sizes of arrays to write necessary to update progress bar
diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py
index 805835e5..845b21d5 100644
--- a/gprMax/hash_cmds_file.py
+++ b/gprMax/hash_cmds_file.py
@@ -178,14 +178,14 @@ def write_processed_file(processedlines):
         for item in processedlines:
             f.write(f'{item}')
 
-    logger.info(f'Written input commands, after processing any Python code and ' +
-                f'include commands, to file: {processedfile}\n')
+    logger.info(f'Written input commands, after processing any Python ' +
+                f'code and include commands, to file: {processedfile}\n')
 
 
 def check_cmd_names(processedlines, checkessential=True):
     """Checks the validity of commands, i.e. are they gprMax commands,
         and that all essential commands are present.
-
+        
     Args:
         processedlines: list of input commands after Python processing.
         checkessential: boolean to check for essential commands or not.
@@ -212,6 +212,7 @@ def check_cmd_names(processedlines, checkessential=True):
     # - these will be lists within the dictionary
     multiplecmds = {key: [] for key in ['#geometry_view', 
                                         '#geometry_objects_write', '#material',
+                                        '#material_range', '#material_list',
                                         '#soil_peplinski', 
                                         '#add_dispersion_debye', 
                                         '#add_dispersion_lorentz', 
@@ -219,14 +220,14 @@ def check_cmd_names(processedlines, checkessential=True):
                                         '#waveform', '#voltage_source', 
                                         '#hertzian_dipole', '#magnetic_dipole', 
                                         '#transmission_line', '#rx', '#rx_array', 
-                                        '#snapshot', '#include_file', '#material_range', '#material_list']}
+                                        '#snapshot', '#include_file']}
 
     # Geometry object building commands that there can be multiple instances
     # of in a model - these will be lists within the dictionary
     geometrycmds = ['#geometry_objects_read', '#edge', '#plate', '#triangle', 
-                    '#box', '#sphere', '#ellipsoid', '#cone', '#cylinder', '#cylindrical_sector', 
-                    '#fractal_box', '#add_surface_roughness', 
-                    '#add_surface_water', '#add_grass']
+                    '#box', '#sphere', '#ellipsoid', '#cone', '#cylinder', 
+                    '#cylindrical_sector', '#fractal_box', 
+                    '#add_surface_roughness', '#add_surface_water', '#add_grass']
     # List to store all geometry object commands in order from input file
     geometry = []
 
diff --git a/gprMax/hash_cmds_geometry.py b/gprMax/hash_cmds_geometry.py
index dbe8f9c8..b8c7da89 100644
--- a/gprMax/hash_cmds_geometry.py
+++ b/gprMax/hash_cmds_geometry.py
@@ -17,21 +17,22 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import logging
+
 import numpy as np
 
 from .cmds_geometry.add_grass import AddGrass
 from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness
 from .cmds_geometry.add_surface_water import AddSurfaceWater
 from .cmds_geometry.box import Box
-from .cmds_geometry.cylinder import Cylinder
 from .cmds_geometry.cone import Cone
+from .cmds_geometry.cylinder import Cylinder
 from .cmds_geometry.cylindrical_sector import CylindricalSector
 from .cmds_geometry.edge import Edge
+from .cmds_geometry.ellipsoid import Ellipsoid
 from .cmds_geometry.fractal_box import FractalBox
 from .cmds_geometry.plate import Plate
 from .cmds_geometry.sphere import Sphere
 from .cmds_geometry.triangle import Triangle
-from .cmds_geometry.ellipsoid import Ellipsoid
 from .utilities.utilities import round_value
 
 logger = logging.getLogger(__name__)
@@ -76,7 +77,8 @@ def process_geometrycmds(geometry):
         tmp = object.split()
 
         if tmp[0] == '#geometry_objects_read:':
-            from .cmds_geometry.geometry_objects_read import GeometryObjectsRead
+            from .cmds_geometry.geometry_objects_read import \
+                GeometryObjectsRead
             
             if len(tmp) != 6:
                 logger.exception("'" + ' '.join(tmp) + "'" + 
@@ -219,7 +221,6 @@ def process_geometrycmds(geometry):
 
             scene_objects.append(cylinder)
 
-
         elif tmp[0] == '#cone:':
             if len(tmp) < 10:
                 logger.exception("'" + ' '.join(tmp) + "'" + 
@@ -239,7 +240,7 @@ def process_geometrycmds(geometry):
             elif len(tmp) == 11:
                 averaging = check_averaging(tmp[10].lower())
                 cone = Cone(p1=p1, p2=p2, r1=r1, r2=r2, material_id=tmp[9], 
-                                    averaging=averaging)
+                            averaging=averaging)
 
             # Uniaxial anisotropic case
             elif len(tmp) == 12:
@@ -252,7 +253,6 @@ def process_geometrycmds(geometry):
 
             scene_objects.append(cone)
 
-
         elif tmp[0] == '#cylindrical_sector:':
             if len(tmp) < 10:
                 logger.exception("'" + ' '.join(tmp) + "'" + 
@@ -331,7 +331,6 @@ def process_geometrycmds(geometry):
 
             scene_objects.append(sphere)
 
-
         elif tmp[0] == '#ellipsoid:':
             if len(tmp) < 8:
                 logger.exception("'" + ' '.join(tmp) + "'" + 
@@ -345,17 +344,19 @@ def process_geometrycmds(geometry):
 
             # Isotropic case with no user specified averaging
             if len(tmp) == 8:
-                ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, material_id=tmp[7])
+                ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, 
+                                      material_id=tmp[7])
 
             # Isotropic case with user specified averaging
             elif len(tmp) == 9:
                 averaging = check_averaging(tmp[8].lower())
-                ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, material_id=tmp[7], 
-                                averaging=averaging)
+                ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, 
+                                      material_id=tmp[7], averaging=averaging)
 
             # Uniaxial anisotropic case
             elif len(tmp) == 8:
-                ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, material_id=tmp[7:])
+                ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, 
+                                      material_id=tmp[7:])
 
             else:
                 logger.exception("'" + ' '.join(tmp) + "'" + 
@@ -364,8 +365,6 @@ def process_geometrycmds(geometry):
 
             scene_objects.append(ellipsoid)
 
-    
-
         elif tmp[0] == '#fractal_box:':
             # Default is no dielectric smoothing for a fractal box
 
@@ -479,4 +478,4 @@ def process_geometrycmds(geometry):
 
                     scene_objects.append(grass)
 
-    return scene_objects
+    return scene_objects
\ No newline at end of file
diff --git a/gprMax/hash_cmds_multiuse.py b/gprMax/hash_cmds_multiuse.py
index a46431a4..6925ae8c 100644
--- a/gprMax/hash_cmds_multiuse.py
+++ b/gprMax/hash_cmds_multiuse.py
@@ -21,8 +21,9 @@ import logging
 from .cmds_multiuse import (AddDebyeDispersion, AddDrudeDispersion,
                             AddLorentzDispersion, GeometryObjectsWrite,
                             GeometryView, HertzianDipole, MagneticDipole,
-                            Material, Rx, RxArray, Snapshot, SoilPeplinski,
-                            TransmissionLine, VoltageSource, Waveform, MaterialRange, MaterialList)
+                            Material, MaterialList, MaterialRange, Rx, RxArray,
+                            Snapshot, SoilPeplinski, TransmissionLine,
+                            VoltageSource, Waveform)
 
 logger = logging.getLogger(__name__)
 
@@ -352,17 +353,16 @@ def process_multicmds(multicmds):
                                  ' requires at exactly nine parameters')
                 raise ValueError
             material_range = MaterialRange(er_lower=float(tmp[0]),
-                                 er_upper=float(tmp[1]),
-                                 sigma_lower=float(tmp[2]),
-                                 sigma_upper=float(tmp[3]),
-                                 mr_lower=float(tmp[4]),
-                                 mr_upper=float(tmp[5]),
-                                 ro_lower=float(tmp[6]),
-                                 ro_upper=float(tmp[7]),
-                                 id=tmp[8])
+                                           er_upper=float(tmp[1]),
+                                           sigma_lower=float(tmp[2]),
+                                           sigma_upper=float(tmp[3]),
+                                           mr_lower=float(tmp[4]),
+                                           mr_upper=float(tmp[5]),
+                                           ro_lower=float(tmp[6]),
+                                           ro_upper=float(tmp[7]),
+                                           id=tmp[8])
             scene_objects.append(material_range)
 
-
     cmdname = '#material_list'
     if multicmds[cmdname] is not None:
         for cmdinstance in multicmds[cmdname]:
@@ -375,14 +375,11 @@ def process_multicmds(multicmds):
             
             tokens = len(tmp)
             lmats = []
-            for iter in range(0,tokens-1):
+            for iter in range(tokens-1):
                 lmats.append(tmp[iter])
 
-
-
-            material_list = MaterialList(list_of_materials=lmats,
-                               id=tmp[tokens-1])
+            material_list = MaterialList(list_of_materials=lmats, 
+                                         id=tmp[tokens-1])
             scene_objects.append(material_list)
 
-
-    return scene_objects
+    return scene_objects
\ No newline at end of file
diff --git a/gprMax/hash_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py
index 4ca0c223..86ac37cc 100644
--- a/gprMax/hash_cmds_singleuse.py
+++ b/gprMax/hash_cmds_singleuse.py
@@ -55,8 +55,8 @@ def process_singlecmds(singlecmds):
     if singlecmds[cmd] is not None:
         tmp = tuple(int(x) for x in singlecmds[cmd].split())
         if len(tmp) != 1:
-            logger.exception(cmd + ' requires exactly one parameter to specify ' +
-                             'the number of CPU OpenMP threads to use')
+            logger.exception(f'{cmd} requires exactly one parameter to specify ' +
+                             f'the number of CPU OpenMP threads to use')
             raise ValueError
 
         omp_threads = OMPThreads(n=tmp[0])
@@ -66,7 +66,7 @@ def process_singlecmds(singlecmds):
     if singlecmds[cmd] is not None:
         tmp = [float(x) for x in singlecmds[cmd].split()]
         if len(tmp) != 3:
-            logger.exception(cmd + ' requires exactly three parameters')
+            logger.exception(f'{cmd} requires exactly three parameters')
             raise ValueError
 
         dl = (tmp[0], tmp[1], tmp[2])
@@ -77,7 +77,7 @@ def process_singlecmds(singlecmds):
     if singlecmds[cmd] is not None:
         tmp = [float(x) for x in singlecmds[cmd].split()]
         if len(tmp) != 3:
-            logger.exception(cmd + ' requires exactly three parameters')
+            logger.exception(f'{cmd} requires exactly three parameters')
             raise ValueError
 
         p1 = (tmp[0], tmp[1], tmp[2])
@@ -94,8 +94,8 @@ def process_singlecmds(singlecmds):
     if singlecmds[cmd] is not None:
         tmp = singlecmds[cmd].split()
         if len(tmp) != 1:
-            logger.exception(cmd + ' requires exactly one parameter to specify ' +
-                             'the time window. Either in seconds or number of iterations.')
+            logger.exception(f'{cmd} requires exactly one parameter to specify the ' +
+                             f'time window. Either in seconds or number of iterations.')
             raise ValueError
         tmp = tmp[0].lower()
 
@@ -113,8 +113,8 @@ def process_singlecmds(singlecmds):
     cmd = '#pml_cells'
     if singlecmds[cmd] is not None:
         tmp = singlecmds[cmd].split()
-        if len(tmp) != 1 and len(tmp) != 6:
-            logger.exception(cmd + ' requires either one or six parameter(s)')
+        if len(tmp) not in [1, 6]:
+            logger.exception(f'{cmd} requires either one or six parameter(s)')
             raise ValueError
         if len(tmp) == 1:
             pml_cells = PMLProps(thickness=int(tmp[0]))
@@ -132,7 +132,7 @@ def process_singlecmds(singlecmds):
     if singlecmds[cmd] is not None:
         tmp = singlecmds[cmd].split()
         if len(tmp) != 3:
-            logger.exception(cmd + ' requires exactly three parameters')
+            logger.exception(f'{cmd} requires exactly three parameters')
             raise ValueError
 
         p1 = (float(tmp[0]), float(tmp[1]), float(tmp[2]))
@@ -143,7 +143,7 @@ def process_singlecmds(singlecmds):
     if singlecmds[cmd] is not None:
         tmp = singlecmds[cmd].split()
         if len(tmp) != 3:
-            logger.exception(cmd + ' requires exactly three parameters')
+            logger.exception(f'{cmd} requires exactly three parameters')
             raise ValueError
 
         p1 = (float(tmp[0]), float(tmp[1]), float(tmp[2]))
@@ -154,8 +154,8 @@ def process_singlecmds(singlecmds):
     cmd = '#excitation_file'
     if singlecmds[cmd] is not None:
         tmp = singlecmds[cmd].split()
-        if len(tmp) != 1 and len(tmp) != 3:
-            logger.exception(cmd + ' requires either one or three parameter(s)')
+        if len(tmp) not in [1, 3]:
+            logger.exception(f'{cmd} requires either one or three parameter(s)')
             raise ValueError
 
         if len(tmp) > 1:
diff --git a/gprMax/materials.py b/gprMax/materials.py
index abe3d310..d065f67e 100644
--- a/gprMax/materials.py
+++ b/gprMax/materials.py
@@ -213,78 +213,7 @@ class DispersiveMaterial(Material):
                 er -= ersum
 
         return er
-
-
-def process_materials(G):
-    """Processes complete list of materials - calculates update coefficients,
-        stores in arrays, and builds text list of materials/properties
-
-    Args:
-        G: FDTDGrid class describing a grid in a model.
-
-    Returns:
-        materialsdata: list of material IDs, names, and properties to
-                        print a table.
-    """
-
-    if config.get_model_config().materials['maxpoles'] == 0:
-        materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
-                          '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
-    else:
-        materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
-                          'Delta\neps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]',
-                          'gamma\n[Hz]', '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
-
-    for material in G.materials:
-        # Calculate update coefficients for specific material
-        material.calculate_update_coeffsE(G)
-        material.calculate_update_coeffsH(G)
-
-        # Add update coefficients to overall storage for all materials
-        G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
-        G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm
-
-        # Add update coefficients to overall storage for dispersive materials
-        if hasattr(material, 'poles'):
-            z = 0
-            for pole in range(config.get_model_config().materials['maxpoles']):
-                G.updatecoeffsdispersive[material.numID, z:z + 3] = (config.sim_config.em_consts['e0'] *
-                                                                     material.eqt2[pole], material.eqt[pole], material.zt[pole])
-                z += 3
-
-        # Construct information on material properties for printing table
-        materialtext = []
-        materialtext.append(str(material.numID))
-        materialtext.append(material.ID[:50] if len(material.ID) > 50 else material.ID)
-        materialtext.append(material.type)
-        materialtext.append(f'{material.er:g}')
-        materialtext.append(f'{material.se:g}')
-        if config.get_model_config().materials['maxpoles'] > 0:
-            if 'debye' in material.type:
-                materialtext.append('\n'.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
-                materialtext.append('\n'.join('{:g}'.format(tau) for tau in material.tau))
-                materialtext.extend(['', '', ''])
-            elif 'lorentz' in material.type:
-                materialtext.append(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
-                materialtext.append('')
-                materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
-                materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha))
-                materialtext.append('')
-            elif 'drude' in material.type:
-                materialtext.extend(['', ''])
-                materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
-                materialtext.append('')
-                materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha))
-            else:
-                materialtext.extend(['', '', '', '', ''])
-
-        materialtext.append(f'{material.mr:g}')
-        materialtext.append(f'{material.sm:g}')
-        materialtext.append(material.averagable)
-        materialsdata.append(materialtext)
-
-    return materialsdata
-
+    
 
 class PeplinskiSoil:
     """Soil objects that are characterised according to a mixing model
@@ -350,12 +279,11 @@ class PeplinskiSoil:
         # The limiting values of the ranges are not included in this.
 
         #mubins = np.linspace(self.mu[0], self.mu[1], nbins)
-        mubins = np.linspace(self.mu[0], self.mu[1], nbins+1)
+        mubins = np.linspace(self.mu[0], self.mu[1], nbins + 1)
         # Generate a range of volumetric water fraction values the mid-point of 
         # each bin to make materials from
         #mumaterials = mubins + (mubins[1] - mubins[0]) / 2
-        mumaterials = 0.5*(mubins[1:nbins+1] + mubins[0:nbins])
-
+        mumaterials = 0.5 * (mubins[1:nbins+1] + mubins[0:nbins])
 
         # Create an iterator
         muiter = np.nditer(mumaterials, flags=['c_index'])
@@ -399,7 +327,6 @@ class PeplinskiSoil:
             muiter.iternext()
       
 
-
 class RangeMaterial:
     """Material objects defined by a given range of their parameters to be used for 
        factal spatial disttibutions. 
@@ -425,7 +352,6 @@ class RangeMaterial:
         # and assume that all must be sequentially numbered. This allows for more general mixing models
         self.matID = [] 
        
-
     def calculate_properties(self, nbins, G):
         """Calculates the specific properties of each of the materials. 
 
@@ -436,50 +362,50 @@ class RangeMaterial:
 
         # Generate a set of relative permittivity bins based on the given range
         erbins = np.linspace(self.er[0], self.er[1], nbins+1)
+
         # Generate a range of relative permittivity values the mid-point of 
         # each bin to make materials from
         #ermaterials = erbins + np.abs((erbins[1] - erbins[0])) / 2
-        ermaterials = 0.5*(erbins[1:nbins+1] + erbins[0:nbins])
+        ermaterials = 0.5 * (erbins[1:nbins+1] + erbins[0:nbins])
 
         # Generate a set of conductivity bins based on the given range
-        sigmabins = np.linspace(self.sig[0], self.sig[1], nbins+1)
+        sigmabins = np.linspace(self.sig[0], self.sig[1], nbins + 1)
+
         # Generate a range of conductivity values the mid-point of 
         # each bin to make materials from
         #sigmamaterials = sigmabins + (sigmabins[1] - sigmabins[0]) / 2
-        sigmamaterials = 0.5*(sigmabins[1:nbins+1] + sigmabins[0:nbins])
+        sigmamaterials = 0.5 * (sigmabins[1:nbins+1] + sigmabins[0:nbins])
 
         # Generate a set of magnetic permeability bins based on the given range
-        mubins = np.linspace(self.mu[0], self.mu[1], nbins+1)
+        mubins = np.linspace(self.mu[0], self.mu[1], nbins + 1)
+
         # Generate a range of magnetic permeability values the mid-point of 
         # each bin to make materials from
         #mumaterials = mubins + np.abs((mubins[1] - mubins[0])) / 2
-        mumaterials = 0.5*(mubins[1:nbins+1] + mubins[0:nbins])
+        mumaterials = 0.5 * (mubins[1:nbins+1] + mubins[0:nbins])
         
         # Generate a set of magnetic loss bins based on the given range
-        robins = np.linspace(self.ro[0], self.ro[1], nbins+1)
-        # Generate a range of magnetic loss values the mid-point of 
-        # each bin to make materials from
+        robins = np.linspace(self.ro[0], self.ro[1], nbins + 1)
+
+        # Generate a range of magnetic loss values the mid-point of each bin to 
+        # make materials from
         #romaterials = robins + np.abs((robins[1] - robins[0])) / 2
-        romaterials = 0.5*(robins[1:nbins+1] + robins[0:nbins])
+        romaterials = 0.5 * (robins[1:nbins+1] + robins[0:nbins])
 
 
         # Iterate over the bins
-        for iter in np.arange(0,nbins):
-        
+        for iter in np.arange(nbins):
             # Relative permittivity 
             er = ermaterials[iter]
-
             # Effective conductivity
             se = sigmamaterials[iter] 
-
-            # magnetic permeability
+            # Magnetic permeability
             mr = mumaterials[iter]
-
-            # magnetic loss
+            # Magnetic loss
             sm = romaterials[iter]
 
             # Check to see if the material already exists before creating a new one
-            requiredID = '|{:.4f}+{:.4f}+{:.4f}+{:.4f}|'.format(float(er),float(se),float(mr),float(sm))
+            requiredID = f'|{float(er):.4f}+{float(se):.4f}+{float(mr):.4f}+{float(sm):.4f}|'
             material = next((x for x in G.materials if x.ID == requiredID), None)
             if iter == 0:
                 if material:
@@ -499,11 +425,10 @@ class RangeMaterial:
                 self.matID.append(m.numID)
 
 
-
 class ListMaterial:
-    """A list of predefined materials that are to be used for 
-       factal spatial disttibutions. No new materials are created but the ones specified are
-       grouped together to be used in fractal spatial distributions.
+    """A list of predefined materials to be used for 
+       factal spatial disttibutions. This command does not create new materials but collects them to be used in a 
+       stochastic distribution by a fractal box.
     """
 
     def __init__(self, ID, listofmaterials):
@@ -524,21 +449,18 @@ class ListMaterial:
         
 
     def calculate_properties(self, nbins, G):
-        """Calculates the properties of the materials. No Debye is used but name kept the same as used in other 
-           class that needs Debye
+        """Calculates the properties of the materials.
 
         Args:
             nbins: int for number of bins to use to create the different materials.
             G: FDTDGrid class describing a grid in a model.
         """
 
-       
         # Iterate over the bins
-        for iter in np.arange(0,nbins):
-        
-            # Check to see if the material already exists before creating a new one
+        for iter in np.arange(nbins):
             #requiredID = '|{:}_in_{:}|'.format((self.mat[iter]),(self.ID))
             requiredID = self.mat[iter]
+            # Check if the material already exists before creating a new one
             material = next((x for x in G.materials if x.ID == requiredID), None)
             self.matID.append(material.numID)
  
@@ -555,15 +477,11 @@ class ListMaterial:
             #    m.numID = len(G.materials)
             #    G.materials.append(m)
        
-            
             if not material:
                 logger.exception(self.__str__() + f' material(s) {material} do not exist')
                 raise ValueError
             
             
-
-
-
 def create_built_in_materials(G):
     """Creates pre-defined (built-in) materials.
 
@@ -571,8 +489,6 @@ def create_built_in_materials(G):
         G: FDTDGrid class describing a grid in a model.
     """
 
-    G.n_built_in_materials = len(G.materials)
-
     m = Material(0, 'pec')
     m.se = float('inf')
     m.type = 'builtin'
@@ -583,8 +499,6 @@ def create_built_in_materials(G):
     m.type = 'builtin'
     G.materials.append(m)
 
-    G.n_built_in_materials = len(G.materials)
-
 
 def calculate_water_properties(T=25, S=0):
     """Get extended Debye model properties for water.
@@ -627,8 +541,6 @@ def create_water(G, T=25, S=0):
 
     eri, er, tau, sig = calculate_water_properties(T, S)
     
-    G.n_built_in_materials = len(G.materials)
-
     m = DispersiveMaterial(len(G.materials), 'water')
     m.averagable = False
     m.type = 'builtin, debye'
@@ -641,8 +553,6 @@ def create_water(G, T=25, S=0):
     if config.get_model_config().materials['maxpoles'] == 0:
         config.get_model_config().materials['maxpoles'] = 1
 
-    G.n_built_in_materials = len(G.materials)
-
 
 def create_grass(G):
     """Creates single-pole Debye model for grass
@@ -657,8 +567,6 @@ def create_grass(G):
     tau = 1.0793e-11
     sig = 0
 
-    G.n_built_in_materials = len(G.materials)
-
     m = DispersiveMaterial(len(G.materials), 'grass')
     m.averagable = False
     m.type = 'builtin, debye'
@@ -671,4 +579,70 @@ def create_grass(G):
     if config.get_model_config().materials['maxpoles'] == 0:
         config.get_model_config().materials['maxpoles'] = 1
 
-    G.n_built_in_materials = len(G.materials)
\ No newline at end of file
+
+def process_materials(G):
+    """Processes complete list of materials - calculates update coefficients,
+        stores in arrays, and builds text list of materials/properties
+
+    Args:
+        G: FDTDGrid class describing a grid in a model.
+
+    Returns:
+        materialsdata: list of material IDs, names, and properties to
+                        print a table.
+    """
+
+    if config.get_model_config().materials['maxpoles'] == 0:
+        materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
+                          '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
+    else:
+        materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
+                          'Delta\neps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]',
+                          'gamma\n[Hz]', '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
+
+    for material in G.materials:
+        # Calculate update coefficients for specific material
+        material.calculate_update_coeffsE(G)
+        material.calculate_update_coeffsH(G)
+
+        # Add update coefficients to overall storage for all materials
+        G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
+        G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm
+
+        # Add update coefficients to overall storage for dispersive materials
+        if hasattr(material, 'poles'):
+            z = 0
+            for pole in range(config.get_model_config().materials['maxpoles']):
+                G.updatecoeffsdispersive[material.numID, z:z + 3] = (config.sim_config.em_consts['e0'] *
+                                                                     material.eqt2[pole], material.eqt[pole], material.zt[pole])
+                z += 3
+
+        # Construct information on material properties for printing table
+        materialtext = [str(material.numID),
+                        material.ID[:50] if len(material.ID) > 50 else material.ID,
+                        material.type,
+                        f'{material.er:g}',
+                        f'{material.se:g}']
+        if config.get_model_config().materials['maxpoles'] > 0:
+            if 'debye' in material.type:
+                materialtext.append('\n'.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
+                materialtext.append('\n'.join('{:g}'.format(tau) for tau in material.tau))
+                materialtext.extend(['', '', ''])
+            elif 'lorentz' in material.type:
+                materialtext.append(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
+                materialtext.append('')
+                materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
+                materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha))
+                materialtext.append('')
+            elif 'drude' in material.type:
+                materialtext.extend(['', ''])
+                materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
+                materialtext.append('')
+                materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha))
+            else:
+                materialtext.extend(['', '', '', '', ''])
+
+        materialtext.extend((f'{material.mr:g}', f'{material.sm:g}', material.averagable))
+        materialsdata.append(materialtext)
+
+    return materialsdata
\ No newline at end of file
diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py
index 4fd1191d..19a4753d 100644
--- a/gprMax/model_build_run.py
+++ b/gprMax/model_build_run.py
@@ -73,7 +73,7 @@ class ModelBuildRun:
         self.p = psutil.Process()
 
         # Normal model reading/building process; bypassed if geometry information to be reused
-        self.build_geometry() if not config.get_model_config().reuse_geometry else self.reuse_geometry()
+        self.reuse_geometry() if config.get_model_config().reuse_geometry else self.build_geometry()
 
         logger.info(f'\nOutput directory: {config.get_model_config().output_file_path.parent.resolve()}')
 
@@ -109,11 +109,11 @@ class ModelBuildRun:
 
         # Write files for any geometry views and geometry object outputs
         gvs = G.geometryviews + [gv for sg in G.subgrids for gv in sg.geometryviews]
-        if not (gvs or G.geometryobjectswrite) and config.sim_config.args.geometry_only:
+        if (not gvs and not G.geometryobjectswrite and config.sim_config.args.geometry_only):
             logger.exception('\nNo geometry views or geometry objects found.')
             raise ValueError
         save_geometry_views(gvs)
-        
+
         if G.geometryobjectswrite:
             logger.info('')
             for i, go in enumerate(G.geometryobjectswrite):
@@ -271,8 +271,7 @@ class ModelBuildRun:
                 logger.warning(f"You have specified more threads ({config.get_model_config().ompthreads}) "
                                f"than available physical CPU cores ({config.sim_config.hostinfo['physicalcores']}). "
                                f"This may lead to degraded performance.")
-        # Print information about any compute device, e.g. GPU, in use
-        elif config.sim_config.general['solver'] == 'cuda' or config.sim_config.general['solver'] == 'opencl':           
+        elif config.sim_config.general['solver'] in ['cuda', 'opencl']:           
             if config.sim_config.general['solver'] == 'opencl':
                 solvername = 'OpenCL'
                 platformname = ' on ' + ' '.join(config.get_model_config().device['dev'].platform.name.split())
@@ -338,11 +337,11 @@ class GridBuilder:
         pbar.close()
 
     def tm_grid_update(self):
-        if '2D TMx' == config.get_model_config().mode:
+        if config.get_model_config().mode == '2D TMx':
             self.grid.tmx()
-        elif '2D TMy' == config.get_model_config().mode:
+        elif config.get_model_config().mode == '2D TMy':
             self.grid.tmy()
-        elif '2D TMz' == config.get_model_config().mode:
+        elif config.get_model_config().mode == '2D TMz':
             self.grid.tmz()
 
     def update_voltage_source_materials(self):
diff --git a/gprMax/mpi.py b/gprMax/mpi.py
index 9e307b34..2a0dfbff 100644
--- a/gprMax/mpi.py
+++ b/gprMax/mpi.py
@@ -220,26 +220,26 @@ class MPIExecutor(object):
 
     def join(self):
         """Joins the workers."""
-        if self.is_master():
+        if not self.is_master():
+            return
+        logger.debug(f'({self.comm.name}) - Terminating. Sending sentinel to all workers.')
+        # Send sentinel to all workers
+        for worker in self.workers:
+            self.comm.send(None, dest=worker, tag=Tags.EXIT)
 
-            logger.debug(f'({self.comm.name}) - Terminating. Sending sentinel to all workers.')
-            # Send sentinel to all workers
-            for worker in self.workers:
-                self.comm.send(None, dest=worker, tag=Tags.EXIT)
+        logger.debug(f'({self.comm.name}) - Waiting for all workers to terminate.')
 
-            logger.debug(f'({self.comm.name}) - Waiting for all workers to terminate.')
+        down = [False] * len(self.workers)
+        while True:
+            for i, worker in enumerate(self.workers):
+                if self.comm.Iprobe(source=worker, tag=Tags.EXIT):
+                    self.comm.recv(source=worker, tag=Tags.EXIT)
+                    down[i] = True
+            if all(down):
+                break
 
-            down = [False] * len(self.workers)
-            while True:
-                for i, worker in enumerate(self.workers):
-                    if self.comm.Iprobe(source=worker, tag=Tags.EXIT):
-                        self.comm.recv(source=worker, tag=Tags.EXIT)
-                        down[i] = True
-                if all(down):
-                    break
-
-            self._up = False
-            logger.debug(f'({self.comm.name}) - All workers terminated.')
+        self._up = False
+        logger.debug(f'({self.comm.name}) - All workers terminated.')
 
     def submit(self, jobs, sleep=0.0):
         """Submits a list of jobs to the workers and returns the results.
diff --git a/gprMax/pml.py b/gprMax/pml.py
index d92cb69e..a320c271 100644
--- a/gprMax/pml.py
+++ b/gprMax/pml.py
@@ -22,6 +22,8 @@ import numpy as np
 
 import gprMax.config as config
 
+from .cython.pml_build import pml_average_er_mr
+
 
 class CFSParameter:
     """Individual CFS parameter (e.g. alpha, kappa, or sigma)."""
@@ -609,9 +611,6 @@ def build_pml(G, key, value):
     elif config.sim_config.general['solver'] == 'opencl':
         pml_type = OpenCLPML
 
-    sumer = 0  # Sum of relative permittivities in PML slab
-    summr = 0  # Sum of relative permeabilities in PML slab
-
     if key[0] == 'x':
         if key == 'x0':
             pml = pml_type(G, ID=key, direction='xminus', 
@@ -621,14 +620,7 @@ def build_pml(G, key, value):
                            xs=G.nx - value, xf=G.nx, yf=G.ny, zf=G.nz)
         pml.CFS = G.pmls['cfs']
         G.pmls['slabs'].append(pml)
-        for j in range(G.ny):
-            for k in range(G.nz):
-                numID = G.solid[pml.xs, j, k]
-                material = next(x for x in G.materials if x.numID == numID)
-                sumer += material.er
-                summr += material.mr
-        averageer = sumer / (G.ny * G.nz)
-        averagemr = summr / (G.ny * G.nz)
+        averageer, averagemr = pml_average_er_mr('x', pml.xs, G)
 
     elif key[0] == 'y':
         if key == 'y0':
@@ -639,14 +631,7 @@ def build_pml(G, key, value):
                            ys=G.ny - value, xf=G.nx, yf=G.ny, zf=G.nz)
         pml.CFS = G.pmls['cfs']
         G.pmls['slabs'].append(pml)
-        for i in range(G.nx):
-            for k in range(G.nz):
-                numID = G.solid[i, pml.ys, k]
-                material = next(x for x in G.materials if x.numID == numID)
-                sumer += material.er
-                summr += material.mr
-        averageer = sumer / (G.nx * G.nz)
-        averagemr = summr / (G.nx * G.nz)
+        averageer, averagemr = pml_average_er_mr('y', pml.ys, G)
 
     elif key[0] == 'z':
         if key == 'z0':
@@ -657,13 +642,6 @@ def build_pml(G, key, value):
                            zs=G.nz - value, xf=G.nx, yf=G.ny, zf=G.nz)
         pml.CFS = G.pmls['cfs']
         G.pmls['slabs'].append(pml)
-        for i in range(G.nx):
-            for j in range(G.ny):
-                numID = G.solid[i, j, pml.zs]
-                material = next(x for x in G.materials if x.numID == numID)
-                sumer += material.er
-                summr += material.mr
-        averageer = sumer / (G.nx * G.ny)
-        averagemr = summr / (G.nx * G.ny)
+        averageer, averagemr = pml_average_er_mr('z', pml.zs, G)
 
     pml.calculate_update_coeffs(averageer, averagemr)
diff --git a/gprMax/scene.py b/gprMax/scene.py
index 9f6f37a4..363cad5e 100644
--- a/gprMax/scene.py
+++ b/gprMax/scene.py
@@ -105,7 +105,7 @@ class Scene:
 
         # Check essential commands and warn user if missing
         for cmd_type in self.essential_cmds:
-            d = any([isinstance(cmd, cmd_type) for cmd in cmds_unique])
+            d = any(isinstance(cmd, cmd_type) for cmd in cmds_unique)
             if not d:
                 logger.exception('Your input file is missing essential commands ' +
                                  'required to run a model. Essential commands ' +
diff --git a/gprMax/sources.py b/gprMax/sources.py
index d44e1434..88d3ac6b 100644
--- a/gprMax/sources.py
+++ b/gprMax/sources.py
@@ -96,7 +96,7 @@ class VoltageSource(Source):
             i = self.xcoord
             j = self.ycoord
             k = self.zcoord
-            componentID = 'E' + self.polarisation
+            componentID = f'E{self.polarisation}'
 
             if self.polarisation == 'x':
                 if self.resistance != 0:
@@ -130,30 +130,31 @@ class VoltageSource(Source):
             G: FDTDGrid class describing a grid in a model.
         """
 
-        if self.resistance != 0:
-            i = self.xcoord
-            j = self.ycoord
-            k = self.zcoord
+        if self.resistance == 0:
+            return
+        i = self.xcoord
+        j = self.ycoord
+        k = self.zcoord
 
-            componentID = 'E' + self.polarisation
-            requirednumID = G.ID[G.IDlookup[componentID], i, j, k]
-            material = next(x for x in G.materials if x.numID == requirednumID)
-            newmaterial = deepcopy(material)
-            newmaterial.ID = material.ID + '+' + self.ID
-            newmaterial.numID = len(G.materials)
-            newmaterial.averagable = False
-            newmaterial.type += ',\nvoltage-source' if newmaterial.type else 'voltage-source'
+        componentID = f'E{self.polarisation}'
+        requirednumID = G.ID[G.IDlookup[componentID], i, j, k]
+        material = next(x for x in G.materials if x.numID == requirednumID)
+        newmaterial = deepcopy(material)
+        newmaterial.ID = f'{material.ID}+{self.ID}'
+        newmaterial.numID = len(G.materials)
+        newmaterial.averagable = False
+        newmaterial.type += ',\nvoltage-source' if newmaterial.type else 'voltage-source'
 
-            # Add conductivity of voltage source to underlying conductivity
-            if self.polarisation == 'x':
-                newmaterial.se += G.dx / (self.resistance * G.dy * G.dz)
-            elif self.polarisation == 'y':
-                newmaterial.se += G.dy / (self.resistance * G.dx * G.dz)
-            elif self.polarisation == 'z':
-                newmaterial.se += G.dz / (self.resistance * G.dx * G.dy)
+        # Add conductivity of voltage source to underlying conductivity
+        if self.polarisation == 'x':
+            newmaterial.se += G.dx / (self.resistance * G.dy * G.dz)
+        elif self.polarisation == 'y':
+            newmaterial.se += G.dy / (self.resistance * G.dx * G.dz)
+        elif self.polarisation == 'z':
+            newmaterial.se += G.dz / (self.resistance * G.dx * G.dy)
 
-            G.ID[G.IDlookup[componentID], i, j, k] = newmaterial.numID
-            G.materials.append(newmaterial)
+        G.ID[G.IDlookup[componentID], i, j, k] = newmaterial.numID
+        G.materials.append(newmaterial)
 
 
 class HertzianDipole(Source):
@@ -180,7 +181,7 @@ class HertzianDipole(Source):
             i = self.xcoord
             j = self.ycoord
             k = self.zcoord
-            componentID = 'E' + self.polarisation
+            componentID = f'E{self.polarisation}'
             if self.polarisation == 'x':
                 Ex[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] *
                                 self.waveformvalues_halfdt[iteration] * self.dl *
@@ -217,7 +218,7 @@ class MagneticDipole(Source):
             i = self.xcoord
             j = self.ycoord
             k = self.zcoord
-            componentID = 'H' + self.polarisation
+            componentID = f'H{self.polarisation}'
 
             if self.polarisation == 'x':
                 Hx[i, j, k] -= (updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] *
diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py
index ab679638..f5e85ece 100644
--- a/gprMax/subgrids/precursor_nodes.py
+++ b/gprMax/subgrids/precursor_nodes.py
@@ -174,7 +174,7 @@ class PrecursorNodesBase:
 
         for f in field_names:
             try:
-                val = c1 * getattr(self, f + '_0') + c2 * getattr(self, f + '_1')
+                val = c1 * getattr(self, f'{f}_0') + c2 * getattr(self, f'{f}_1')
             except ValueError:
                 raise
             setattr(self, f, val)
@@ -184,7 +184,7 @@ class PrecursorNodesBase:
         current main time step, i.e. ey_left = copy.ey_left_1
         """
         for f in field_names:
-            val = np.copy(getattr(self, f + '_1'))
+            val = np.copy(getattr(self, f'{f}_1'))
             setattr(self, f, val)
 
     def calc_exact_magnetic_in_time(self):
@@ -218,9 +218,9 @@ class PrecursorNodesBase:
 
     def update_previous_timestep_fields(self, field_names):
         for fn in field_names:
-            val = getattr(self, fn + '_1')
+            val = getattr(self, f'{fn}_1')
             val_c = np.copy(val)
-            setattr(self, fn + '_0', val_c)
+            setattr(self, f'{fn}_0', val_c)
 
     def interpolate_to_sub_grid(self, field, coords):
         x, z, x_sg, z_sg = coords
diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py
index 07b1a9c5..fc7f5d69 100644
--- a/gprMax/subgrids/updates.py
+++ b/gprMax/subgrids/updates.py
@@ -33,10 +33,10 @@ def create_updates(G):
         sg_type = type(sg)
         if sg_type == SubGridHSG and sg.filter:
             precursors = PrecursorNodesFiltered(G, sg)
-        elif sg_type == SubGridHSG and not sg.filter:
+        elif sg_type == SubGridHSG:
             precursors = PrecursorNodes(G, sg)
         else:
-            logger.exception(str(sg) + ' is not a subgrid type')
+            logger.exception(f'{str(sg)} is not a subgrid type')
             raise ValueError
 
         sgu = SubgridUpdater(sg, precursors, G)
diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py
index 1136ea01..2ea9cd64 100644
--- a/gprMax/subgrids/user_objects.py
+++ b/gprMax/subgrids/user_objects.py
@@ -45,7 +45,7 @@ class SubGridBase(UserObjectMulti):
         elif isinstance(node, UserObjectGeometry):
             self.children_geometry.append(node)
         else:
-            logger.exception(str(node) + ' this Object can not be added to a sub grid')
+            logger.exception(f'{str(node)} this Object can not be added to a sub grid')
             raise ValueError
 
     def set_discretisation(self, sg, grid):
@@ -123,12 +123,12 @@ class SubGridBase(UserObjectMulti):
         self.subgrid = sg
 
         # Copy over built in materials
-        sg.materials = [copy(m) for m in grid.materials if m.numID in range(0, grid.n_built_in_materials + 1)]
+        sg.materials = [copy(m) for m in grid.materials if m.type == 'builtin']
 
         # Don't mix and match different subgrid types
         for sg_made in grid.subgrids:
             if type(sg) != type(sg_made):
-                logger.exception(self.__str__() + ' please only use one type of subgrid')
+                logger.exception(f'{self.__str__()} please only use one type of subgrid')
                 raise ValueError
 
         # Reference the subgrid under the main grid to which it belongs
diff --git a/gprMax/updates.py b/gprMax/updates.py
index cd9bab3c..8c05f560 100644
--- a/gprMax/updates.py
+++ b/gprMax/updates.py
@@ -378,7 +378,7 @@ class CUDAUpdates:
         for pml in self.grid.pmls['slabs']:
             pml.htod_field_arrays()
             pml.set_blocks_per_grid()
-            knl_name = 'order' + str(len(pml.CFS)) + '_' + pml.direction
+            knl_name = f'order{len(pml.CFS)}_{pml.direction}'
             self.subs_name_args['FUNC'] = knl_name
 
             knl_electric = getattr(knl_pml_updates_electric, knl_name)
@@ -880,21 +880,21 @@ class OpenCLUpdates:
         for pml in self.grid.pmls['slabs']:
             pml.set_queue(self.queue)
             pml.htod_field_arrays()
-            knl_name = 'order' + str(len(pml.CFS)) + '_' + pml.direction
+            knl_name = f'order{len(pml.CFS)}_{pml.direction}'
             knl_electric_name = getattr(knl_pml_updates_electric, knl_name)
             knl_magnetic_name = getattr(knl_pml_updates_magnetic, knl_name)   
 
-            pml.update_electric_dev = self.elwise(self.ctx, 
-                                        knl_electric_name['args_opencl'].substitute({'REAL': config.sim_config.dtypes['C_float_or_double']}), 
+            pml.update_electric_dev = self.elwise(self.ctx,
+                                        knl_electric_name['args_opencl'].substitute({'REAL': config.sim_config.dtypes['C_float_or_double']}),
                                         knl_electric_name['func'].substitute(subs),
-                                        'pml_updates_electric_' + knl_name, 
+                                        f'pml_updates_electric_{knl_name}',
                                         preamble=self.knl_common,
-                                        options=config.sim_config.devices['compiler_opts'])
-            
+                                        options=config.sim_config.devices['compiler_opts'],)
+
             pml.update_magnetic_dev = self.elwise(self.ctx, 
                                         knl_magnetic_name['args_opencl'].substitute({'REAL': config.sim_config.dtypes['C_float_or_double']}), 
                                         knl_magnetic_name['func'].substitute(subs),
-                                        'pml_updates_magnetic_' + knl_name, 
+                                        f'pml_updates_magnetic_{knl_name}', 
                                         preamble=self.knl_common,
                                         options=config.sim_config.devices['compiler_opts'])
 
diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py
index 096ed74d..c920f5e7 100644
--- a/gprMax/user_inputs.py
+++ b/gprMax/user_inputs.py
@@ -65,7 +65,7 @@ class UserInput:
         except ValueError as err:
             v = ['x', 'y', 'z']
             # Discretisation
-            dl = getattr(self.grid, 'd' + err.args[0])
+            dl = getattr(self.grid, f'd{err.args[0]}')
             # Incorrect index
             i = p[v.index(err.args[0])]
             if name:
diff --git a/gprMax/utilities/host_info.py b/gprMax/utilities/host_info.py
index 9d1ecb33..ebaca218 100644
--- a/gprMax/utilities/host_info.py
+++ b/gprMax/utilities/host_info.py
@@ -508,14 +508,14 @@ def print_opencl_info(devs):
 
     import pyopencl as cl
 
-    logger.basic('|--->OpenCL:') 
+    logger.basic('|--->OpenCL:')
     logger.debug(f'PyOpenCL: {cl.VERSION_TEXT}')
 
     for i, (ID, dev) in enumerate(devs.items()):
         if i == 0:
             platform = dev.platform.name
             logger.basic(f'     |--->Platform: {platform}')
-        if not platform == dev.platform.name:
+        if platform != dev.platform.name:
             logger.basic(f'     |--->Platform: {dev.platform.name}')
         types = cl.device_type.to_string(dev.type)
         if 'CPU' in types:
diff --git a/gprMax/utilities/utilities.py b/gprMax/utilities/utilities.py
index 048aefad..c8b53bdc 100644
--- a/gprMax/utilities/utilities.py
+++ b/gprMax/utilities/utilities.py
@@ -91,7 +91,7 @@ def logo(version):
                      v""" + version + '\n\n'
 
     str = f"{description} {'=' * (get_terminal_width() - len(description) - 1)}\n\n"
-    str += Fore.CYAN + f'{logo}'
+    str += f'{Fore.CYAN}{logo}'
     str += Style.RESET_ALL + textwrap.fill(copyright, 
                                            width=get_terminal_width() - 1, 
                                            initial_indent=' ') + '\n'
diff --git a/gprMax/waveforms.py b/gprMax/waveforms.py
index 501778a3..65e52d40 100644
--- a/gprMax/waveforms.py
+++ b/gprMax/waveforms.py
@@ -54,18 +54,19 @@ class Waveform:
             waveforms.
         """
 
-        if (self.type == 'gaussian' or self.type == 'gaussiandot' or 
-            self.type == 'gaussiandotnorm' or self.type == 'gaussianprime' or 
-            self.type == 'gaussiandoubleprime'):
+        if self.type in ['gaussian',
+                         'gaussiandot',
+                         'gaussiandotnorm',
+                         'gaussianprime',
+                         'gaussiandoubleprime']:
             self.chi = 1 / self.freq
             self.zeta = 2 * np.pi**2 * self.freq**2
-        elif (self.type == 'gaussiandotdot' or 
-              self.type == 'gaussiandotdotnorm' or self.type == 'ricker'):
+        elif self.type in ['gaussiandotdot', 'gaussiandotdotnorm', 'ricker']:
             self.chi = np.sqrt(2) / self.freq
             self.zeta = np.pi**2 * self.freq**2
 
     def calculate_value(self, time, dt):
-        """Calculates value of the waveform at a specific time.
+        """Calculates the value of the waveform at a specific time.
 
         Args:
             time: float for absolute time.
@@ -82,7 +83,7 @@ class Waveform:
             delay = time - self.chi
             ampvalue = np.exp(-self.zeta * delay**2)
 
-        elif self.type == 'gaussiandot' or self.type == 'gaussianprime':
+        elif self.type in ['gaussiandot', 'gaussianprime']:
             delay = time - self.chi
             ampvalue = -2 * self.zeta * delay * np.exp(-self.zeta * delay**2)
 
@@ -91,7 +92,7 @@ class Waveform:
             normalise = np.sqrt(np.exp(1) / (2 * self.zeta))
             ampvalue = -2 * self.zeta * delay * np.exp(-self.zeta * delay**2) * normalise
 
-        elif self.type == 'gaussiandotdot' or self.type == 'gaussiandoubleprime':
+        elif self.type in ['gaussiandotdot', 'gaussiandoubleprime']:
             delay = time - self.chi
             ampvalue = (2 * self.zeta * (2 * self.zeta * delay**2 - 1) *
                         np.exp(-self.zeta * delay**2))
@@ -116,8 +117,8 @@ class Waveform:
         elif self.type == 'contsine':
             rampamp = 0.25
             ramp = rampamp * time * self.freq
-            if ramp > 1:
-                ramp = 1
+            ramp = min(ramp, 1)
+
             ampvalue = ramp * np.sin(2 * np.pi * self.freq * time)
 
         elif self.type == 'impulse':
@@ -126,7 +127,7 @@ class Waveform:
                 ampvalue = 1
             elif time >= dt:
                 ampvalue = 0
-
+                
         elif self.type == 'user':
             ampvalue = self.userfunc(time)
 
diff --git a/setup.py b/setup.py
index 3f7d9c52..8f3acf4f 100644
--- a/setup.py
+++ b/setup.py
@@ -37,8 +37,9 @@ if sys.version_info[:2] < MIN_PYTHON_VERSION:
 
 # Importing gprMax _version__.py before building can cause issues.
 with open('gprMax/_version.py', 'r') as fd:
-    version = re.search(r'^__version__\s*=\s*[\'"]([^\'"]*)[\'"]',
-                        fd.read(), re.MULTILINE).group(1)
+    version = re.search(
+        r'^__version__\s*=\s*[\'"]([^\'"]*)[\'"]', fd.read(), re.MULTILINE
+    )[1]
 
 def build_dispersive_material_templates():
     """Function to generate Cython .pyx file for dispersive media update.
@@ -46,7 +47,7 @@ def build_dispersive_material_templates():
         functions.
     """
 
-    iswin = True if sys.platform == 'win32' else False
+    iswin = (sys.platform == 'win32')
 
     env = Environment(loader = FileSystemLoader(os.path.join('gprMax', 'cython')), )
 
@@ -133,12 +134,12 @@ if 'cleanall' in sys.argv:
     for file in cythonfiles:
         filebase = os.path.splitext(file)[0]
         # Remove Cython C files
-        if os.path.isfile(filebase + '.c'):
+        if os.path.isfile(f'{filebase}.c'):
             try:
-                os.remove(filebase + '.c')
-                print(f'Removed: {filebase + ".c"}')
+                os.remove(f'{filebase}.c')
+                print(f'Removed: {filebase}.c')
             except OSError:
-                print(f'Could not remove: {filebase + ".c"}')
+                print(f'Could not remove: {filebase}.c')
         # Remove compiled Cython modules
         libfile = (glob.glob(os.path.join(os.getcwd(), 
                              os.path.splitext(file)[0]) + '*.pyd') + 
@@ -151,7 +152,7 @@ if 'cleanall' in sys.argv:
                 print(f'Removed: {os.path.abspath(libfile)}')
             except OSError:
                 print(f'Could not remove: {os.path.abspath(libfile)}')
-    
+
     # Remove build, dist, egg and __pycache__ directories
     shutil.rmtree(Path.cwd().joinpath('build'), ignore_errors=True)
     shutil.rmtree(Path.cwd().joinpath('dist'), ignore_errors=True)
@@ -163,7 +164,7 @@ if 'cleanall' in sys.argv:
     # Remove 'gprMax/cython/fields_updates_dispersive.jinja' if its there
     if os.path.isfile(cython_disp_file):
         os.remove(cython_disp_file)
-        
+
     # Now do a normal clean
     sys.argv[1] = 'clean'  # this is what distutils understands
 
@@ -176,10 +177,6 @@ else:
         linker_args = []
         libraries = []
 
-    # Compiler options - macOS - needs gcc (usually via HomeBrew) because the 
-    #                               default compiler LLVM (clang) does not 
-    #                               support OpenMP. With gcc -fopenmp option 
-    #                               implies -pthread
     elif sys.platform == 'darwin':
         # Check for Intel or Apple M series CPU
         cpuID = subprocess.check_output("sysctl -n machdep.cpu.brand_string", 
@@ -209,7 +206,7 @@ else:
                       'to be installed - easily done through the Homebrew package ' +
                       'manager (http://brew.sh). Note: gcc with OpenMP support ' +
                       'is required.')
-        
+
         # Minimum supported macOS deployment target
         MIN_MACOS_VERSION = '10.13'
         try:
@@ -219,12 +216,14 @@ else:
             pass
         os.environ['MIN_SUPPORTED_MACOSX_DEPLOYMENT_TARGET'] = MIN_MACOS_VERSION
         # Sometimes worth testing with '-fstrict-aliasing', '-fno-common'
-        compile_args = ['-O3', '-w', '-fopenmp', '-march=native', 
-                        '-mmacosx-version-min=' + MIN_MACOS_VERSION]
-        linker_args = ['-fopenmp', '-mmacosx-version-min=' + MIN_MACOS_VERSION]
+        compile_args = ['-O3',
+                        '-w',
+                        '-fopenmp',
+                        '-march=native',
+                        f'-mmacosx-version-min={MIN_MACOS_VERSION}']
+        linker_args = ['-fopenmp', f'-mmacosx-version-min={MIN_MACOS_VERSION}']
         libraries=['gomp']
 
-    # Compiler options - Linux
     elif sys.platform == 'linux':
         compile_args = ['-O3', '-w', '-fopenmp', '-march=native']
         linker_args = ['-fopenmp']
@@ -256,17 +255,17 @@ else:
     # Parse long_description from README.rst file.
     with open('README.rst','r') as fd:
         long_description = fd.read()
-        
+
     setup(name='gprMax',
           version=version,
           author='Craig Warren, Antonis Giannopoulos, and John Hartley',
           url='http://www.gprmax.com',
-          description='Electromagnetic Modelling Software based on the ' +
-                      'Finite-Difference Time-Domain (FDTD) method',
+          description='Electromagnetic Modelling Software based on the '
+          + 'Finite-Difference Time-Domain (FDTD) method',
           long_description=long_description,
           long_description_content_type="text/x-rst",
           license='GPLv3+',
-          python_requires='>' + str(MIN_PYTHON_VERSION[0]) + '.' + str(MIN_PYTHON_VERSION[1]),
+          python_requires=f'>{str(MIN_PYTHON_VERSION[0])}.{str(MIN_PYTHON_VERSION[1])}',
           install_requires=['colorama',
                             'cython',
                             'h5py',
@@ -289,4 +288,5 @@ else:
                      'Operating System :: POSIX :: Linux',
                      'Programming Language :: Cython',
                      'Programming Language :: Python :: 3',
-                     'Topic :: Scientific/Engineering'])
+                     'Topic :: Scientific/Engineering'],
+    )
diff --git a/testing/models_pmls/plot_pml_comparison.py b/testing/models_pmls/plot_pml_comparison.py
index 2acf1ecf..a713ebf4 100644
--- a/testing/models_pmls/plot_pml_comparison.py
+++ b/testing/models_pmls/plot_pml_comparison.py
@@ -93,8 +93,8 @@ for x, model in enumerate(testmodels):
     # Diffs
     datadiffs = np.zeros(datatest.shape, dtype=np.float64)
     for i in range(len(outputstest)):
-        max = np.amax(np.abs(dataref[:, i]))
-        datadiffs[:, i] = np.divide(np.abs(datatest[:, i] - dataref[:, i]), max, out=np.zeros_like(dataref[:, i]), where=max != 0)  # Replace any division by zero with zero
+        maxi = np.amax(np.abs(dataref[:, i]))
+        datadiffs[:, i] = np.divide(np.abs(datatest[:, i] - dataref[:, i]), maxi, out=np.zeros_like(dataref[:, i]), where=maxi != 0)  # Replace any division by zero with zero
 
         # Calculate power (ignore warning from taking a log of any zero values)
         with np.errstate(divide='ignore'):
diff --git a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py
index 9e080890..e87751e5 100755
--- a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py
+++ b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py
@@ -29,10 +29,16 @@ rx = gprMax.Rx(p1=(0.038, 0.114, 0.013))
 plate = gprMax.Plate(p1=(0.013, 0.013, 0.013), 
                      p2=(0.038, 0.113, 0.013), material_id='pec')
 
-gv1 = gprMax.GeometryView(p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl),
-                          filename=Path(*parts[:-1], parts[-1] + '_n'), output_type='n')
-gv2 = gprMax.GeometryView(p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl),
-                          filename=Path(*parts[:-1], parts[-1] + '_f'), output_type='f')
+gv1 = gprMax.GeometryView(p1=(0, 0, 0),
+                          p2=(x, y, z),
+                         dl=(dl, dl, dl),
+                         filename=Path(*parts[:-1], f'{parts[-1]}_n'),
+                         output_type='n',)
+gv2 = gprMax.GeometryView(p1=(0, 0, 0),
+                          p2=(x, y, z),
+                          dl=(dl, dl, dl),
+                          filename=Path(*parts[:-1], f'{parts[-1]}_f'),
+                          output_type='f',)
 
 pmls = {'CFS-PML': {'pml': gprMax.PMLProps(formulation='HORIPML', thickness=10),
                     # Parameters from http://dx.doi.org/10.1109/TAP.2018.2823864
diff --git a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py
index 1b0c8d85..f9f42c9f 100755
--- a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py
+++ b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py
@@ -29,10 +29,16 @@ rx = gprMax.Rx(p1=(0.113, 0.189, 0.088))
 plate = gprMax.Plate(p1=(0.088, 0.088, 0.088), 
                      p2=(0.113, 0.188, 0.088), material_id='pec')
 
-gv1 = gprMax.GeometryView(p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl),
-                          filename=Path(*parts[:-1], parts[-1] + '_n'), output_type='n')
-gv2 = gprMax.GeometryView(p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl),
-                          filename=Path(*parts[:-1], parts[-1] + '_f'), output_type='f')
+gv1 = gprMax.GeometryView(p1=(0, 0, 0),
+                          p2=(x, y, z),
+                          dl=(dl, dl, dl),
+                          filename=Path(*parts[:-1], f'{parts[-1]}_n'),
+                          output_type='n',)
+gv2 = gprMax.GeometryView(p1=(0, 0, 0),
+                          p2=(x, y, z),
+                          dl=(dl, dl, dl),
+                          filename=Path(*parts[:-1], f'{parts[-1]}_f'),
+                          output_type='f',)
 
 pml = gprMax.PMLProps(formulation='HORIPML', thickness=10)
 
diff --git a/testing/test_experimental.py b/testing/test_experimental.py
index e954cd3e..1146f935 100644
--- a/testing/test_experimental.py
+++ b/testing/test_experimental.py
@@ -81,8 +81,10 @@ realmax = np.where(np.abs(real[:, 1]) == 1)[0][0]
 difftime = - (timemodel[modelmax] - real[realmax, 0])
 
 # Plot modelled and real data
-fig, ax = plt.subplots(num=modelfile.stem + '_vs_' + realfile.stem, 
-                       figsize=(20, 10), facecolor='w', edgecolor='w')
+fig, ax = plt.subplots(num=f'{modelfile.stem}_vs_{realfile.stem}',
+                       figsize=(20, 10),
+                       facecolor='w',
+                       edgecolor='w',)
 ax.plot(timemodel + difftime, model, 'r', lw=2, label='Model')
 ax.plot(real[:, 0], real[:, 1], 'r', ls='--', lw=2, label='Experiment')
 ax.set_xlabel('Time [s]')
@@ -93,7 +95,7 @@ ax.legend()
 ax.grid()
 
 # Save a PDF/PNG of the figure
-savename =  modelfile.stem + '_vs_' + realfile.stem
+savename = f'{modelfile.stem}_vs_{realfile.stem}'
 savename = modelfile.parent / savename
 # fig.savefig(savename.with_suffix('.pdf'), dpi=None, format='pdf', 
 #             bbox_inches='tight', pad_inches=0.1)
diff --git a/testing/test_models.py b/testing/test_models.py
index e7e57bf2..f21bc144 100644
--- a/testing/test_models.py
+++ b/testing/test_models.py
@@ -112,11 +112,9 @@ for i, model in enumerate(testmodels):
                                      filetest.attrs['dt'], 
                                      filetest.attrs['dx_dy_dz'], rxposrelative)
 
-        filetest.close()
-
     else:
         # Get output for model and reference files
-        fileref = file.stem + '_ref'
+        fileref = f'{file.stem}_ref'
         fileref = file.parent / Path(fileref)
         fileref = h5py.File(fileref.with_suffix('.h5'), 'r')
         filetest = h5py.File(file.with_suffix('.h5'), 'r')
@@ -159,15 +157,15 @@ for i, model in enumerate(testmodels):
                 raise ValueError
 
         fileref.close()
-        filetest.close()
+    filetest.close()
 
     # Diffs
     datadiffs = np.zeros(datatest.shape, dtype=np.float64)
     for i in range(len(outputstest)):
-        max = np.amax(np.abs(dataref[:, i]))
-        datadiffs[:, i] = np.divide(np.abs(dataref[:, i] - datatest[:, i]), max, 
+        maxi = np.amax(np.abs(dataref[:, i]))
+        datadiffs[:, i] = np.divide(np.abs(dataref[:, i] - datatest[:, i]), maxi, 
                                     out=np.zeros_like(dataref[:, i]), 
-                                    where=max != 0)  # Replace any division by zero with zero
+                                    where=maxi != 0)  # Replace any division by zero with zero
 
         # Calculate power (ignore warning from taking a log of any zero values)
         with np.errstate(divide='ignore'):
@@ -188,17 +186,17 @@ for i, model in enumerate(testmodels):
                                                               facecolor='w', 
                                                               edgecolor='w')
     ex1.plot(timetest, datatest[:, 0], 'r', lw=2, label=model)
-    ex1.plot(timeref, dataref[:, 0], 'g', lw=2, ls='--', label=model + '(Ref)')
+    ex1.plot(timeref, dataref[:, 0], 'g', lw=2, ls='--', label=f'{model}(Ref)')
     ey1.plot(timetest, datatest[:, 1], 'r', lw=2, label=model)
-    ey1.plot(timeref, dataref[:, 1], 'g', lw=2, ls='--', label=model + '(Ref)')
+    ey1.plot(timeref, dataref[:, 1], 'g', lw=2, ls='--', label=f'{model}(Ref)')
     ez1.plot(timetest, datatest[:, 2], 'r', lw=2, label=model)
-    ez1.plot(timeref, dataref[:, 2], 'g', lw=2, ls='--', label=model + '(Ref)')
+    ez1.plot(timeref, dataref[:, 2], 'g', lw=2, ls='--', label=f'{model}(Ref)')
     hx1.plot(timetest, datatest[:, 3], 'r', lw=2, label=model)
-    hx1.plot(timeref, dataref[:, 3], 'g', lw=2, ls='--', label=model + '(Ref)')
+    hx1.plot(timeref, dataref[:, 3], 'g', lw=2, ls='--', label=f'{model}(Ref)')
     hy1.plot(timetest, datatest[:, 4], 'r', lw=2, label=model)
-    hy1.plot(timeref, dataref[:, 4], 'g', lw=2, ls='--', label=model + '(Ref)')
+    hy1.plot(timeref, dataref[:, 4], 'g', lw=2, ls='--', label=f'{model}(Ref)')
     hz1.plot(timetest, datatest[:, 5], 'r', lw=2, label=model)
-    hz1.plot(timeref, dataref[:, 5], 'g', lw=2, ls='--', label=model + '(Ref)')
+    hz1.plot(timeref, dataref[:, 5], 'g', lw=2, ls='--', label=f'{model}(Ref)')
     ylabels = ['$E_x$, field strength [V/m]', '$H_x$, field strength [A/m]', 
                '$E_y$, field strength [V/m]', '$H_y$, field strength [A/m]', 
                '$E_z$, field strength [V/m]', '$H_z$, field strength [A/m]']
@@ -232,7 +230,7 @@ for i, model in enumerate(testmodels):
         ax.grid()
 
     # Save a PDF/PNG of the figure
-    filediffs = file.stem + '_diffs'
+    filediffs = f'{file.stem}_diffs'
     filediffs = file.parent / Path(filediffs)
     # fig1.savefig(file.with_suffix('.pdf'), dpi=None, format='pdf', 
     #              bbox_inches='tight', pad_inches=0.1)