diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py
index cdaab566..32d2387f 100644
--- a/gprMax/cmds_geometry/add_grass.py
+++ b/gprMax/cmds_geometry/add_grass.py
@@ -16,8 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
+import logging
+
 import numpy as np
-from tqdm import tqdm
 
 import gprMax.config as config
 from .cmds_geometry import UserObjectGeometry
@@ -28,6 +29,8 @@ from ..materials import Material
 from ..utilities import round_value
 
 
+log = logging.getLogger(__name__)
+
 class AddGrass(UserObjectGeometry):
     """Allows you to add grass with roots to a :class:`gprMax.cmds_geometry.fractal_box.FractalBox` in the model.
 
@@ -68,12 +71,12 @@ class AddGrass(UserObjectGeometry):
         except KeyError:
             seed = None
 
-        # grab the correct fractal volume
+        # Get the correct fractal volume
         volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
         if volumes:
             volume = volumes[0]
         else:
-            raise CmdInputError(self.__str__() + ' Cant find FractalBox {}'.format(fractal_box_id))
+            raise CmdInputError(self.__str__() + f' cannot find FractalBox {fractal_box_id}')
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
         xs, ys, zs = p1
@@ -165,7 +168,7 @@ class AddGrass(UserObjectGeometry):
         # Set the fractal surface using the pre-calculated spatial distribution and a random height
         surface.fractalsurface = np.zeros((surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))
         for i in range(len(bladesindex[0])):
-                surface.fractalsurface[bladesindex[0][i], bladesindex[1][i]] = R.randint(surface.fractalrange[0], surface.fractalrange[1], size=1)
+            surface.fractalsurface[bladesindex[0][i], bladesindex[1][i]] = R.randint(surface.fractalrange[0], surface.fractalrange[1], size=1)
 
         # Create grass geometry parameters
         g = Grass(n_blades)
@@ -193,4 +196,4 @@ class AddGrass(UserObjectGeometry):
         volume.fractalsurfaces.append(surface)
 
         if config.is_messages():
-            tqdm.write('{} blades of grass on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(n_blades, xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, surface.dimension, surface.seed, limits[0], limits[1], surface.operatingonID))
+            tqdm.write(f'{n_blades} blades of grass on surface from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m with fractal dimension {surface.dimension:g}, fractal seeding {surface.seed}, and range {limits[0]:g}m to {limits[1]:g}m, added to {surface.operatingonID}.')
diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py
index 0e8d815e..04c99532 100644
--- a/gprMax/cmds_geometry/add_surface_roughness.py
+++ b/gprMax/cmds_geometry/add_surface_roughness.py
@@ -16,8 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
+import logging
+
 import numpy as np
-from tqdm import tqdm
 
 import gprMax.config as config
 from .cmds_geometry import UserObjectGeometry
@@ -26,6 +27,8 @@ from ..exceptions import CmdInputError
 from ..fractals import FractalSurface
 from ..utilities import round_value
 
+log = logging.getLogger(__name__)
+
 
 class AddSurfaceRoughness(UserObjectGeometry):
     """Allows you to add grass with roots to a :class:`gprMax.cmds_geometry.fractal_box.FractalBox` in the model.
@@ -62,19 +65,19 @@ class AddSurfaceRoughness(UserObjectGeometry):
             limits = np.array(self.kwargs['limits'])
             fractal_box_id = self.kwargs['fractal_box_id']
         except KeyError:
-            raise CmdInputError(self.__str__() + ' Incorrect parameters')
+            raise CmdInputError(self.__str__() + ' incorrect parameters')
 
         try:
             seed = self.kwargs['seed']
         except KeyError:
             seed = None
 
-        # grab the correct fractal volume
+        # Get the correct fractal volume
         volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
         if volumes:
             volume = volumes[0]
         else:
-            raise CmdInputError(self.__str__() + ' Cant find FractalBox {}'.format(fractal_box_id))
+            raise CmdInputError(self.__str__() + ' cannot find FractalBox {fractal_box_id}')
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
         xs, ys, zs = p1
@@ -152,10 +155,10 @@ class AddSurfaceRoughness(UserObjectGeometry):
         # List of existing surfaces IDs
         existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
         if surface.surfaceID in existingsurfaceIDs:
-            raise CmdInputError(self.__str__() + ' has already been used on the {} surface'.format(surface.surfaceID))
+            raise CmdInputError(self.__str__() + f' has already been used on the {surface.surfaceID} surface')
 
         surface.generate_fractal_surface(grid)
         volume.fractalsurfaces.append(surface)
 
         if config.is_messages():
-            tqdm.write('Fractal surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal weightings {:g}, {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, surface.dimension, surface.weighting[0], surface.weighting[1], surface.seed, limits[0], limits[1], surface.operatingonID))
+            log.info(f'Fractal surface from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m with fractal dimension {surface.dimension:g}, fractal weightings {surface.weighting[0]:g}, {surface.weighting[1]:g}, fractal seeding {surface.seed}, and range {limits[0]:g}m to {limits[1]:g}m, added to {surface.operatingonID}.')
diff --git a/gprMax/cmds_geometry/add_surface_water.py b/gprMax/cmds_geometry/add_surface_water.py
index 654e566c..e549b06e 100644
--- a/gprMax/cmds_geometry/add_surface_water.py
+++ b/gprMax/cmds_geometry/add_surface_water.py
@@ -16,7 +16,7 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-from tqdm import tqdm
+import logging
 
 import gprMax.config as config
 from .cmds_geometry import UserObjectGeometry
@@ -25,6 +25,8 @@ from ..materials import Material
 from ..utilities import round_value
 
 
+log = logging.getLogger(__name__)
+
 class AddSurfaceWater(UserObjectGeometry):
     """Allows you to add surface water to a :class:`gprMax.cmds_geometry.fractal_box.FractalBox` in the model.
 
@@ -54,12 +56,12 @@ class AddSurfaceWater(UserObjectGeometry):
         except KeyError:
             raise CmdInputError(self.__str__() + ' requires exactly eight parameters')
 
-        # grab the correct fractal volume
+        # Get the correct fractal volume
         volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
         if volumes:
             volume = volumes[0]
         else:
-            raise CmdInputError(self.__str__() + ' Cant find FractalBox {}'.format(fractal_box_id))
+            raise CmdInputError(self.__str__() + f' cannot find FractalBox {fractal_box_id}')
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
         xs, ys, zs = p1
@@ -143,4 +145,4 @@ class AddSurfaceWater(UserObjectGeometry):
             raise CmdInputError(self.__str__() + ' requires the time step for the model to be less than the relaxation time required to model water.')
 
         if config.is_messages():
-            tqdm.write('Water on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with depth {:g}m, added to {}.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, filldepth, surface.operatingonID))
+            log.info('Water on surface from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m with depth {filldepth:g}m, added to {surface.operatingonID}.')
diff --git a/gprMax/cmds_multiple.py b/gprMax/cmds_multiple.py
index f67ff313..7fe1624a 100644
--- a/gprMax/cmds_multiple.py
+++ b/gprMax/cmds_multiple.py
@@ -15,14 +15,13 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
+import logging
 import sys
+
 import numpy as np
-from tqdm import tqdm
 
 import gprMax.config as config
-from .config import z0
-from .config import dtypes
-
 from .cmds_geometry.cmds_geometry import UserObjectGeometry
 from .exceptions import CmdInputError
 from .geometry_outputs import GeometryObjects as GeometryObjectsUser
@@ -40,18 +39,16 @@ from .subgrids.base import SubGridBase
 from .utilities import round_value
 from .waveforms import Waveform as WaveformUser
 
-floattype = dtypes['float_or_double']
 
+log = logging.getLogger(__name__)
 
 class UserObjectMulti:
     """Object that can occur multiple times in a model."""
 
     def __init__(self, **kwargs):
-        """Constructor."""
         self.kwargs = kwargs
         self.order = None
         self.hash = '#example'
-        # auto translate
         self.autotranslate = True
 
 
@@ -63,7 +60,7 @@ class UserObjectMulti:
                 v = ' '.join([str(el) for el in v])
             s += str(v) + ' '
 
-        return '{}: {}'.format(self.hash, s[:-1])
+        return f'{self.hash}: {s[:-1]}'
 
     def create(self, grid, uip):
         """Create the object and add it to the grid."""
@@ -96,16 +93,15 @@ class Waveform(UserObjectMulti):
             amp = self.kwargs['amp']
             freq = self.kwargs['freq']
             ID = self.kwargs['id']
-
         except KeyError:
-            raise CmdInputError(self.params_str() + ' requires exactly four parameters')
+            raise CmdInputError(f"'{self.params_str()}' requires exactly four parameters")
 
         if wavetype not in WaveformUser.types:
-            raise CmdInputError(self.__str__() + ' must have one of the following types {}'.format(','.join(WaveformUser.types)))
+            raise CmdInputError(f"'{self.params_str()}' must have one of the following types {','.join(WaveformUser.types)}")
         if freq <= 0:
-            raise CmdInputError(self.__str__() + ' requires an excitation frequency value of greater than zero')
+            raise CmdInputError(f"'{self.params_str()}' requires an excitation frequency value of greater than zero")
         if any(x.ID == ID for x in grid.waveforms):
-            raise CmdInputError(self.__str__() + ' with ID {} already exists'.format(ID))
+            raise CmdInputError(f"'{self.params_str()}' with ID {ID} already exists")
 
         w = WaveformUser()
         w.ID = ID
@@ -113,8 +109,7 @@ class Waveform(UserObjectMulti):
         w.amp = amp
         w.freq = freq
 
-        if config.is_messages():
-            print('Waveform {} of type {} with maximum amplitude scaling {:g}, frequency {:g}Hz created.'.format(w.ID, w.type, w.amp, w.freq))
+        log.info(f'Waveform {w.ID} of type {w.type} with maximum amplitude scaling {w.amp:g}, frequency {w.freq:g}Hz created.')
 
         grid.waveforms.append(w)
 
@@ -142,34 +137,32 @@ class VoltageSource(UserObjectMulti):
         self.hash = '#voltage_source'
 
     def create(self, grid, uip):
-        """Create voltage source and add it to the grid."""
         try:
             p1 = self.kwargs['p1']
             polarisation = self.kwargs['polarisation'].lower()
             resistance = self.kwargs['resistance']
             waveform_id = self.kwargs['waveform_id']
-
         except KeyError:
-            raise CmdInputError(self.__str__() + "  requires at least six parameters")
+            raise CmdInputError(f"'{self.__str__()}' requires at least six parameters")
 
         # Check polarity & position parameters
         if polarisation not in ('x', 'y', 'z'):
-            raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be x, y, or z")
         if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be x in 2D TMx mode")
         elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be y in 2D TMy mode")
         elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'):
-            raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be z in 2D TMz mode")
 
         xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__())
 
         if resistance < 0:
-            raise CmdInputError("'{}' requires a source resistance of zero or greater".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' requires a source resistance of zero or greater")
 
         # Check if there is a waveformID in the waveforms list
         if not any(x.ID == waveform_id for x in grid.waveforms):
-            raise CmdInputError("'{}' there is no waveform with the identifier {}'.format(tmp[5]".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' there is no waveform with the identifier {tmp[5]}")
 
         v = VoltageSourceUser()
         v.polarisation = polarisation
@@ -185,17 +178,17 @@ class VoltageSource(UserObjectMulti):
             stop = self.kwargs['stop']
             # Check source start & source remove time parameters
             if start < 0:
-                raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' delay of the initiation of the source should not be less than zero")
             if stop < 0:
-                raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' time to remove the source should not be less than zero")
             if stop - start <= 0:
-                raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' duration of the source should not be zero or less")
             v.start = start
             if stop > grid.timewindow:
                 v.stop = grid.timewindow
             else:
                 v.stop = stop
-            startstop = ' start time {:g} secs, finish time {:g} secs '.format(v.start, v.stop)
+            startstop = f' start time {v.start:g} secs, finish time {v.stop:g} secs '
         except KeyError:
             v.start = 0
             v.stop = grid.timewindow
@@ -203,8 +196,7 @@ class VoltageSource(UserObjectMulti):
 
         v.calculate_waveform_values(grid)
 
-        if config.is_messages():
-            print('Voltage source with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(v.polarisation, v.xcoord * grid.dx, v.ycoord * grid.dy, v.zcoord * grid.dz, v.resistance) + startstop + 'using waveform {} created.'.format(v.waveformID))
+        log.info('Voltage source with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(v.polarisation, v.xcoord * grid.dx, v.ycoord * grid.dy, v.zcoord * grid.dz, v.resistance) + startstop + 'using waveform {} created.'.format(v.waveformID))
 
         grid.voltagesources.append(v)
 
@@ -226,36 +218,33 @@ class HertzianDipole(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 2
         self.hash = '#hertzian_dipole'
 
     def create(self, grid, uip):
-        """Create HertzianDipole and add it to the grid."""
         try:
             polarisation = self.kwargs['polarisation'].lower()
             p1 = self.kwargs['p1']
             waveform_id = self.kwargs['waveform_id']
-
         except KeyError:
-            raise CmdInputError("'{}' requires at least 3 parameters".format(self.params_str()))
+            raise CmdInputError(f"'{self.__str__()}' requires at least 3 parameters")
 
         # Check polarity & position parameters
         if polarisation not in ('x', 'y', 'z'):
-            raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be x, y, or z")
         if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be x in 2D TMx mode")
         elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be y in 2D TMy mode")
         elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'):
-            raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be z in 2D TMz mode")
 
         xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__())
 
         # Check if there is a waveformID in the waveforms list
         if not any(x.ID == waveform_id for x in grid.waveforms):
-            raise CmdInputError("'{}' there is no waveform with the identifier {}'.format(tmp[4]".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' there is no waveform with the identifier {tmp[4]}')
 
         h = HertzianDipoleUser()
         h.polarisation = polarisation
@@ -282,17 +271,17 @@ class HertzianDipole(UserObjectMulti):
             start = self.kwargs['start']
             stop = self.kwargs['stop']
             if start < 0:
-                raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' delay of the initiation of the source should not be less than zero")
             if stop < 0:
-                raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()self.__str__()}' time to remove the source should not be less than zero")
             if stop - start <= 0:
-                raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' duration of the source should not be zero or less")
             h.start = start
             if stop > grid.timewindow:
                 h.stop = grid.timewindow
             else:
                 h.stop = stop
-            startstop = ' start time {:g} secs, finish time {:g} secs '.format(h.start, h.stop)
+            startstop = f' start time {h.start:g} secs, finish time {h.stop:g} secs '
         except KeyError:
             h.start = 0
             h.stop = grid.timewindow
@@ -300,11 +289,10 @@ class HertzianDipole(UserObjectMulti):
 
         h.calculate_waveform_values(grid)
 
-        if config.is_messages():
-            if grid.mode == '2D':
-                print('Hertzian dipole is a line source in 2D with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * grid.dx, h.ycoord * grid.dy, h.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))
-            else:
-                print('Hertzian dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * grid.dx, h.ycoord * grid.dy, h.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))
+        if grid.mode == '2D':
+            log.info(f'Hertzian dipole is a line source in 2D with polarity {h.polarisation} at {h.xcoord * grid.dx:g}m, {h.ycoord * grid.dy:g}m, {h.zcoord * grid.dz:g}m,' + startstop + f'using waveform {h.waveformID} created.')
+        else:
+            log.info(f'Hertzian dipole with polarity {h.polarisation} at {h.xcoord * grid.dx:g}m, {h.ycoord * grid.dy:g}m, {h.zcoord * grid.dz:g}m,' + startstop + f'using waveform {h.waveformID} created.')
 
         grid.hertziandipoles.append(h)
 
@@ -326,35 +314,33 @@ class MagneticDipole(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 3
         self.hash = '#magnetic_dipole'
 
     def create(self, grid, uip):
-        """Create Magnetic Dipole and add it the grid."""
         try:
             polarisation = self.kwargs['polarisation'].lower()
             p1 = self.kwargs['p1']
             waveform_id = self.kwargs['waveform_id']
         except KeyError:
-            raise CmdInputError("'{}' requires at least five parameters".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' requires at least five parameters")
 
         # Check polarity & position parameters
         if polarisation not in ('x', 'y', 'z'):
-            raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be x, y, or z")
         if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be x in 2D TMx mode")
         elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be y in 2D TMy mode")
         elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'):
-            raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' polarisation must be z in 2D TMz mode")
 
         xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__())
 
         # Check if there is a waveformID in the waveforms list
         if not any(x.ID == waveform_id for x in grid.waveforms):
-            raise CmdInputError("'{}' there is no waveform with the identifier {}".format(self.__str__(), waveform_id))
+            raise CmdInputError(f"'{self.__str__()}' there is no waveform with the identifier {waveform_id}")
 
         m = MagneticDipoleUser()
         m.polarisation = polarisation
@@ -372,17 +358,17 @@ class MagneticDipole(UserObjectMulti):
             start = self.kwargs['start']
             stop = self.kwargs['stop']
             if start < 0:
-                raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' delay of the initiation of the source should not be less than zero")
             if stop < 0:
-                raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' time to remove the source should not be less than zero")
             if stop - start <= 0:
-                raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' duration of the source should not be zero or less")
             m.start = start
             if stop > grid.timewindow:
                 m.stop = grid.timewindow
             else:
                 m.stop = stop
-            startstop = ' start time {:g} secs, finish time {:g} secs '.format(m.start, m.stop)
+            startstop = f' start time {m.start:g} secs, finish time {m.stop:g} secs '
         except KeyError:
             m.start = 0
             m.stop = grid.timewindow
@@ -390,8 +376,7 @@ class MagneticDipole(UserObjectMulti):
 
         m.calculate_waveform_values(grid)
 
-        if config.is_messages():
-            print('Magnetic dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(m.polarisation, m.xcoord * grid.dx, m.ycoord * grid.dy, m.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(m.waveformID))
+        log.info(f'Magnetic dipole with polarity {m.polarisation} at {m.xcoord * grid.dx:g}m, {m.ycoord * grid.dy:g}m, {m.zcoord * grid.dz:g}m,' + startstop + f'using waveform {m.waveformID} created.')
 
         grid.magneticdipoles.append(m)
 
@@ -415,44 +400,41 @@ class TransmissionLine(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 4
         self.hash = '#transmission_line'
 
     def create(self, grid, uip):
-
         try:
             polarisation = self.kwargs['polarisation'].lower()
             p1 = self.kwargs['p1']
             waveform_id = self.kwargs['waveform_id']
             resistance = self.kwargs['resistance']
-
         except KeyError:
-            raise CmdInputError("'{}' requires at least six parameters".format(self.params_str()))
+            raise CmdInputError(f"'{self.params_str()}' requires at least six parameters")
 
         # Warn about using a transmission line on GPU
         if grid.gpu is not None:
-            raise CmdInputError("'{}' A #transmission_line cannot currently be used with GPU solving. Consider using a #voltage_source instead.".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' A #transmission_line cannot currently be used with GPU solving. Consider using a #voltage_source instead.")
 
         # Check polarity & position parameters
         if polarisation not in ('x', 'y', 'z'):
-            raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' polarisation must be x, y, or z")
         if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' polarisation must be x in 2D TMx mode")
         elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'):
-            raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' polarisation must be y in 2D TMy mode")
         elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'):
-            raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' polarisation must be z in 2D TMz mode")
 
         xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__())
 
         if resistance <= 0 or resistance >= z0:
-            raise CmdInputError("'{}' requires a resistance greater than zero and less than the impedance of free space, i.e. 376.73 Ohms".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' requires a resistance greater than zero and less than the impedance of free space, i.e. 376.73 Ohms")
 
         # Check if there is a waveformID in the waveforms list
         if not any(x.ID == waveform_id for x in grid.waveforms):
-            raise CmdInputError("'{}' there is no waveform with the identifier {}'.format(tmp[5]".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' there is no waveform with the identifier {tmp[5]}")
 
         t = TransmissionLineUser(grid)
         t.polarisation = polarisation
@@ -468,17 +450,17 @@ class TransmissionLine(UserObjectMulti):
             start = self.kwargs['start']
             stop = self.kwargs['stop']
             if start < 0:
-                raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.params_str()}' delay of the initiation of the source should not be less than zero")
             if stop < 0:
-                raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.params_str()}' time to remove the source should not be less than zero")
             if stop - start <= 0:
-                raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__()))
+                raise CmdInputError(f"'{self.params_str()}' duration of the source should not be zero or less")
             t.start = start
             if stop > grid.timewindow:
                 t.stop = grid.timewindow
             else:
                 t.stop = stop
-            startstop = ' start time {:g} secs, finish time {:g} secs '.format(t.start, t.stop)
+            startstop = f' start time {t.start:g} secs, finish time {t.stop:g} secs '
         except KeyError:
             t.start = 0
             t.stop = grid.timewindow
@@ -487,8 +469,7 @@ class TransmissionLine(UserObjectMulti):
         t.calculate_waveform_values(grid)
         t.calculate_incident_V_I(grid)
 
-        if config.is_messages():
-            print('Transmission line with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(t.polarisation, t.xcoord * grid.dx, t.ycoord * grid.dy, t.zcoord * grid.dz, t.resistance) + startstop + 'using waveform {} created.'.format(t.waveformID))
+        log.info(f'Transmission line with polarity {t.polarisation} at {t.xcoord * grid.dx:g}m, {t.ycoord * grid.dy:g}m, {t.zcoord * grid.dz:g}m, resistance {t.resistance:.1f} Ohms,' + startstop + f'using waveform {t.waveformID} created.')
 
         grid.transmissionlines.append(t)
 
@@ -508,7 +489,6 @@ class Rx(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 5
         self.hash = '#rx'
@@ -518,7 +498,7 @@ class Rx(UserObjectMulti):
         try:
             p1 = self.kwargs['p1']
         except KeyError:
-            raise CmdInputError("'{}' has an incorrect number of parameters".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' has an incorrect number of parameters")
 
         p = uip.check_src_rx_point(p1, self.__str__())
 
@@ -537,17 +517,16 @@ class Rx(UserObjectMulti):
             # Check and add field output names
             for field in outputs:
                 if field in allowableoutputs:
-                    r.outputs[field] = np.zeros(grid.iterations, dtype=floattype)
+                    r.outputs[field] = np.zeros(grid.iterations, dtype=config.dtypes['float_or_double'])
                 else:
-                    raise CmdInputError("{} contains an output type that is not allowable. Allowable outputs in current context are {}".format(self.__str__(), allowableoutputs))
-
+                    raise CmdInputError(f"'{self.__str__()}' contains an output type that is not allowable. Allowable outputs in current context are {allowableoutputs}")
         # If no ID or outputs are specified, use default
         except KeyError:
             r.ID = 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=floattype)
-        if config.is_messages():
-            print('Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs)))
+                r.outputs[key] = np.zeros(grid.iterations, dtype=config.dtypes['float_or_double'])
+
+        log.info(f"Receiver at {r.xcoord * grid.dx:g}m, {r.ycoord * grid.dy:g}m, {r.zcoord * grid.dz:g}m with output component(s) {', '.join(r.outputs)} created.")
 
         grid.rxs.append(r)
 
@@ -566,47 +545,43 @@ class RxArray(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 6
         self.hash = '#rx_array'
 
     def create(self, grid, uip):
-
         try:
             p1 = self.kwargs['p1']
             p2 = self.kwargs['p2']
             dl = self.kwargs['dl']
-
         except KeyError:
-            raise CmdInputError("'{}' requires exactly 9 parameters".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' requires exactly 9 parameters")
 
         xs, ys, zs = uip.check_src_rx_point(p1, self.__str__(), 'lower')
         xf, yf, zf = uip.check_src_rx_point(p2, self.__str__(), 'upper')
         dx, dy, dz = uip.discretise_point(dl)
 
         if xs > xf or ys > yf or zs > zf:
-            raise CmdInputError("'{}' the lower coordinates should be less than the upper coordinates".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' the lower coordinates should be less than the upper coordinates")
         if dx < 0 or dy < 0 or dz < 0:
-            raise CmdInputError("'{}' the step size should not be less than zero".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' the step size should not be less than zero")
         if dx < 1:
             if dx == 0:
                 dx = 1
             else:
-                raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation")
         if dy < 1:
             if dy == 0:
                 dy = 1
             else:
-                raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation")
         if dz < 1:
             if dz == 0:
                 dz = 1
             else:
-                raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation")
 
-        if config.is_messages():
-            print('Receiver array {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g}m'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz))
+        log.info(f'Receiver array {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m with steps {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):
@@ -620,9 +595,8 @@ class RxArray(UserObjectMulti):
                     r.zcoordorigin = z
                     r.ID = r.__class__.__name__ + '(' + str(x) + ',' + str(y) + ',' + str(z) + ')'
                     for key in RxUser.defaultoutputs:
-                        r.outputs[key] = np.zeros(grid.iterations, dtype=floattype)
-                    if config.is_messages():
-                        print('  Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs)))
+                        r.outputs[key] = np.zeros(grid.iterations, dtype=config.dtypes['float_or_double'])
+                    log.info(f"  Receiver at {r.xcoord * grid.dx:g}m, {r.ycoord * grid.dy:g}m, {r.zcoord * grid.dz:g}m with output component(s) {', '.join(r.outputs)} created.")
                     grid.rxs.append(r)
 
 
@@ -644,22 +618,20 @@ class Snapshot(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 19
         self.hash = '#snapshot'
 
     def create(self, grid, uip):
-
         if isinstance(grid, SubGridBase):
-            raise CmdInputError("'{}' Do not add Snapshots to Subgrids.".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' do not add Snapshots to Subgrids.")
         try:
             p1 = self.kwargs['p1']
             p2 = self.kwargs['p2']
             dl = self.kwargs['dl']
             filename = self.kwargs['filename']
         except KeyError:
-            raise CmdInputError("'{}' requires exactly 11 parameters".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' requires exactly 11 parameters")
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
         xs, ys, zs = p1
@@ -674,18 +646,18 @@ class Snapshot(UserObjectMulti):
             try:
                 time = self.kwargs['time']
             except KeyError:
-                raise CmdInputError("'{}' requires exactly 5 parameters".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' requires exactly 5 parameters")
             if time > 0:
                 iterations = round_value((time / grid.dt)) + 1
             else:
-                raise CmdInputError("'{}' time value must be greater than zero".format(self.__str__()))
+                raise CmdInputError(f"'{self.__str__()}' time value must be greater than zero")
 
         if dx < 0 or dy < 0 or dz < 0:
-            raise CmdInputError("'{}' the step size should not be less than zero".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' the step size should not be less than zero")
         if dx < 1 or dy < 1 or dz < 1:
-            raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation")
         if iterations <= 0 or iterations > grid.iterations:
-            raise CmdInputError("'{}' time value is not valid".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' time value is not valid")
 
         # Replace with old style snapshots if there are subgrids
         #if grid.subgrids:
@@ -694,8 +666,7 @@ class Snapshot(UserObjectMulti):
         #else:
         s = SnapshotUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, iterations, filename)
 
-        if config.is_messages():
-            print('Snapshot from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, at {:g} secs with filename {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, s.time * grid.dt, s.basefilename))
+        log.info(f'Snapshot from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m, discretisation {dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m, at {s.time * grid.dt:g} secs with filename {s.basefilename} created.')
 
         grid.snapshots.append(s)
 
@@ -714,7 +685,6 @@ class Material(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 8
         self.hash = '#material'
@@ -727,22 +697,22 @@ class Material(UserObjectMulti):
             sm = self.kwargs['sm']
             material_id = self.kwargs['id']
         except KeyError:
-            raise CmdInputError('{} requires exactly five parameters'.format(self.params_str()))
+            raise CmdInputError(f"'{self.params_str()}' requires exactly five parameters")
 
         if er < 1:
-            raise CmdInputError('{} requires a positive value of one or greater for static (DC) permittivity'.format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value of one or greater for static (DC) permittivity")
         if se != 'inf':
             se = float(se)
             if se < 0:
-                raise CmdInputError('{} requires a positive value for conductivity'.format(self.__str__()))
+                raise CmdInputError(f"'{self.params_str()}' requires a positive value for conductivity")
         else:
             se = float('inf')
         if mr < 1:
-            raise CmdInputError('{} requires a positive value of one or greater for permeability'.format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value of one or greater for permeability")
         if sm < 0:
-            raise CmdInputError('{} requires a positive value for magnetic conductivity'.format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for magnetic conductivity")
         if any(x.ID == material_id for x in grid.materials):
-            raise CmdInputError('{} with ID {} already exists'.format(material_id).format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' with ID {material_id} already exists")
 
         # Create a new instance of the Material class material (start index after pec & free_space)
         m = MaterialUser(len(grid.materials), material_id)
@@ -755,8 +725,7 @@ class Material(UserObjectMulti):
         if m.se == float('inf'):
             m.averagable = False
 
-        if config.is_messages():
-            tqdm.write('Material {} with eps_r={:g}, sigma={:g} S/m; mu_r={:g}, sigma*={:g} Ohm/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm))
+        log.info(f'Material {m.ID} with eps_r={m.er:g}, sigma={m.se:g} S/m; mu_r={m.mr:g}, sigma*={m.sm:g} Ohm/m created.')
 
         # Append the new material object to the materials list
         grid.materials.append(m)
@@ -776,31 +745,28 @@ class AddDebyeDispersion(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 9
         self.hash = '#add_dispersion_debye'
 
     def create(self, grid, uip):
-
         try:
             poles = self.kwargs['n_poles']
             er_delta = self.kwargs['er_delta']
             tau = self.kwargs['tau']
             material_ids = self.kwargs['material_ids']
-
         except KeyError:
-            raise CmdInputError(self.__str__() + ' requires at least four parameters')
+            raise CmdInputError(f"'{self.params_str()}' requires at least four parameters")
 
         if poles < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for number of poles')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for number of poles")
 
         # Look up requested materials in existing list of material instances
         materials = [y for x in material_ids for y in grid.materials if y.ID == x]
 
         if len(materials) != len(material_ids):
             notfound = [x for x in material_ids if x not in materials]
-            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+            raise CmdInputError(f"'{self.params_str()}' material(s) {notfound} do not exist")
 
         for material in materials:
             material.type = 'debye'
@@ -812,12 +778,11 @@ class AddDebyeDispersion(UserObjectMulti):
                     material.deltaer.append(er_delta[i])
                     material.tau.append(tau[i])
                 else:
-                    raise CmdInputError(self.__str__() + ' requires positive values for the permittivity difference.')
+                    raise CmdInputError(f"'{self.params_str()}' requires positive values for the permittivity difference.")
             if material.poles > MaterialUser.maxpoles:
                 MaterialUser.maxpoles = material.poles
 
-            if config.is_messages():
-                tqdm.write('Debye disperion added to {} with delta_eps_r={}, and tau={} secs created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau)))
+            log.info(f"Debye disperion added to {material.ID} with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in material.deltaer)}, and tau={', '.join('%4.3e' % tau for tau in material.tau)} secs created.")
 
 
 class AddLorentzDispersion(UserObjectMulti):
@@ -836,13 +801,11 @@ class AddLorentzDispersion(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 10
         self.hash = '#add_dispersion_lorentz'
 
     def create(self, grid, uip):
-
         try:
             poles = self.kwargs['n_poles']
             er_delta = self.kwargs['er_delta']
@@ -850,17 +813,17 @@ class AddLorentzDispersion(UserObjectMulti):
             alpha = self.kwargs['delta']
             material_ids = self.kwargs['material_ids']
         except KeyError:
-            raise CmdInputError(self.__str__() + ' requires at least five parameters')
+            raise CmdInputError(f"'{self.params_str()}' requires at least five parameters")
 
         if poles < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for number of poles')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for number of poles")
 
         # Look up requested materials in existing list of material instances
         materials = [y for x in material_ids for y in grid.materials if y.ID == x]
 
         if len(materials) != len(material_ids):
             notfound = [x for x in material_ids if x not in materials]
-            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+            raise CmdInputError(f"'{self.params_str()}' material(s) {notfound} do not exist")
 
         for material in materials:
             material.type = 'lorentz'
@@ -872,12 +835,11 @@ class AddLorentzDispersion(UserObjectMulti):
                     material.tau.append(tau[i])
                     material.alpha.append(alpha[i])
                 else:
-                    raise CmdInputError(self.__str__() + ' requires positive values for the permittivity difference and frequencies, and associated times that are greater than the time step for the model.')
+                    raise CmdInputError(f"'{self.params_str()}' requires positive values for the permittivity difference and frequencies, and associated times that are greater than the time step for the model.")
             if material.poles > MaterialUser.maxpoles:
                 MaterialUser.maxpoles = material.poles
 
-            if config.is_messages():
-                tqdm.write('Lorentz disperion added to {} with delta_eps_r={}, omega={} secs, and gamma={} created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))
+            log.info(f"Lorentz disperion added to {material.ID} with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in material.deltaer)}, omega={', '.join('%4.3e' % tau for tau in material.tau)} secs, and gamma={', '.join('%4.3e' % alpha for alpha in material.alpha)} created.")
 
 
 class AddDrudeDispersion(UserObjectMulti):
@@ -894,29 +856,28 @@ class AddDrudeDispersion(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 11
         self.hash = '#add_dispersion_Drude'
 
     def create(self, grid, uip):
-
         try:
             poles = self.kwargs['n_poles']
             tau = self.kwargs['tau']
             alpha = self.kwargs['alpha']
             material_ids = self.kwargs['material_ids']
         except KeyError:
-            raise CmdInputError(self.__str__() + ' requires at least four parameters')
+            raise CmdInputError(f"'{self.params_str()}' requires at least four parameters")
+
         if poles < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for number of poles')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for number of poles")
 
         # Look up requested materials in existing list of material instances
         materials = [y for x in material_ids for y in grid.materials if y.ID == x]
 
         if len(materials) != len(material_ids):
             notfound = [x for x in material_ids if x not in materials]
-            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+            raise CmdInputError(f"'{self.params_str()}' material(s) {notfound} do not exist")
 
         for material in materials:
             material.type = 'drude'
@@ -927,12 +888,11 @@ class AddDrudeDispersion(UserObjectMulti):
                     material.tau.append(tau[i])
                     material.alpha.append(alpha[i])
                 else:
-                    raise CmdInputError(self.__str__() + ' requires positive values for the frequencies, and associated times that are greater than the time step for the model.')
+                    raise CmdInputError(f"'{self.params_str()}' requires positive values for the frequencies, and associated times that are greater than the time step for the model.")
             if material.poles > MaterialUser.maxpoles:
                 MaterialUser.maxpoles = material.poles
 
-            if config.is_messages():
-                tqdm.write('Drude disperion added to {} with omega={} secs, and gamma={} secs created.'.format(material.ID, ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))
+            log.info(f"Drude disperion added to {material.ID} with omega={', '.join('%4.3e' % tau for tau in material.tau)} secs, and gamma={', '.join('%4.3e' % alpha for alpha in material.alpha)} secs created.")
 
 
 class SoilPeplinski(UserObjectMulti):
@@ -953,13 +913,11 @@ class SoilPeplinski(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 12
         self.hash = '#soil_peplinski'
 
     def create(self, grid, uip):
-
         try:
             sand_fraction = self.kwargs['sand_fraction']
             clay_fraction = self.kwargs['clay_fraction']
@@ -968,30 +926,28 @@ class SoilPeplinski(UserObjectMulti):
             water_fraction_lower = self.kwargs['water_fraction_lower']
             water_fraction_upper = self.kwargs['water_fraction_upper']
             ID = self.kwargs['id']
-
         except KeyError:
-            raise CmdInputError(self.__str__() + ' requires at exactly seven parameters')
+            raise CmdInputError(f"'{self.params_str()}' requires at exactly seven parameters")
 
         if sand_fraction < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for the sand fraction')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for the sand fraction")
         if clay_fraction < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for the clay fraction')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for the clay fraction")
         if bulk_density < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for the bulk density')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for the bulk density")
         if sand_density < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for the sand particle density')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for the sand particle density")
         if water_fraction_lower < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for the lower limit of the water volumetric fraction')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for the lower limit of the water volumetric fraction")
         if water_fraction_upper < 0:
-            raise CmdInputError(self.__str__() + ' requires a positive value for the upper limit of the water volumetric fraction')
+            raise CmdInputError(f"'{self.params_str()}' requires a positive value for the upper limit of the water volumetric fraction")
         if any(x.ID == ID for x in grid.mixingmodels):
-            raise CmdInputError(self.__str__() + ' with ID {} already exists'.format(ID))
+            raise CmdInputError(f"'{self.params_str()}' with ID {ID} already exists")
 
         # Create a new instance of the Material class material (start index after pec & free_space)
         s = PeplinskiSoilUser(ID, sand_fraction, clay_fraction, bulk_density, sand_density, (water_fraction_lower, water_fraction_upper))
 
-        if config.is_messages():
-            print('Mixing model (Peplinski) used to create {} with sand fraction {:g}, clay fraction {:g}, bulk density {:g}g/cm3, sand particle density {:g}g/cm3, and water volumetric fraction {:g} to {:g} created.'.format(s.ID, s.S, s.C, s.rb, s.rs, s.mu[0], s.mu[1]))
+        log.info(f'Mixing model (Peplinski) used to create {s.ID} with sand fraction {s.S:g}, clay fraction {s.C:g}, bulk density {s.rb:g}g/cm3, sand particle density {s.rs:g}g/cm3, and water volumetric fraction {s.mu[0]:g} to {s.mu[1]:g} created.')
 
         # Append the new material object to the materials list
         grid.mixingmodels.append(s)
@@ -1014,7 +970,6 @@ class GeometryView(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 18
         self.hash = '#geometry_view'
@@ -1026,9 +981,9 @@ class GeometryView(UserObjectMulti):
             self.kwargs['multi_grid']
             # there is no voxel output for multi grid output
             if isinstance(grid, SubGridBase):
-                    raise CmdInputError("'{}' Do not add multi_grid output to subgrid user object. Please add to Scene".format(self.__str__()))
+                    raise CmdInputError(f"'{self.params_str()}' do not add multi_grid output to subgrid user object. Please add to Scene")
             if output_type == 'n':
-                raise CmdInputError("'{}' Voxel output type (n) is not supported for multigrid output :(".format(self.__str__()))
+                raise CmdInputError(f"'{self.params_str()}' voxel output type (n) is not supported for multigrid output.")
             # Change constructor to the multi grid output
             from .geometry_outputs import GeometryViewFineMultiGrid as GeometryViewUser
             self.multi_grid = True
@@ -1046,7 +1001,7 @@ class GeometryView(UserObjectMulti):
             output_type = self.kwargs['output_type'].lower()
             filename = self.kwargs['filename']
         except KeyError:
-            raise CmdInputError("'{}'  requires exactly eleven parameters".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' requires exactly eleven parameters")
 
         GeometryViewUser = self.geometry_view_constructor(grid, output_type)
 
@@ -1057,15 +1012,15 @@ class GeometryView(UserObjectMulti):
         dx, dy, dz = uip.discretise_point(dl)
 
         if dx < 0 or dy < 0 or dz < 0:
-            raise CmdInputError("'{}' the step size should not be less than zero".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' the step size should not be less than zero")
         if dx > grid.nx or dy > grid.ny or dz > grid.nz:
-            raise CmdInputError("'{}' the step size should be less than the domain size".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' the step size should be less than the domain size")
         if dx < 1 or dy < 1 or dz < 1:
-            raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' the step size should not be less than the spatial discretisation")
         if output_type != 'n' and output_type != 'f':
-            raise CmdInputError("'{}' requires type to be either n (normal) or f (fine)".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' requires type to be either n (normal) or f (fine)")
         if output_type == 'f' and (dx * grid.dx != grid.dx or dy * grid.dy != grid.dy or dz * grid.dz != grid.dz):
-            raise CmdInputError("'{}' requires the spatial discretisation for the geometry view to be the same as the model for geometry view of type f (fine)".format(self.__str__()))
+            raise CmdInputError(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)")
 
         # Set type of geometry file
         if output_type == 'n':
@@ -1075,8 +1030,7 @@ class GeometryView(UserObjectMulti):
 
         g = GeometryViewUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, grid)
 
-        if config.is_messages():
-            print('Geometry view from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, multi_grid {}, grid={}, with filename base {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, self.multi_grid, grid.name, g.basefilename))
+        log.info(f'Geometry view from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m, discretisation {dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m, multi_grid={self.multi_grid}, grid={grid.name}, with filename base {g.basefilename} created.')
 
         # Append the new GeometryView object to the geometry views list
         grid.geometryviews.append(g)
@@ -1095,7 +1049,6 @@ class GeometryObjectsWrite(UserObjectMulti):
     """
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 14
         self.hash = '#geometry_objects_write'
@@ -1106,7 +1059,7 @@ class GeometryObjectsWrite(UserObjectMulti):
             p2 = self.kwargs['p2']
             filename = self.kwargs['filename']
         except KeyError:
-            raise CmdInputError("'{}' requires exactly seven parameters".format(self.__str__()))
+            raise CmdInputError(f"'{self.params_str()}' requires exactly seven parameters")
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
         x0, y0, z0 = p1
@@ -1114,8 +1067,7 @@ class GeometryObjectsWrite(UserObjectMulti):
 
         g = GeometryObjectsUser(x0, y0, z0, x1, y1, z1, filename)
 
-        if config.is_messages():
-            print('Geometry objects in the volume from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, will be written to {}, with materials written to {}'.format(p1[0] * grid.dx, p1[1] * grid.dy, p1[2] * grid.dz, p2[0] * grid.dx, p2[1] * grid.dy, p2[2] * grid.dz, g.filename, g.materialsfilename))
+        log.info(f'Geometry objects in the volume from {p1[0] * grid.dx:g}m, {p1[1] * grid.dy:g}m, {p1[2] * grid.dz:g}m, to {p2[0] * grid.dx:g}m, {p2[1] * grid.dy:g}m, {p2[2] * grid.dz:g}m, will be written to {g.filename}, with materials written to {g.materialsfilename}')
 
         # Append the new GeometryView object to the geometry objects to write list
         grid.geometryobjectswrite.append(g)
@@ -1153,16 +1105,14 @@ class PMLCFS(UserObjectMulti):
     count = 0
 
     def __init__(self, **kwargs):
-        """Constructor."""
         super().__init__(**kwargs)
         self.order = 15
         self.hash = '#pml_cfs'
         PMLCFS.count += 1
         if PMLCFS.count == 2:
-            raise CmdInputError(self.__str__() + ' can only be used up to two times, for up to a 2nd order PML')
+            raise CmdInputError(f"'{self.params_str()}' can only be used up to two times, for up to a 2nd order PML")
 
     def create(self, grid, uip):
-
         try:
             alphascalingprofile = self.kwargs['alphascalingprofile']
             alphascalingdirection = self.kwargs['alphascalingdirection']
@@ -1176,18 +1126,17 @@ class PMLCFS(UserObjectMulti):
             sigmascalingdirection = self.kwargs['sigmascalingdirection']
             sigmamin = self.kwargs['sigmamin']
             sigmamax = self.kwargs['sigmamax']
-
         except KeyError:
-            raise CmdInputError(self.__str__() + ' requires exactly twelve parameters')
+            raise CmdInputError(f"'{self.params_str()}' requires exactly twelve parameters")
 
         if alphascalingprofile not in CFSParameter.scalingprofiles.keys() or kappascalingprofile not in CFSParameter.scalingprofiles.keys() or sigmascalingprofile not in CFSParameter.scalingprofiles.keys():
-            raise CmdInputError(self.__str__() + ' must have scaling type {}'.format(','.join(CFSParameter.scalingprofiles.keys())))
+            raise CmdInputError(f"'{self.params_str()}' must have scaling type {','.join(CFSParameter.scalingprofiles.keys())}")
         if alphascalingdirection not in CFSParameter.scalingdirections or kappascalingdirection not in CFSParameter.scalingdirections or sigmascalingdirection not in CFSParameter.scalingdirections:
-            raise CmdInputError(self.__str__() + ' must have scaling type {}'.format(','.join(CFSParameter.scalingdirections)))
+            raise CmdInputError(f"'{self.params_str()}' must have scaling type {','.join(CFSParameter.scalingdirections)}")
         if float(alphamin) < 0 or float(alphamax) < 0 or float(kappamin) < 0 or float(kappamax) < 0 or float(sigmamin) < 0:
-            raise CmdInputError(self.__str__() + ' minimum and maximum scaling values must be greater than zero')
+            raise CmdInputError(f"'{self.params_str()}' minimum and maximum scaling values must be greater than zero")
         if float(kappamin) < 1:
-            raise CmdInputError(self.__str__() + ' minimum scaling value for kappa must be greater than or equal to one')
+            raise CmdInputError(f"'{self.params_str()}' minimum scaling value for kappa must be greater than or equal to one")
 
         cfsalpha = CFSParameter()
         cfsalpha.ID = 'alpha'
@@ -1215,8 +1164,7 @@ class PMLCFS(UserObjectMulti):
         cfs.kappa = cfskappa
         cfs.sigma = cfssigma
 
-        if config.is_messages():
-            print('PML CFS parameters: alpha (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), kappa (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), sigma (scaling: {}, scaling direction: {}, min: {:g}, max: {}) created.'.format(cfsalpha.scalingprofile, cfsalpha.scalingdirection, cfsalpha.min, cfsalpha.max, cfskappa.scalingprofile, cfskappa.scalingdirection, cfskappa.min, cfskappa.max, cfssigma.scalingprofile, cfssigma.scalingdirection, cfssigma.min, cfssigma.max))
+        log.info(f'PML CFS parameters: alpha (scaling: {cfsalpha.scalingprofile}, scaling direction: {cfsalpha.scalingdirection}, min: {cfsalpha.min:g}, max: {cfsalpha.max:g}), kappa (scaling: {cfskappa.scalingprofile}, scaling direction: {cfskappa.scalingdirection}, min: {cfskappa.min:g}, max: {cfskappa.max:g}), sigma (scaling: {cfssigma.scalingprofile}, scaling direction: {cfssigma.scalingdirection}, min: {cfssigma.min:g}, max: {cfssigma.max:g}) created.')
 
         grid.cfs.append(cfs)
 
@@ -1234,7 +1182,7 @@ class Subgrid(UserObjectMulti):
         elif isinstance(node, UserObjectGeometry):
             self.children_geometry.append(node)
         else:
-            raise Exception('This Object is Unknown to gprMax')
+            raise GeneralError('This object is unknown to gprMax')
 
 
 class SubgridHSG(UserObjectMulti):
diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py
index 82de288d..f92b8762 100644
--- a/gprMax/cmds_single_use.py
+++ b/gprMax/cmds_single_use.py
@@ -18,6 +18,7 @@
 
 import decimal as d
 import inspect
+import logging
 import os
 import sys
 
@@ -26,17 +27,16 @@ from colorama import Fore
 from colorama import Style
 init()
 import numpy as np
+from scipy.constants import c
 from scipy import interpolate
 
 import gprMax.config as config
-from .config import c
-from .config import dtypes
-from .config import hostinfo
 from .exceptions import CmdInputError
 from .waveforms import Waveform
 from .utilities import round_value
 
-floattype = dtypes['float_or_double']
+
+log = logging.getLogger(__name__)
 
 
 class Properties:
@@ -82,7 +82,7 @@ class Domain(UserObjectSingle):
                                            self.kwargs['p1'][1],
                                            self.kwargs['p1'][2])
         except KeyError:
-            print('error message')
+            log.warning('error message')
 
         return s
 
@@ -90,40 +90,35 @@ class Domain(UserObjectSingle):
         try:
             G.nx, G.ny, G.nz = uip.discretise_point(self.kwargs['p1'])
         except KeyError:
-            raise CmdInputError(self.__str__ + ' Please specify a point')
+            raise CmdInputError(f"'{self.params_str()}' please specify a point")
 
         if G.nx == 0 or G.ny == 0 or G.nz == 0:
-            raise CmdInputError(self.__str__ + ' requires at least one cell in every dimension')
-        if config.is_messages():
-            print('Domain size: {:g} x {:g} x {:g}m ({:d} x {:d} x {:d} = {:g} cells)'.format(self.kwargs['p1'][0], self.kwargs['p1'][1], self.kwargs['p1'][2], G.nx, G.ny, G.nz, (G.nx * G.ny * G.nz)))
+            raise CmdInputError(f"'{self.params_str()}' requires at least one cell in every dimension")
 
-        # Time step CFL limit (either 2D or 3D); switch off appropriate PMLs for 2D
+        log.info(f'Domain size: {self.kwargs['p1'][0]:g} x {self.kwargs['p1'][1]:g} x {self.kwargs['p1'][2]:g}m ({G.nx:d} x {G.ny:d} x {G.nz:d} = {(G.nx * G.ny * G.nz):g} cells)')
+
+        # Calculate time step at CFL limit; switch off appropriate PMLs for 2D
         if G.nx == 1:
-            G.dt = 1 / (c * np.sqrt((1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
+            G.calculate_dt()
             G.mode = '2D TMx'
             G.pmlthickness['x0'] = 0
             G.pmlthickness['xmax'] = 0
         elif G.ny == 1:
-            G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dz) * (1 / G.dz)))
+            G.calculate_dt()
             G.mode = '2D TMy'
             G.pmlthickness['y0'] = 0
             G.pmlthickness['ymax'] = 0
         elif G.nz == 1:
-            G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy)))
+            G.calculate_dt()
             G.mode = '2D TMz'
             G.pmlthickness['z0'] = 0
             G.pmlthickness['zmax'] = 0
         else:
-            G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
+            G.calculate_dt()
             G.mode = '3D'
 
-        # Round down time step to nearest float with precision one less than hardware maximum.
-        # Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
-        G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1)
-
-        if config.is_messages():
-            print('Mode: {}'.format(G.mode))
-            print('Time step (at CFL limit): {:g} secs'.format(G.dt))
+        log.info(f'Mode: {G.mode}')
+        log.info(f'Time step (at CFL limit): {G.dt:g} secs')
 
         # Number of threads (OpenMP) to use
         if sys.platform == 'darwin':
@@ -134,7 +129,7 @@ class Domain(UserObjectSingle):
         # os.environ['OMP_DISPLAY_ENV'] = 'TRUE' # Prints OMP version and environment variables (useful for debug)
 
         # Catch bug with Windows Subsystem for Linux (https://github.com/Microsoft/BashOnWindows/issues/785)
-        if 'Microsoft' in hostinfo['osversion']:
+        if 'Microsoft' in config.hostinfo['osversion']:
             os.environ['KMP_AFFINITY'] = 'disabled'
             del os.environ['OMP_PLACES']
             del os.environ['OMP_PROC_BIND']
@@ -143,13 +138,12 @@ class Domain(UserObjectSingle):
             G.nthreads = int(os.environ.get('OMP_NUM_THREADS'))
         else:
             # Set number of threads to number of physical CPU cores
-            G.nthreads = hostinfo['physicalcores']
+            G.nthreads = config.hostinfo['physicalcores']
             os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
 
-        if config.is_messages():
-            print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
-        if G.nthreads > hostinfo['physicalcores']:
-            print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
+        log.info(f'Number of CPU (OpenMP) threads: {G.nthreads}')
+        if G.nthreads > config.hostinfo['physicalcores']:
+            log.warning(Fore.RED + f'You have specified more threads ({G.nthreads}) than available physical CPU cores ({config.hostinfo['physicalcores']}). This may lead to degraded performance.' + Style.RESET_ALL)
 
 
 class Discretisation(UserObjectSingle):
@@ -169,7 +163,7 @@ class Discretisation(UserObjectSingle):
                                              self.kwargs['p1'][1],
                                              self.kwargs['p1'][2])
         except KeyError:
-            print('error message')
+            log.info('error message')
 
         return s
 
@@ -188,7 +182,7 @@ class Discretisation(UserObjectSingle):
             raise CmdInputError('Discretisation requires the z-direction spatial step to be greater than zero')
 
         if config.is_messages():
-            print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(*G.dl))
+            log.info('Spatial discretisation: {:g} x {:g} x {:g}m'.format(*G.dl))
 
 
 class TimeWindow(UserObjectSingle):
@@ -211,7 +205,7 @@ class TimeWindow(UserObjectSingle):
             try:
                 s = '#time_window: {}'.format(self.kwargs['iterations'])
             except KeyError:
-                print('time window error')
+                log.info('time window error')
 
         return s
 
@@ -242,7 +236,7 @@ class TimeWindow(UserObjectSingle):
             raise CmdInputError('TimeWindow: Specify a time or number of iterations')
 
         if config.is_messages():
-            print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
+            log.info('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
 
 
 class Messages(UserObjectSingle):
@@ -260,7 +254,7 @@ class Messages(UserObjectSingle):
         try:
             s = '#messages: {}'.format(self.kwargs['yn'])
         except KeyError:
-            print('messages problem')
+            log.info('messages problem')
 
     def create(self, G, uip):
         try:
@@ -296,7 +290,7 @@ class Title(UserObjectSingle):
             pass
 
         if config.is_messages():
-            print('Model title: {}'.format(G.title))
+            log.info('Model title: {}'.format(G.title))
 
 
 class NumThreads(UserObjectSingle):
@@ -332,14 +326,14 @@ class NumThreads(UserObjectSingle):
         os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
 
         if config.is_messages():
-            print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
-        if G.nthreads > hostinfo['physicalcores']:
-            print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
+            log.info('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
+        if G.nthreads > config.hostinfo['physicalcores']:
+            log.info(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, config.hostinfo['physicalcores']) + Style.RESET_ALL)
 
         # Print information about any GPU in use
         if config.is_messages():
             if G.gpu is not None:
-                print('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name))
+                log.info('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name))
 
 
 class TimeStepStabilityFactor(UserObjectSingle):
@@ -369,7 +363,7 @@ class TimeStepStabilityFactor(UserObjectSingle):
             raise CmdInputError(self.__str__() + ' requires the value of the time step stability factor to be between zero and one')
         G.dt = G.dt * f
         if config.is_messages():
-            print('Time step (modified): {:g} secs'.format(G.dt))
+            log.info('Time step (modified): {:g} secs'.format(G.dt))
 
 
 class PMLCells(UserObjectSingle):
@@ -442,7 +436,7 @@ class SrcSteps(UserObjectSingle):
             raise CmdInputError('#src_steps: requires exactly three parameters')
 
         if config.is_messages():
-            print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcsteps[0] * G.dx, G.srcsteps[1] * G.dy, G.srcsteps[2] * G.dz))
+            log.info('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcsteps[0] * G.dx, G.srcsteps[1] * G.dy, G.srcsteps[2] * G.dz))
 
 
 class RxSteps(UserObjectSingle):
@@ -463,7 +457,7 @@ class RxSteps(UserObjectSingle):
             raise CmdInputError('#rx_steps: requires exactly three parameters')
 
         if config.is_messages():
-            print('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxsteps[0] * G.dx, G.rxsteps[1] * G.dy, G.rxsteps[2] * G.dz))
+            log.info('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxsteps[0] * G.dx, G.rxsteps[1] * G.dy, G.rxsteps[2] * G.dz))
 
 
 class ExcitationFile(UserObjectSingle):
@@ -498,7 +492,7 @@ class ExcitationFile(UserObjectSingle):
                 excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
 
             if config.is_messages():
-                print('\nExcitation file: {}'.format(excitationfile))
+                log.info('\nExcitation file: {}'.format(excitationfile))
 
             # Get waveform names
             with open(excitationfile, 'r') as f:
@@ -538,7 +532,7 @@ class ExcitationFile(UserObjectSingle):
                 w.userfunc = interpolate.interp1d(waveformtime, singlewaveformvalues, **kwargs)
 
                 if config.is_messages():
-                    print('User waveform {} created using {} and, if required, interpolation parameters (kind: {}, fill value: {}).'.format(w.ID, timestr, kwargs['kind'], kwargs['fill_value']))
+                    log.info('User waveform {} created using {} and, if required, interpolation parameters (kind: {}, fill value: {}).'.format(w.ID, timestr, kwargs['kind'], kwargs['fill_value']))
 
                 G.waveforms.append(w)
 
diff --git a/gprMax/config.pxd b/gprMax/config.pxd
index 35f85fb5..5d7d64c2 100644
--- a/gprMax/config.pxd
+++ b/gprMax/config.pxd
@@ -27,9 +27,3 @@ ctypedef fused float_or_double:
 ctypedef fused float_or_double_complex:
     float complex
     double complex
-
-ctypedef fused real_or_complex:
-    float
-    double
-    float complex
-    double complex
diff --git a/gprMax/config.py b/gprMax/config.py
index 7bce848b..6783a992 100644
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -16,9 +16,7 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-import os
 from pathlib import Path
-import sys
 
 from colorama import init
 from colorama import Fore
@@ -27,31 +25,37 @@ init()
 import cython
 import numpy as np
 from scipy.constants import c
-from scipy.constants import mu_0 as m0
 from scipy.constants import epsilon_0 as e0
+from scipy.constants import mu_0 as m0
 
 from .utilities import get_host_info
 from .utilities import get_terminal_width
 
 
-# Impedance of free space (Ohms)
-z0 = np.sqrt(m0 / e0)
-
-# General setting for the simulation
+# General settings for the simulation
 #   inputfilepath: path to inputfile location
 #   outputfilepath: path to outputfile location
 #   messages: whether to print all messages as output to stdout or not
 #   progressbars: whether to show progress bars on stdoout or not
 #   mode: 2D TMx, 2D TMy, 2D TMz, or 3D
 #   cpu, cuda, opencl: solver type
+#   precision: data type for electromagnetic field output (single/double)
 #   autotranslate: auto translate objects with main grid coordinates
 #   to their equivalent local grid coordinate within the subgrid. If this option is off
 #   users must specify sub-grid object point within the global subgrid space.
-general = {'inputfilepath': 'gprMax', 'outputfilepath': 'gprMax', 'messages': True,
-           'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False, 'autotranslate': False}
+#   z0: impedance of free space (Ohms)
+general = {'messages': True,
+           'progressbars': True,
+           'mode': '3D',
+           'cpu': True,
+           'cuda': False,
+           'opencl': False,
+           'precision': 'single',
+           'autotranslate': False,
+           'z0': np.sqrt(m0 / e0)}
 
 def is_messages():
-    """Function to return messages."""
+    """Return messages."""
     return general['messages']
 
 # Store information about host machine
@@ -64,37 +68,30 @@ hostinfo = get_host_info()
 #     on the memory of the GPU. If True this will slow performance significantly
 cuda = {'gpus': None, 'snapsgpu2cpu': False}
 
-# Numerical dispersion analysis parameters
-#   highestfreqthres: threshold (dB) down from maximum power (0dB) of main frequency used
-#     to calculate highest frequency for numerical dispersion analysis
-#   maxnumericaldisp: maximum allowable percentage physical phase-velocity phase error
-#   mingridsampling: minimum grid sampling of smallest wavelength for physical wave propagation
-numdispersion = {'highestfreqthres': 40, 'maxnumericaldisp': 2, 'mingridsampling': 3}
-
-# Materials
-#   maxpoles: Maximum number of dispersive material poles in a model
-materials = {'maxpoles': 0, 'dispersivedtype': None, 'dispersiveCdtype': None}
-
-# Data types
+# Data type (precision) for electromagnetic field output
 #   Solid and ID arrays use 32-bit integers (0 to 4294967295)
 #   Rigid arrays use 8-bit integers (the smallest available type to store true/false)
-#   Fractal and dispersive coefficient arrays use complex numbers (complex)
-#                    which are represented as two floats
-#   Main field arrays use floats (float_or_double) and complex numbers (complex)
-#   Precision of float_or_double and complex: single or double for numpy and C (CUDA) arrays
-precision = 'double'
+#   Fractal arrays use complex numbers
+#   Dispersive coefficient arrays use either float or complex numbers
+#   Main field arrays use floats
 
-if precision == 'single':
-    dtypes = {'float_or_double': np.float32, 'complex': np.complex64,
-              'cython_float_or_double': cython.float, 'cython_complex': cython.floatcomplex,
-              'C_float_or_double': 'float', 'C_complex': 'pycuda::complex<float>'}
-elif precision == 'double':
-    dtypes = {'float_or_double': np.float64, 'complex': np.complex128,
-              'cython_float_or_double': cython.double, 'cython_complex': cython.doublecomplex,
-              'C_float_or_double': 'double', 'C_complex': 'pycuda::complex<double>'}
+if general['precision'] == 'single':
+    dtypes = {'float_or_double': np.float32,
+              'complex': np.complex64,
+              'cython_float_or_double': cython.float,
+              'cython_complex': cython.floatcomplex,
+              'C_float_or_double': 'float',
+              'C_complex': 'pycuda::complex<float>'}
+elif general['precision'] == 'double':
+    dtypes = {'float_or_double': np.float64,
+              'complex': np.complex128,
+              'cython_float_or_double': cython.double,
+              'cython_complex': cython.doublecomplex,
+              'C_float_or_double': 'double',
+              'C_complex': 'pycuda::complex<double>'}
 
 
-class ModelConfig():
+class ModelConfig:
     """Configuration parameters for a model.
         N.B. Multiple models can exist within a simulation
     """
@@ -107,8 +104,8 @@ class ModelConfig():
         """
 
         self.sim_config = sim_config
-        self.reuse_geometry = False
         self.i = i # Indexed from 0
+        self.reuse_geometry = False
 
         if not sim_config.single_model:
             self.appendmodelnumber = str(self.i + 1) # Indexed from 1
@@ -128,16 +125,33 @@ class ModelConfig():
         self.inputfilestr = inputfilestr_f.format(self.i + 1, self.sim_config.model_end, self.sim_config.input_file_path)
         self.next_model = Fore.GREEN + '{} {}\n'.format(self.inputfilestr, '-' * (get_terminal_width() - 1 - len(self.inputfilestr))) + Style.RESET_ALL
 
+        # Numerical dispersion analysis parameters
+        #   highestfreqthres: threshold (dB) down from maximum power (0dB) of main frequency used
+        #     to calculate highest frequency for numerical dispersion analysis
+        #   maxnumericaldisp: maximum allowable percentage physical phase-velocity phase error
+        #   mingridsampling: minimum grid sampling of smallest wavelength for physical wave propagation
+        self.numdispersion = {'highestfreqthres': 40,
+                              'maxnumericaldisp': 2,
+                              'mingridsampling': 3}
+
+        # General information to configure materials
+        #   maxpoles: Maximum number of dispersive material poles in a model
+        #   dispersivedtype: Data type for dispersive materials
+        #   dispersiveCdtype: Data type for dispersive materials in Cython
+        self.materials = {'maxpoles': 0,
+                          'dispersivedtype': None,
+                          'dispersiveCdtype': None}
+
     def get_scene(self):
         if self.sim_config.scenes:
             return self.sim_config.scenes[self.i]
         else: return None
 
     def get_usernamespace(self):
-        return {'c': c,
-                'e0': e0,
-                'm0': m0,
-                'z0': z0,
+        return {'c': c, # Speed of light in free space (m/s)
+                'e0': e0, # Permittivity of free space (F/m)
+                'm0': m0, # Permeability of free space (H/m)
+                'z0': general['z0'], # Impedance of free space (Ohms)
                 'number_model_runs': self.sim_config.model_end + 1,
                 'current_model_run': self.i + 1,
                 'inputfile': self.sim_config.input_file_path.resolve()}
@@ -155,16 +169,10 @@ class SimulationConfig:
         """
 
         self.args = args
-        self.n_models = args.n
-        self.inputfile = args.inputfile
-        self.gpu = args.gpu
-        self.mpi = args.mpi
-        self.mpi_no_spawn = args.mpi_no_spawn
-        self.general = {}
-        self.general['messages'] = general['messages']
-        self.geometry_fixed = args.geometry_fixed
-        self.geometry_only = args.geometry_only
-        self.write_processed = args.write_processed
+        self.general = general
+        self.hostinfo = hostinfo
+        self.cuda = cuda
+        self.dtypes = dtypes
 
         # Subgrid parameter may not exist if user enters via CLI
         try:
@@ -191,7 +199,7 @@ class SimulationConfig:
             self.single_model = False
 
     def set_model_start_end(self):
-        # Set range for number of models to run (internally 0 index)
+        """Set range for number of models to run (internally 0 index)."""
         if self.args.task:
             # Job array feeds args.n number of single tasks
             modelstart = self.args.task - 1
@@ -206,23 +214,20 @@ class SimulationConfig:
         self.model_start = modelstart
         self.model_end = modelend
 
-    def set_precision(self):
-        pass
-
     def set_input_file_path(self):
-        # If the API is in use an id for the simulation must be provided.
+        """If the API is in use an id for the simulation must be provided."""
         if self.args.inputfile is None:
             self.input_file_path = Path(self.args.outputfile)
         else:
             self.input_file_path = Path(self.args.inputfile)
 
     def set_output_file_path(self):
-        # Output file path can be provided by the user. If they havent provided one
-        # use the inputfile file path instead
+        """Output file path can be provided by the user. If they havent provided one
+            use the inputfile file path instead."""
         try:
             self.output_file_path = Path(self.args.outputfile)
         except AttributeError:
-            self.output_file_path = Path(self.args.inputfile)
+            self.output_file_path = Path(self.input_file_path)
 
 
 class SimulationConfigMPI(SimulationConfig):
diff --git a/gprMax/contexts.py b/gprMax/contexts.py
index ed2e7f1e..7ca31ba0 100644
--- a/gprMax/contexts.py
+++ b/gprMax/contexts.py
@@ -18,15 +18,17 @@
 
 import datetime
 
+from ._version import __version__, codename
 from .config import create_model_config
 from .model_build_run import ModelBuildRun
 from .solvers import create_solver
 from .solvers import create_G
 from .utilities import get_terminal_width
+from .utilities import logo
 from .utilities import timer
 
 
-class Context():
+class Context:
     """Generic context for the model to run in. Sub-class with specific contexts
     e.g. an MPI context.
     """
@@ -43,18 +45,23 @@ class Context():
         self.tsimstart = 1
 
     def run(self):
-        """Function to run the simulation in the correct context."""
+        """Run the simulation in the correct context."""
+        self.print_logo_copyright()
         self.tsimstart = timer()
         self._run()
         self.tsimend = timer()
 
+    def print_logo_copyright(self):
+        """Print gprMax logo, version, and copyright/licencing information."""
+        logo(__version__ + ' (' + codename + ')')
+
     def print_time_report(self):
-        """Function to print the total simulation time based on context."""
-        s = self.make_time_report(sim_time)
-        print(s)
+        """Print the total simulation time based on context."""
+        s = self.make_time_report()
+        log.info(s)
 
     def make_time_report(self):
-        """Function to generate a string for the total simulation time."""
+        """Generate a string for the total simulation time."""
         pass
 
 
@@ -63,6 +70,9 @@ class NoMPIContext(Context):
         is parallelised using either OpenMP (CPU) or CUDA (GPU).
     """
 
+    def __init__(self):
+        super().__init__()
+
     def _run(self):
         """Specialise how the models are farmed out."""
 
@@ -83,7 +93,7 @@ class NoMPIContext(Context):
             solver = create_solver(G, self.sim_config)
 
             if not self.sim_config.geometry_only:
-                model.run_model(solver)
+                model.solve(solver)
 
     def make_time_report(self):
         """Function to specialise the time reporting for the standard Simulation
@@ -91,9 +101,8 @@ class NoMPIContext(Context):
         """
 
         sim_time = datetime.timedelta(seconds=self.tsimend - self.tsimstart)
-        s = '\n=== Simulation on {} completed in [HH:MM:SS]: {}'
-        s = s.format(self.simconfig.hostinfo['hostname'], sim_time)
-        return '{} {}\n'.format(s, '=' * (get_terminal_width() - 1 - len(s)))
+        s = f'\n=== Simulation on {self.simconfig.hostinfo['hostname']} completed in [HH:MM:SS]: {sim_time}'
+        return f'{s} {'=' * (get_terminal_width() - 1 - len(s))}\n'
 
 
 class MPIContext(Context):
@@ -116,6 +125,9 @@ class MPIContext(Context):
 
 class MPINoSpawnContext(Context):
 
+    def __init__(self):
+        super().__init__()
+
     def _run(self):
         pass
 
@@ -133,9 +145,9 @@ def create_context(sim_config):
         context (Context): Context for the model to run in.
     """
 
-    if sim_config.mpi:
+    if sim_config.args.mpi:
         context = MPIContext(sim_config)
-    elif sim_config.mpi_no_spawn:
+    elif sim_config.args.mpi_no_spawn:
         context = MPINoSpawnContext(sim_config)
     else:
         context = NoMPIContext(sim_config)
diff --git a/gprMax/cython/fields_updates_hsg.pyx b/gprMax/cython/fields_updates_hsg.pyx
index ea0379ff..fa58bc05 100644
--- a/gprMax/cython/fields_updates_hsg.pyx
+++ b/gprMax/cython/fields_updates_hsg.pyx
@@ -1,4 +1,4 @@
-# Copyright (C) 2015-2017: The University of Edinburgh
+# Copyright (C) 2015-2019: The University of Edinburgh
 #                 Authors: Craig Warren and Antonis Giannopoulos
 #
 # This file is part of gprMax.
@@ -15,10 +15,13 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
 import numpy as np
 cimport numpy as np
+
 from cython.parallel import prange
 
+
 cpdef void cython_update_electric_os(
         np.float64_t[:, :] updatecoeffsE,
         np.uint32_t[:, :, :, :] ID,
diff --git a/gprMax/cython/fields_updates_normal.pyx b/gprMax/cython/fields_updates_normal.pyx
index ccd65ee1..771b7f0d 100644
--- a/gprMax/cython/fields_updates_normal.pyx
+++ b/gprMax/cython/fields_updates_normal.pyx
@@ -1,4 +1,4 @@
-cdef float phi# Copyright (C) 2015-2019: The University of Edinburgh
+# Copyright (C) 2015-2019: The University of Edinburgh
 #                 Authors: Craig Warren and Antonis Giannopoulos
 #
 # This file is part of gprMax.
diff --git a/gprMax/cython/geometry_outputs.pyx b/gprMax/cython/geometry_outputs.pyx
index 21490909..994afcd6 100644
--- a/gprMax/cython/geometry_outputs.pyx
+++ b/gprMax/cython/geometry_outputs.pyx
@@ -16,6 +16,7 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
+import numpy as np
 cimport numpy as np
 
 
diff --git a/gprMax/cython/geometry_primitives.pyx b/gprMax/cython/geometry_primitives.pyx
index d6e65532..586ff1db 100644
--- a/gprMax/cython/geometry_primitives.pyx
+++ b/gprMax/cython/geometry_primitives.pyx
@@ -18,8 +18,8 @@
 
 import numpy as np
 cimport numpy as np
+np.seterr(divide='raise')
 
-from gprMax.utilities import round_value
 from gprMax.cython.yee_cell_setget_rigid cimport set_rigid_Ex
 from gprMax.cython.yee_cell_setget_rigid cimport set_rigid_Ey
 from gprMax.cython.yee_cell_setget_rigid cimport set_rigid_Ez
@@ -30,8 +30,7 @@ from gprMax.cython.yee_cell_setget_rigid cimport set_rigid_E
 from gprMax.cython.yee_cell_setget_rigid cimport unset_rigid_E
 from gprMax.cython.yee_cell_setget_rigid cimport set_rigid_H
 from gprMax.cython.yee_cell_setget_rigid cimport unset_rigid_H
-
-np.seterr(divide='raise')
+from gprMax.utilities import round_value
 
 
 cpdef bint are_clockwise(
diff --git a/gprMax/cython/pml_updates_electric_HORIPML.pyx b/gprMax/cython/pml_updates_electric_HORIPML.pyx
index 5bc9ea9e..3c13021c 100644
--- a/gprMax/cython/pml_updates_electric_HORIPML.pyx
+++ b/gprMax/cython/pml_updates_electric_HORIPML.pyx
@@ -4,16 +4,16 @@
 # This file is part of gprMax.
 #
 # gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU GenRAl Public License as published by
+# 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 GenRAl Public License for more details.
+# GNU General Public License for more details.
 #
-# You should have received a copy of the GNU GenRAl Public License
+# You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import numpy as np
diff --git a/gprMax/cython/pml_updates_electric_MRIPML.pyx b/gprMax/cython/pml_updates_electric_MRIPML.pyx
index 805a36ea..c2ddbfde 100644
--- a/gprMax/cython/pml_updates_electric_MRIPML.pyx
+++ b/gprMax/cython/pml_updates_electric_MRIPML.pyx
@@ -1,19 +1,19 @@
-cdef float_or_double# Copyright (C) 2015-2019: The University of Edinburgh
+# Copyright (C) 2015-2019: The University of Edinburgh
 #                 Authors: Craig Warren and Antonis Giannopoulos
 #
 # This file is part of gprMax.
 #
 # gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU GenRAl Public License as published by
+# 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 GenRAl Public License for more details.
+# GNU General Public License for more details.
 #
-# You should have received a copy of the GNU GenRAl Public License
+# You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import numpy as np
diff --git a/gprMax/cython/snapshots.pyx b/gprMax/cython/snapshots.pyx
index 05a9f516..12332329 100644
--- a/gprMax/cython/snapshots.pyx
+++ b/gprMax/cython/snapshots.pyx
@@ -16,8 +16,6 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-cimport numpy as np
-
 from gprMax.config cimport float_or_double
 
 
diff --git a/gprMax/cython/yee_cell_build.pyx b/gprMax/cython/yee_cell_build.pyx
index 39202e41..5c812b61 100644
--- a/gprMax/cython/yee_cell_build.pyx
+++ b/gprMax/cython/yee_cell_build.pyx
@@ -19,13 +19,13 @@
 import numpy as np
 cimport numpy as np
 
-from gprMax.materials import Material
 from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Ex
 from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Ey
 from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Ez
 from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Hx
 from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Hy
 from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Hz
+from gprMax.materials import Material
 
 
 cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2, int numID3, int numID4, int componentID, G):
diff --git a/gprMax/cython/yee_cell_setget_rigid.pxd b/gprMax/cython/yee_cell_setget_rigid.pxd
index 818113e0..df70ba6a 100644
--- a/gprMax/cython/yee_cell_setget_rigid.pxd
+++ b/gprMax/cython/yee_cell_setget_rigid.pxd
@@ -19,8 +19,9 @@
 import numpy as np
 cimport numpy as np
 
-# Get and set functions for the rigid electric component array. The rigid array is 4D with the 1st dimension holding
-# the 12 electric edge components of a cell - Ex1, Ex2, Ex3, Ex4, Ey1, Ey2, Ey3, Ey4, Ez1, Ez2, Ez3, Ez4
+# Get and set functions for the rigid electric component array. The rigid array
+# is 4D with the 1st dimension holding the 12 electric edge components of a
+# cell - Ex1, Ex2, Ex3, Ex4, Ey1, Ey2, Ey3, Ey4, Ez1, Ez2, Ez3, Ez4
 cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
 cdef bint get_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
 cdef bint get_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
@@ -30,8 +31,9 @@ cdef void set_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
 cdef void set_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
 cdef void unset_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
 
-# Get and set functions for the rigid magnetic component array. The rigid array is 4D with the 1st dimension holding
-# the 6 magnetic edge components - Hx1, Hx2, Hy1, Hy2, Hz1, Hz2
+# Get and set functions for the rigid magnetic component array. The rigid array
+# is 4D with the 1st dimension holding the 6 magnetic edge components - Hx1,
+# Hx2, Hy1, Hy2, Hz1, Hz2
 cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
 cdef bint get_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
 cdef bint get_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
@@ -40,5 +42,3 @@ cdef void set_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
 cdef void set_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
 cdef void set_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
 cdef void unset_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-
-
diff --git a/gprMax/cython/yee_cell_setget_rigid.pyx b/gprMax/cython/yee_cell_setget_rigid.pyx
index 707fff52..6858b7e8 100644
--- a/gprMax/cython/yee_cell_setget_rigid.pyx
+++ b/gprMax/cython/yee_cell_setget_rigid.pyx
@@ -19,8 +19,10 @@
 import numpy as np
 cimport numpy as np
 
-# Get and set functions for the rigid electric component array. The rigid array is 4D with the 1st dimension holding
-# the 12 electric edge components of a cell - Ex1, Ex2, Ex3, Ex4, Ey1, Ey2, Ey3, Ey4, Ez1, Ez2, Ez3, Ez4
+
+# Get and set functions for the rigid electric component array. The rigid array
+# is 4D with the 1st dimension holding the 12 electric edge components of a
+# cell - Ex1, Ex2, Ex3, Ex4, Ey1, Ey2, Ey3, Ey4, Ez1, Ez2, Ez3, Ez4
 cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
     cdef bint result
     result = False
@@ -102,8 +104,9 @@ cdef void set_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
 cdef void unset_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
     rigidE[:, i, j, k] = False
 
-# Get and set functions for the rigid magnetic component array. The rigid array is 4D with the 1st dimension holding
-# the 6 magnetic edge components - Hx1, Hx2, Hy1, Hy2, Hz1, Hz2
+# Get and set functions for the rigid magnetic component array. The rigid array
+# is 4D with the 1st dimension holding the 6 magnetic edge components - Hx1,
+# Hx2, Hy1, Hy2, Hz1, Hz2
 cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
     cdef bint result
     result = False
@@ -154,4 +157,3 @@ cdef void set_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
 
 cdef void unset_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
     rigidH[:, i, j, k] = False
-
diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py
index fcbc6375..7f8d1963 100644
--- a/gprMax/gprMax.py
+++ b/gprMax/gprMax.py
@@ -17,11 +17,14 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import argparse
+import logging
 
 from .config import create_simulation_config
 from .contexts import create_context
 
 
+logging.basicConfig(level=logging.INFO)
+
 def run(
     scenes=None,
     inputfile=None,
@@ -137,3 +140,4 @@ def run_main(args):
     sim_config = create_simulation_config(args)
     context = create_context(sim_config)
     context.run()
+    context.print_time_report()
diff --git a/gprMax/grid.py b/gprMax/grid.py
index 643cf425..d6b9f784 100644
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -36,80 +36,16 @@ from .utilities import human_size
 from .utilities import round_value
 
 
-class Grid(object):
-    """Generic grid/mesh."""
-
-    def __init__(self, grid):
-        self.nx = grid.shape[0]
-        self.ny = grid.shape[1]
-        self.nz = grid.shape[2]
-        self.dx = 1
-        self.dy = 1
-        self.dz = 1
-        self.i_max = self.nx - 1
-        self.j_max = self.ny - 1
-        self.k_max = self.nz - 1
-        self.grid = grid
-
-    def n_edges(self):
-        i = self.nx
-        j = self.ny
-        k = self.nz
-        e = (i * j * (k - 1)) + (j * k * (i - 1)) + (i * k * (j - 1))
-        return e
-
-    def n_nodes(self):
-        return self.nx * self.ny * self.nz
-
-    def n_cells(self):
-        return (self.nx - 1) * (self.ny - 1) * (self.nz - 1)
-
-    def get(self, i, j, k):
-        return self.grid[i, j, k]
-
-    def within_bounds(self, p):
-        if p[0] < 0 or p[0] > self.nx:
-            raise ValueError('x')
-        if p[1] < 0 or p[1] > self.ny:
-            raise ValueError('y')
-        if p[2] < 0 or p[2] > self.nz:
-            raise ValueError('z')
-
-    def discretise_point(self, p):
-        x = round_value(float(p[0]) / self.dx)
-        y = round_value(float(p[1]) / self.dy)
-        z = round_value(float(p[2]) / self.dz)
-        return (x, y, z)
-
-    def round_to_grid(self, p):
-        p = self.discretise_point(p)
-        p_r = (p[0] * self.dx,
-               p[1] * self.dy,
-               p[2] * self.dz)
-        return p_r
-
-
-    def within_pml(self, p):
-        if (p[0] < self.pmlthickness['x0'] or
-            p[0] > self.nx - self.pmlthickness['xmax'] or
-            p[1] < self.pmlthickness['y0'] or
-            p[1] > self.ny - self.pmlthickness['ymax'] or
-            p[2] < self.pmlthickness['z0'] or
-            p[2] > self.nz - self.pmlthickness['zmax']):
-            return True
-        else:
-            return False
-
-
-class FDTDGrid(Grid):
-    """Holds attributes associated with entire grid. A convenient
-        way for accessing regularly used parameters.
+class FDTDGrid:
+    """Holds attributes associated with entire grid. A convenient way for
+        accessing regularly used parameters.
     """
 
     def __init__(self):
         self.title = ''
         self.memoryusage = 0
         self.name = 'Main'
+        self.mode = '' # 2D TMx, 2D TMy, 2D TMz, or 3D
         self.gpu = None
         self.outputdirectory = ''
 
@@ -150,6 +86,50 @@ class FDTDGrid(Grid):
         self.snapshots = []
         self.subgrids = []
 
+    def n_edges(self):
+        i = self.nx
+        j = self.ny
+        k = self.nz
+        e = (i * j * (k - 1)) + (j * k * (i - 1)) + (i * k * (j - 1))
+        return e
+
+    def n_nodes(self):
+        return self.nx * self.ny * self.nz
+
+    def n_cells(self):
+        return (self.nx - 1) * (self.ny - 1) * (self.nz - 1)
+
+    def within_bounds(self, p):
+        if p[0] < 0 or p[0] > self.nx:
+            raise ValueError('x')
+        if p[1] < 0 or p[1] > self.ny:
+            raise ValueError('y')
+        if p[2] < 0 or p[2] > self.nz:
+            raise ValueError('z')
+
+    def discretise_point(self, p):
+        x = round_value(float(p[0]) / self.dx)
+        y = round_value(float(p[1]) / self.dy)
+        z = round_value(float(p[2]) / self.dz)
+        return (x, y, z)
+
+    def round_to_grid(self, p):
+        p = self.discretise_point(p)
+        p_r = (p[0] * self.dx,
+               p[1] * self.dy,
+               p[2] * self.dz)
+        return p_r
+
+    def within_pml(self, p):
+        if (p[0] < self.pmlthickness['x0'] or
+            p[0] > self.nx - self.pmlthickness['xmax'] or
+            p[1] < self.pmlthickness['y0'] or
+            p[1] > self.ny - self.pmlthickness['ymax'] or
+            p[2] < self.pmlthickness['z0'] or
+            p[2] > self.nz - self.pmlthickness['zmax']):
+            return True
+        else:
+            return False
 
     def initialise_geometry_arrays(self):
         """Initialise an array for volumetric material IDs (solid);
@@ -173,9 +153,8 @@ class FDTDGrid(Grid):
         self.Hy = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.dtypes['float_or_double'])
         self.Hz = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.dtypes['float_or_double'])
 
-
     def initialise_grids(self):
-        """Function to call the initialisation of all grids."""
+        """Initialise all grids."""
         for g in [self] + self.subgrids:
             g.initialise_geometry_arrays()
             g.initialise_field_arrays()
@@ -227,32 +206,75 @@ class FDTDGrid(Grid):
 
         self.memoryusage = int(stdoverhead + fieldarrays + solidarray + rigidarrays + pmlarrays)
 
-    def memory_check(self, snapsmemsize=0):
-        """Check if the required amount of memory (RAM) is available on the host and GPU if specified.
-
-        Args:
-            snapsmemsize (int): amount of memory (bytes) required to store all requested snapshots
+    def memory_check(self):
+        """Check if the required amount of memory (RAM) is available to build
+            and/or run model on the host.
         """
-
-        # Check if model can be built and/or run on host
         if self.memoryusage > config.hostinfo['ram']:
-            raise GeneralError('Memory (RAM) required ~{} exceeds {} detected!\n'.format(human_size(self.memoryusage), human_size(config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)))
+            raise GeneralError(f'Memory (RAM) required ~{human_size(self.memoryusage)} exceeds {human_size(config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)} detected!\n')
 
-        # Check if model can be run on specified GPU if required
-        if config.cuda['gpus'] is not None:
-            if self.memoryusage - snapsmemsize > config.cuda['gpus'].totalmem:
-                raise GeneralError('Memory (RAM) required ~{} exceeds {} detected on specified {} - {} GPU!\n'.format(human_size(self.memoryusage), human_size(config.cuda['gpus'].totalmem, a_kilobyte_is_1024_bytes=True), config.cuda['gpus'].deviceID, config.cuda['gpus'].name))
+    def tmx(self):
+        """Add PEC boundaries to invariant direction in 2D TMx mode.
+            N.B. 2D modes are a single cell slice of 3D grid.
+        """
+        # Ey & Ez components
+        self.ID[1, 0, :, :] = 0
+        self.ID[1, 1, :, :] = 0
+        self.ID[2, 0, :, :] = 0
+        self.ID[2, 1, :, :] = 0
 
-            # If the required memory without the snapshots will fit on the GPU then transfer and store snaphots on host
-            if snapsmemsize != 0 and self.memoryusage - snapsmemsize < config.cuda['gpus'].totalmem:
-                config.cuda['snapsgpu2cpu'] = True
+    def tmy(self):
+        """Add PEC boundaries to invariant direction in 2D TMy mode.
+            N.B. 2D modes are a single cell slice of 3D grid.
+        """
+        # Ex & Ez components
+        self.ID[0, :, 0, :] = 0
+        self.ID[0, :, 1, :] = 0
+        self.ID[2, :, 0, :] = 0
+        self.ID[2, :, 1, :] = 0
 
-    def gpu_set_blocks_per_grid(self):
+    def tmz(self):
+        """Add PEC boundaries to invariant direction in 2D TMz mode.
+            N.B. 2D modes are a single cell slice of 3D grid.
+        """
+        # Ex & Ey components
+        self.ID[0, :, :, 0] = 0
+        self.ID[0, :, :, 1] = 0
+        self.ID[1, :, :, 0] = 0
+        self.ID[1, :, :, 1] = 0
+
+    def reset_fields(self):
+        """Clear arrays for field components and PMLs."""
+        # Clear arrays for field components
+        self.initialise_field_arrays()
+
+        # Clear arrays for fields in PML
+        for pml in self.pmls:
+            pml.initialise_field_arrays()
+
+    def calculate_dt(self):
+        """Calculate time step at the CFL limit."""
+        self.dt = (1 / (c * np.sqrt(
+                  (1 / self.dx) * (1 / self.dx) +
+                  (1 / self.dy) * (1 / self.dy) +
+                  (1 / self.dz) * (1 / self.dz))))
+
+        # Round down time step to nearest float with precision one less than
+        # hardware maximum. Avoids inadvertently exceeding the CFL due to
+        # binary representation of floating point number.
+        self.dt = round_value(self.dt, decimalplaces=d.getcontext().prec - 1)
+
+
+class CUDAGrid(FDTDGrid):
+    """Additional grid methods for solving on GPU using CUDA."""
+
+    def set_blocks_per_grid(self):
         """Set the blocks per grid size used for updating the electric and magnetic field arrays on a GPU."""
+
         config.cuda['gpus'].bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / config.cuda['gpus'].tpb[0])), 1, 1)
 
-    def gpu_initialise_arrays(self):
-        """Initialise standard field arrays on GPU."""
+    def initialise_arrays(self):
+        """Initialise geometry and field arrays on GPU."""
 
         import pycuda.gpuarray as gpuarray
 
@@ -264,7 +286,7 @@ class FDTDGrid(Grid):
         self.Hy_gpu = gpuarray.to_gpu(self.Hy)
         self.Hz_gpu = gpuarray.to_gpu(self.Hz)
 
-    def gpu_initialise_dispersive_arrays(self):
+    def initialise_dispersive_arrays(self):
         """Initialise dispersive material coefficient arrays on GPU."""
 
         import pycuda.gpuarray as gpuarray
@@ -274,46 +296,19 @@ class FDTDGrid(Grid):
         self.Tz_gpu = gpuarray.to_gpu(self.Tz)
         self.updatecoeffsdispersive_gpu = gpuarray.to_gpu(self.updatecoeffsdispersive)
 
-    # Add PEC boundaries to invariant direction in 2D modes
-    # N.B. 2D modes are a single cell slice of 3D grid
-    def tmx(self):
-        # Ey & Ez components
-        self.ID[1, 0, :, :] = 0
-        self.ID[1, 1, :, :] = 0
-        self.ID[2, 0, :, :] = 0
-        self.ID[2, 1, :, :] = 0
+    def memory_check(self, snapsmemsize=0):
+        """Check if model can be run on specified GPU."""
 
-    def tmy(self):
-        # Ex & Ez components
-        self.ID[0, :, 0, :] = 0
-        self.ID[0, :, 1, :] = 0
-        self.ID[2, :, 0, :] = 0
-        self.ID[2, :, 1, :] = 0
+        super().memory_check()
 
-    def tmz(self):
-        # Ex & Ey components
-        self.ID[0, :, :, 0] = 0
-        self.ID[0, :, :, 1] = 0
-        self.ID[1, :, :, 0] = 0
-        self.ID[1, :, :, 1] = 0
+        if config.cuda['gpus'] is not None:
+            if self.memoryusage - snapsmemsize > config.cuda['gpus'].totalmem:
+                raise GeneralError(f'Memory (RAM) required ~{human_size(self.memoryusage)} exceeds {human_size(config.cuda['gpus'].totalmem, a_kilobyte_is_1024_bytes=True} detected on specified {config.cuda['gpus'].deviceID} - {config.cuda['gpus'].name} GPU!\n')
 
-
-    def reset_fields(self):
-        # Clear arrays for field components
-        self.initialise_field_arrays()
-
-        # Clear arrays for fields in PML
-        for pml in self.pmls:
-            pml.initialise_field_arrays()
-
-    def calculate_dt(self):
-        self.dt = 1 / (c * np.sqrt((1 / self.dx) * (1 / self.dx) + (1 / self.dy) * (1 / self.dy) + (1 / self.dz) * (1 / self.dz)))
-
-    def round_time_step(self):
-        # Round down time step to nearest float with precision one less than hardware maximum. Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
-        # Round down time step to nearest float with precision one less than hardware maximum.
-        # Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
-        self.dt = round_value(self.dt, decimalplaces=d.getcontext().prec - 1)
+            # If the required memory for the model without the snapshots will
+            # fit on the GPU then transfer and store snaphots on host
+            if snapsmemsize != 0 and self.memoryusage - snapsmemsize < config.cuda['gpus'].totalmem:
+                config.cuda['snapsgpu2cpu'] = True
 
 
 def dispersion_analysis(G):
@@ -487,7 +482,3 @@ def Iz(x, y, z, Hx, Hy, Hz, G):
         Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z])
 
     return Iz
-
-
-class GPUGrid(FDTDGrid):
-    pass
diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py
index e2e32ed4..89bcef9a 100644
--- a/gprMax/input_cmds_file.py
+++ b/gprMax/input_cmds_file.py
@@ -308,7 +308,7 @@ def parse_hash_commands(model_config, G, scene):
 
         # Write a file containing the input commands after Python or include
         # file commands have been processed
-        if sim_config.write_processed:
+        if sim_config.args.write_processed:
             write_processed_file(processedlines, model_config.appendmodelnumber, G)
 
         user_objs = get_user_objects(processedlines, check=True)
diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py
index 8cc65f23..7545af3a 100644
--- a/gprMax/model_build_run.py
+++ b/gprMax/model_build_run.py
@@ -19,6 +19,7 @@
 import datetime
 from importlib import import_module
 import itertools
+import logging
 import os
 import psutil
 import sys
@@ -74,6 +75,8 @@ from .utilities import timer
 from .utilities import Printer
 
 
+log = logging.getLogger(__name__)
+
 class ModelBuildRun:
     """Builds and runs (solves) a model."""
 
@@ -159,7 +162,7 @@ class ModelBuildRun:
             gb.grid.initialise_std_update_coeff_arrays()
 
         # Set datatype for dispersive arrays if there are any dispersive materials.
-        if config.materials['maxpoles'] != 0:
+        if self.model_config.materials['maxpoles'] != 0:
             drudelorentz = any([m for m in G.materials if 'drude' in m.type or 'lorentz' in m.type])
             if drudelorentz:
                 config.materials['dispersivedtype'] = config.dtypes['complex']
@@ -206,8 +209,8 @@ class ModelBuildRun:
         G = self.G
         # Reset iteration number
         G.iteration = 0
-        self.model_config.inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(self.model_config.appendmodelnumber, self.sim_config.model_end, self.sim_config.input_file_path)
-        self.printer.print(Fore.GREEN + '{} {}'.format(self.model_config.inputfilestr, '-' * (get_terminal_width() - 1 - len(self.model_config.inputfilestr))) + Style.RESET_ALL)
+        self.model_config.inputfilestr = f'\n--- Model {self.model_config.appendmodelnumber}/{self.sim_config.model_end}, input file (not re-processed, i.e. geometry fixed): {self.sim_config.input_file_path}'
+        log.info(Fore.GREEN + f'{self.model_config.inputfilestr} {'-' * (get_terminal_width() - 1 - len(self.model_config.inputfilestr))}' + Style.RESET_ALL)
         for grid in [G] + G.subgrids:
             grid.reset_fields()
 
@@ -232,52 +235,52 @@ class ModelBuildRun:
             # Check and set output directory and filename
             try:
                 os.mkdir(self.G.outputdirectory)
-                self.printer.print('\nCreated output directory: {}'.format(self.G.outputdirectory))
+                log.info(f'\nCreated output directory: {self.G.outputdirectory}')
             except FileExistsError:
                 pass
             # Modify the output path (hack)
             self.model_config.output_file_path_ext = Path(self.G.outputdirectory, self.model_config.output_file_path_ext)
 
-    def run_model(self, solver):
-        G = self.G
-        self.create_output_directory()
-        self.printer.print('Output file: {}\n'.format(self.model_config.output_file_path_ext))
+    def write_output_data(self):
+        """Write output data, i.e. field data for receivers and snapshots
+            to file(s).
+        """
 
-        # Run solver
-        tsolve = self.solve(solver)
+        G = self.G
 
         # Write an output file in HDF5 format
         write_hdf5_outputfile(self.model_config.output_file_path_ext, G)
 
         # Write any snapshots to file
         if G.snapshots:
-            # Create directory and construct filename from user-supplied name and model run number
+            # Create directory and construct filename from user-supplied name
+            # and model run number
             snapshotdir = self.model_config.snapshot_dir
             if not os.path.exists(snapshotdir):
                 os.mkdir(snapshotdir)
 
-            self.printer.print('')
+            log.info('')
             for i, snap in enumerate(G.snapshots):
                 fn = snapshotdir / Path(self.model_config.output_file_path.stem + '_' + snap.basefilename)
                 snap.filename = fn.with_suffix('.vti')
                 pbar = tqdm(total=snap.vtkdatawritesize, leave=True, unit='byte', unit_scale=True, desc='Writing snapshot file {} of {}, {}'.format(i + 1, len(G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
                 snap.write_vtk_imagedata(pbar, G)
                 pbar.close()
-            self.printer.print('')
+            log.info('')
+
+    def print_resource_info(self):
+        """Print resource information on runtime and memory usage."""
 
         memGPU = ''
         if config.cuda['gpus']:
-            memGPU = ' host + ~{} GPU'.format(human_size(self.solver.get_memsolve()))
+            memGPU = f' host + ~{human_size(self.solver.get_memsolve())} GPU'
 
-        self.printer.print('\nMemory (RAM) used: ~{}{}'.format(human_size(self.p.memory_full_info().uss), memGPU))
-        self.printer.print('Solving time [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=tsolve)))
-
-        return tsolve
+        log.info(f'\nMemory (RAM) used: ~{human_size(self.p.memory_full_info().uss)}{memGPU}')
+        log.info(f'Solving time [HH:MM:SS]: {datetime.timedelta(seconds=tsolve)}')
 
     def solve(self, solver):
         """
-        Solving using FDTD method on CPU. Parallelised using Cython (OpenMP) for
-        electric and magnetic field updates, and PML updates.
+        Solve using FDTD method.
 
         Args:
             solver (Solver): solver object.
@@ -294,7 +297,17 @@ class ModelBuildRun:
         else:
             iterator = range(0, G.iterations)
 
-        tsolve = solver.solve(iterator)
+        self.create_output_directory()
+        log.info(f'Output file: {self.model_config.output_file_path_ext}\n')
+
+        # Run solver
+        tsolve = self.solve(solver)
+
+        # Write output data, i.e. field data for receivers and snapshots to file(s)
+        self.write_output_data
+
+        # Print resource information on runtime and memory usage
+        self.print_resource_info
 
         return tsolve
 
@@ -318,7 +331,7 @@ class GridBuilder:
     def build_components(self):
         # Build the model, i.e. set the material properties (ID) for every edge
         # of every Yee cell
-        self.printer.print('')
+        log.info('')
         pbar = tqdm(total=2, desc='Building {} Grid'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
         build_electric_components(self.grid.solid, self.grid.rigidE, self.grid.ID, self.grid)
         pbar.update()
@@ -348,5 +361,5 @@ class GridBuilder:
         materialstable.outer_border = False
         materialstable.justify_columns[0] = 'right'
 
-        self.printer.print('\n{} Grid Materials:'.format(self.grid.name))
-        self.printer.print(materialstable.table)
+        log.info(f'\n{self.grid.name} Grid Materials:')
+        log.info(materialstable.table)
diff --git a/gprMax/pml.py b/gprMax/pml.py
index 4484db73..7370a759 100644
--- a/gprMax/pml.py
+++ b/gprMax/pml.py
@@ -81,7 +81,7 @@ class CFS(object):
 
         # Calculation of the maximum value of sigma from http://dx.doi.org/10.1109/8.546249
         m = CFSParameter.scalingprofiles[self.sigma.scalingprofile]
-        self.sigma.max = (0.8 * (m + 1)) / (config.z0 * d * np.sqrt(er * mr))
+        self.sigma.max = (0.8 * (m + 1)) / (config.general['z0'] * d * np.sqrt(er * mr))
 
     def scaling_polynomial(self, order, Evalues, Hvalues):
         """Applies the polynomial to be used for the scaling profile for
diff --git a/gprMax/solvers.py b/gprMax/solvers.py
index 44064072..4e58d7cf 100644
--- a/gprMax/solvers.py
+++ b/gprMax/solvers.py
@@ -18,10 +18,10 @@
 
 import gprMax.config as config
 from .grid import FDTDGrid
-from .grid import GPUGrid
+from .grid import CUDAGrid
 from .subgrids.updates import create_updates as create_subgrid_updates
 from .updates import CPUUpdates
-from .updates import GPUUpdates
+from .updates import CUDAUpdates
 
 
 def create_G(sim_config):
@@ -34,8 +34,8 @@ def create_G(sim_config):
         G (Grid): Holds essential parameters describing the model.
     """
 
-    if sim_config.gpu:
-        G = GPUGrid()
+    if sim_config.general['cuda']:
+        G = CUDAGrid()
     elif sim_config.subgrid:
         G = FDTDGrid()
     else:
@@ -56,7 +56,7 @@ def create_solver(G, sim_config):
     """
 
     if sim_config.gpu:
-        updates = GPUUpdates(G)
+        updates = CUDAUpdates(G)
         solver = Solver(updates)
     elif sim_config.subgrid:
         updates = create_subgrid_updates(G)
diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py
index 7a3bd666..86176d63 100644
--- a/gprMax/subgrids/base.py
+++ b/gprMax/subgrids/base.py
@@ -63,19 +63,9 @@ class SubGridBase(FDTDGrid):
 
         self.interpolation = kwargs['interpolation']
 
-    def calculate_dt(self):
-        """Calculate dt at the CFL limit."""
-        self.dt = (1 / (c * np.sqrt(
-                  (1 / self.dx) * (1 / self.dx) +
-                  (1 / self.dy) * (1 / self.dy) +
-                  (1 / self.dz) * (1 / self.dz))))
-
     def main_grid_index_to_subgrid_index(self, i, j, k):
         """Calculate local subgrid index from global main grid index."""
         i_s = self.n_boundary_cells_x + (i - self.i0) * self.ratio
         j_s = self.n_boundary_cells_y + (j - self.j0) * self.ratio
         k_s = self.n_boundary_cells_z + (k - self.k0) * self.ratio
         return (i_s, j_s, k_s)
-
-    def initialise_geometry_arrays(self):
-        super().initialise_geometry_arrays()
diff --git a/gprMax/subgrids/multi.py b/gprMax/subgrids/multi.py
index 9d92825e..0e95b59b 100644
--- a/gprMax/subgrids/multi.py
+++ b/gprMax/subgrids/multi.py
@@ -27,10 +27,6 @@ class ReferenceRx(Rx):
         position as the coarse grid.
     """
 
-    def __init__(self):
-        """Constructor."""
-        super().__init__()
-
     def get_field(self, str_id, field):
         """Return the field value at the equivalent coarse yee cell.
 
diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py
index a14a06b4..a3e83ac1 100644
--- a/gprMax/subgrids/precursor_nodes.py
+++ b/gprMax/subgrids/precursor_nodes.py
@@ -28,7 +28,7 @@ def calculate_weighting_coefficients(x1, x):
     return (c1, c2)
 
 
-class PrecusorNodesBase(object):
+class PrecusorNodesBase:
 
     def __init__(self, fdtd_grid, sub_grid):
         self.G = fdtd_grid
diff --git a/gprMax/templates/fields_updates_dispersive_template b/gprMax/templates/fields_updates_dispersive_template
index 9906ced5..05363ffd 100644
--- a/gprMax/templates/fields_updates_dispersive_template
+++ b/gprMax/templates/fields_updates_dispersive_template
@@ -18,6 +18,7 @@
 
 import numpy as np
 cimport numpy as np
+
 from cython.parallel import prange
 
 cdef extern from "complex.h" nogil:
diff --git a/gprMax/updates.py b/gprMax/updates.py
index bf7cb0d7..a30ac8a0 100644
--- a/gprMax/updates.py
+++ b/gprMax/updates.py
@@ -17,6 +17,8 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 from importlib import import_module
+import logging
+import sys
 
 import gprMax.config as config
 from .cuda.fields_updates import kernel_template_fields
@@ -34,6 +36,8 @@ from .sources import gpu_initialise_src_arrays
 from .utilities import timer
 
 
+log = logging.getLogger(__name__)
+
 class CPUUpdates:
     """Defines update functions for CPU-based solver."""
 
@@ -226,7 +230,7 @@ class CPUUpdates:
         pass
 
 
-class GPUUpdates:
+class CUDAUpdates:
     """Defines update functions for GPU-based (CUDA) solver."""
 
     def __init__(self, G):
@@ -235,6 +239,10 @@ class GPUUpdates:
             G (FDTDGrid): FDTD grid object
         """
 
+        self.grid = G
+        self.dispersive_update_a = None
+        self.dispersive_update_b = None
+
         import pycuda.driver as drv
         from pycuda.compiler import SourceModule
         drv.init()
@@ -249,10 +257,6 @@ class GPUUpdates:
         self.dev = drv.Device(self.grid.gpu.deviceID)
         self.ctx = dev.make_context()
 
-        self.grid = G
-        self.dispersive_update_a = None
-        self.dispersive_update_b = None
-
         # Initialise arrays on GPU, prepare kernels, and get kernel functions
         self.set_field_kernels()
         self.set_pml_kernels()
@@ -261,7 +265,9 @@ class GPUUpdates:
         self.set_snapshot_kernel()
 
     def set_field_kernels(self):
-        # Electric and magnetic field updates - prepare kernels, and get kernel functions
+        """Electric and magnetic field updates - prepare kernels, and
+            get kernel functions.
+        """
         if config.materials['maxpoles'] > 0:
             kernels_fields = SourceModule(kernels_template_fields.substitute(REAL=cudafloattype, COMPLEX=cudacomplextype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=self.grid.updatecoeffsdispersive.shape[1], NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3], NX_T=self.grid.Tx.shape[1], NY_T=self.grid.Tx.shape[2], NZ_T=self.grid.Tx.shape[3]), options=self.compiler_opts)
         else: # Set to one any substitutions for dispersive materials
@@ -269,7 +275,7 @@ class GPUUpdates:
         self.update_electric = kernels_fields.get_function("update_electric")
         self.update_magnetic = kernels_fields.get_function("update_magnetic")
         if self.grid.updatecoeffsE.nbytes + self.grid.updatecoeffsH.nbytes > self.grid.gpu.constmem:
-            raise GeneralError('Too many materials in the model to fit onto constant memory of size {} on {} - {} GPU'.format(human_size(self.grid.gpu.constmem), self.grid.gpu.deviceID, self.grid.gpu.name))
+            raise GeneralError(log.exception(f'Too many materials in the model to fit onto constant memory of size {human_size(self.grid.gpu.constmem)} on {self.grid.gpu.deviceID} - {self.grid.gpu.name} GPU'))
         self.copy_mat_coeffs()
 
         # Electric and magnetic field updates - dispersive materials - get kernel functions and initialise array on GPU
@@ -283,8 +289,8 @@ class GPUUpdates:
         self.grid.gpu_initialise_arrays()
 
     def set_pml_kernels(self):
+        """PMLS - prepare kernels and get kernel functions."""
         if self.grid.pmls:
-            # PMLS - prepare kernels and get kernel functions
             pmlmodulelectric = 'gprMax.cuda.pml_updates_electric_' + self.grid.pmlformulation
             kernelelectricfunc = getattr(import_module(pmlmodulelectric), 'kernels_template_pml_electric_' + self.grid.pmlformulation)
             pmlmodulemagnetic = 'gprMax.cuda.pml_updates_magnetic_' + self.grid.pmlformulation
@@ -299,14 +305,18 @@ class GPUUpdates:
                 pml.gpu_set_blocks_per_grid(self.grid)
 
     def set_rx_kernel(self):
-        # Receivers - initialise arrays on GPU, prepare kernel and get kernel function
+        """Receivers - initialise arrays on GPU, prepare kernel and
+            get kernel function.
+        """
         if self.grid.rxs:
             rxcoords_gpu, rxs_gpu = gpu_initialise_rx_arrays(self.grid)
             kernel_store_outputs = SourceModule(kernel_template_store_outputs.substitute(REAL=cudafloattype, NY_RXCOORDS=3, NX_RXS=6, NY_RXS=self.grid.iterations, NZ_RXS=len(self.grid.rxs), NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1), options=self.compiler_opts)
             self.store_outputs = kernel_store_outputs.get_function("store_outputs")
 
     def set_src_kernels(self):
-        # Sources - initialise arrays on GPU, prepare kernel and get kernel function
+        """Sources - initialise arrays on GPU, prepare kernel and
+            get kernel function.
+        """
         if self.grid.voltagesources + self.grid.hertziandipoles + self.grid.magneticdipoles:
             kernels_sources = SourceModule(kernels_template_sources.substitute(REAL=cudafloattype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_SRCINFO=4, NY_SRCWAVES=self.grid.iterations, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3]), options=self.compiler_opts)
             self.copy_mat_coeffs()
@@ -321,14 +331,18 @@ class GPUUpdates:
                 self.update_voltage_source_gpu = kernels_sources.get_function("update_voltage_source")
 
     def set_snapshot_kernel(self):
+        """Snapshots - initialise arrays on GPU, prepare kernel and
+            get kernel function.
+        """
         if self.grid.snapshots:
-            # Snapshots - initialise arrays on GPU, prepare kernel and get kernel function
             self.snapEx_gpu, self.snapEy_gpu, self.snapEz_gpu, self.snapHx_gpu, self.snapHy_gpu, self.snapHz_gpu = gpu_initialise_snapshot_array(self.grid)
             kernel_store_snapshot = SourceModule(kernel_template_store_snapshot.substitute(REAL=cudafloattype, NX_SNAPS=Snapshot.nx_max, NY_SNAPS=Snapshot.ny_max, NZ_SNAPS=Snapshot.nz_max, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1), options=self.compiler_opts)
             self.store_snapshot_gpu = kernel_store_snapshot.get_function("store_snapshot")
 
     def copy_mat_coeffs(self):
-        # Copy material coefficient arrays to constant memory of GPU (must be <64KB)
+        """Copy material coefficient arrays to constant memory of GPU
+            (must be <64KB).
+        """
         updatecoeffsE = kernels_sources.get_global('updatecoeffsE')[0]
         updatecoeffsH = kernels_sources.get_global('updatecoeffsH')[0]
         self.drv.memcpy_htod(updatecoeffsE, self.grid.updatecoeffsE)
@@ -396,7 +410,6 @@ class GPUUpdates:
 
     def update_magnetic(self):
         """Update magnetic field components."""
-
         self.update_magnetic(np.int32(self.grid.nx),
                         np.int32(self.grid.ny),
                         np.int32(self.grid.nz),
@@ -412,13 +425,11 @@ class GPUUpdates:
 
     def update_magnetic_pml(self):
         """Update magnetic field components with the PML correction."""
-
         for pml in self.grid.pmls:
             pml.gpu_update_magnetic(self.grid)
 
     def update_magnetic_sources(self):
         """Update magnetic field components from sources."""
-
         if self.grid.magneticdipoles:
             self.update_magnetic_dipole_gpu(np.int32(len(self.grid.magneticdipoles)),
                                               np.int32(self.grid.iteration),
@@ -437,7 +448,6 @@ class GPUUpdates:
 
     def update_electric_a(self):
         """Update electric field components."""
-
         # All materials are non-dispersive so do standard update.
         if config.materials['maxpoles'] == 0:
             self.update_electric(np.int32(self.grid.nx),
@@ -476,7 +486,6 @@ class GPUUpdates:
 
     def update_electric_pml(self):
         """Update electric field components with the PML correction."""
-
         for pml in self.grid.pmls:
             pml.gpu_update_electric(self.grid)
 
@@ -484,7 +493,6 @@ class GPUUpdates:
         """Update electric field components from sources -
             update any Hertzian dipole sources last.
         """
-
         if self.grid.voltagesources:
             self.update_voltage_source_gpu(np.int32(len(self.grid.voltagesources)),
                                            np.int32(self.grid.iteration),
@@ -525,7 +533,6 @@ class GPUUpdates:
             field values. Therefore it can only be completely updated after the
             electric field has been updated by the PML and source updates.
         """
-
         if config.materials['maxpoles'] != 0:
             self.dispersive_update_b(np.int32(self.grid.nx),
                                      np.int32(self.grid.ny),
@@ -543,11 +550,14 @@ class GPUUpdates:
                                      grid=self.grid.bpg)
 
     def time_start(self):
+        """Start event timers used to calculate solving time for model."""
         self.iterstart = self.drv.Event()
         self.iterend = self.drv.Event()
         self.iterstart.record()
+        self.iterstart.synchronize()
 
     def calculate_tsolve(self):
+        """Calculate solving time for model."""
         self.iterend.record()
         self.iterend.synchronize()
         tsolve = self.iterstart.time_till(self.iterend) * 1e-3
@@ -555,6 +565,7 @@ class GPUUpdates:
         return tsolve
 
     def finalise(self):
+        """Copy data from GPU back to CPU to save to file(s)."""
         # Copy output from receivers array back to correct receiver objects
         if self.grid.rxs:
             gpu_get_rx_array(self.rxs_gpu.get(),
@@ -574,6 +585,7 @@ class GPUUpdates:
                                        snap)
 
     def cleanup(self):
+        """Cleanup GPU context."""
         # Remove context from top of stack and delete
         self.ctx.pop()
         del self.ctx
diff --git a/gprMax/utilities.py b/gprMax/utilities.py
index 71630ee0..6f2fdfb5 100644
--- a/gprMax/utilities.py
+++ b/gprMax/utilities.py
@@ -27,7 +27,7 @@ import subprocess
 from shutil import get_terminal_size
 import sys
 import textwrap
-from time import perf_counter
+from time import process_time
 
 from colorama import init
 from colorama import Fore
@@ -38,7 +38,6 @@ import numpy as np
 from .exceptions import GeneralError
 
 
-
 def get_terminal_width():
     """Get/set width of terminal being used.
 
@@ -75,14 +74,14 @@ def logo(version):
     |___/|_|
                      v""" + version
 
-    print('{} {}\n'.format(description, '=' * (get_terminal_width() - len(description) - 1)))
-    print(Fore.CYAN + '{}\n'.format(logo))
-    print(Style.RESET_ALL + textwrap.fill(copyright, width=get_terminal_width() - 1, initial_indent=' '))
-    print(textwrap.fill(authors, width=get_terminal_width() - 1, initial_indent=' '))
-    print()
-    print(textwrap.fill(licenseinfo1, width=get_terminal_width() - 1, initial_indent=' ', subsequent_indent='  '))
-    print(textwrap.fill(licenseinfo2, width=get_terminal_width() - 1, initial_indent=' ', subsequent_indent='  '))
-    print(textwrap.fill(licenseinfo3, width=get_terminal_width() - 1, initial_indent=' ', subsequent_indent='  '))
+    log.info(f'{description} {'=' * (get_terminal_width() - len(description) - 1)}\n')
+    log.info(Fore.CYAN + f'{logo}\n')
+    log.info(Style.RESET_ALL + textwrap.fill(copyright, width=get_terminal_width() - 1, initial_indent=' '))
+    log.info(textwrap.fill(authors, width=get_terminal_width() - 1, initial_indent=' '))
+    log.info('')
+    log.info(textwrap.fill(licenseinfo1, width=get_terminal_width() - 1, initial_indent=' ', subsequent_indent='  '))
+    log.info(textwrap.fill(licenseinfo2, width=get_terminal_width() - 1, initial_indent=' ', subsequent_indent='  '))
+    log.info(textwrap.fill(licenseinfo3, width=get_terminal_width() - 1, initial_indent=' ', subsequent_indent='  '))
 
 
 @contextmanager
@@ -401,7 +400,7 @@ def detect_check_gpus(deviceIDs):
     # Check if requested device ID(s) exist
     for ID in deviceIDs:
         if ID not in deviceIDsavail:
-            raise GeneralError('GPU with device ID {} does not exist'.format(ID))
+            raise GeneralError(f'GPU with device ID {ID} does not exist')
 
     # Gather information about selected/detected GPUs
     gpus = []
@@ -411,22 +410,11 @@ def detect_check_gpus(deviceIDs):
         gpu.get_gpu_info(drv)
         if ID in deviceIDs:
             gpus.append(gpu)
-        allgpustext.append('{} - {}, {}'.format(gpu.deviceID, gpu.name, human_size(gpu.totalmem, a_kilobyte_is_1024_bytes=True)))
+        allgpustext.append(f'{gpu.deviceID} - {gpu.name}, {human_size(gpu.totalmem, a_kilobyte_is_1024_bytes=True)}')
 
     return gpus, allgpustext
 
 
 def timer():
     """Function to return the current process wide time in fractional seconds."""
-    return perf_counter()
-
-
-class Printer():
-    """Printing information messages."""
-    
-    def __init__(self, config):
-        self.printing = config.is_messages()
-
-    def print(self, str):
-        if self.printing:
-            print(str)
+    return process_time()
diff --git a/setup.py b/setup.py
index 34dcc9b0..0f8a2217 100644
--- a/setup.py
+++ b/setup.py
@@ -33,14 +33,18 @@ import pathlib
 import re
 import shutil
 import sys
-from jinja2 import Environment, PackageLoader, select_autoescape
+
+from jinja2 import Environment
+from jinja2 import PackageLoader
+from jinja2 import select_autoescape
 
 
 def build_dispersive_material_templates():
+    """Function to generate Cython .pyx files for dispersive media update.
+        Jinja2 templates are used to render the various dispersive update
+        functions.
     """
-    Function to generate Cython .pyx files for dispersive media update.
-    Jinja2 templates are used to render the various dispersive update functions.
-    """
+
     env = Environment(
         loader=PackageLoader(__name__, 'gprMax/templates'),
     )
@@ -152,25 +156,25 @@ if 'cleanall' in sys.argv:
         if os.path.isfile(filebase + '.c'):
             try:
                 os.remove(filebase + '.c')
-                print('Removed: {}'.format(filebase + '.c'))
+                print(f'Removed: {filebase} + ".c"')
             except OSError:
-                print('Could not remove: {}'.format(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') + glob.glob(os.path.join(os.getcwd(), os.path.splitext(file)[0]) + '*.so')
         if libfile:
             libfile = libfile[0]
             try:
                 os.remove(libfile)
-                print('Removed: {}'.format(os.path.abspath(libfile)))
+                print(f'Removed: {os.path.abspath(libfile)}')
             except OSError:
-                print('Could not remove: {}'.format(os.path.abspath(libfile)))
+                print(f'Could not remove: {os.path.abspath(libfile)}')
     # Remove build, dist, egg and __pycache__ directories
     shutil.rmtree(os.path.join(os.getcwd(), 'build'), ignore_errors=True)
     shutil.rmtree(os.path.join(os.getcwd(), 'dist'), ignore_errors=True)
     shutil.rmtree(os.path.join(os.getcwd(), 'gprMax.egg-info'), ignore_errors=True)
     for p in pathlib.Path(os.getcwd()).rglob('__pycache__'):
         shutil.rmtree(p, ignore_errors=True)
-        print('Removed: {}'.format(p))
+        print(f'Removed: {p}')
     # Now do a normal clean
     sys.argv[1] = 'clean'  # this is what distutils understands