Merge branch 'devel' of https://github.com/gprMax/gprMax into f-strings

2025-08-07 15:10:13 +08:00 · 2023-09-08 15:58:24 +05:30
--- a/README.rst
+++ b/README.rst
@@ -117,7 +117,7 @@ Microsoft Windows
 ^^^^^^^^^^^^^^^^^

 * Download and install Microsoft `Build Tools for Visual Studio 2022 <https://aka.ms/vs/17/release/vs_BuildTools.exe>`_ (direct link). You can also find it on the `Microsoft Visual Studio downloads page <https://visualstudio.microsoft.com/downloads/>`_ by scrolling down to the 'All Downloads' section, clicking the disclosure triangle by 'Tools for Visual Studio 2022', then clicking the download button next to 'Build Tools for Visual Studio 2022'. When installing, choose the 'Desktop development with C++' Workload and select only 'MSVC v143' and 'Windows 10 SDK' or 'Windows 11 SDK options.
-* Set the Path and Environment Variables - this can be done by following the `instructions from Microsoft <https://docs.microsoft.com/en-us/cpp/build/building-on-the-command-line?view=msvc-160#developer_command_file_locations>`_, or manually by adding a form of :code:`C:\Program Files (x86)\Microsoft Visual Studio\2019\BuildTools\VC\Tools\MSVC\14.23.28105\bin\Hostx64\x64` (this may vary according to your exact machine and installation) to your system Path environment variable.
+* Set the Path and Environment Variables - this can be done by following the `instructions from Microsoft <https://docs.microsoft.com/en-us/cpp/build/building-on-the-command-line?view=msvc-160#developer_command_file_locations>`_, or manually by adding a form of :code:``C:\Program Files (x86)\Microsoft Visual Studio\2019\BuildTools\VC\Tools\MSVC\14.23.28105\bin\Hostx64\x64`` (this may vary according to your exact machine and installation) to your system Path environment variable.

 Alternatively, if you are using Windows 10/11 you can install the `Windows Subsystem for Linux <https://docs.microsoft.com/en-gb/windows/wsl/about>`_ and then follow the Linux install instructions for gprMax. Note however that currently, WSL does not aim to support GUI desktops or applications, e.g. Gnome, KDE, etc....

--- a/docs/source/examples_simple_2D.rst
+++ b/docs/source/examples_simple_2D.rst
@@ -112,7 +112,7 @@ You should have produced an output file ``cylinder_Ascan_2D.h5``. You can view t

 .. code-block:: none

-    python -m tools.plot_Ascan examples/cylinder_Ascan_2D.h5
+    python -m toolboxes.Plotting.plot_Ascan examples/cylinder_Ascan_2D.h5

 :numref:`cylinder_Ascan_results` shows the time history of the electric and magnetic field components and currents at the receiver location. The :math:`E_z` field component can be converted to a voltage that represents the A-scan (trace). The initial part of the signal (~0.5-1.5 ns) represents the direct wave from transmitter to receiver. Then comes the reflected wavelet (~1.8-2.6 ns), which has opposite polarity, from the metal cylinder.

@@ -153,7 +153,7 @@ You should have produced 60 output files, one for each A-scan, with names ``cyli

 .. code-block:: none

-    python -m tools.outputfiles_merge examples/cylinder_Bscan_2D
+    python -m toolboxes.Utilities.outputfiles_merge examples/cylinder_Bscan_2D

 You should see a combined output file ``cylinder_Bscan_2D_merged.h5``. You can add the optional argument ``--remove-files`` if you want to automatically delete the original single A-scan output files.

@@ -161,7 +161,7 @@ You can now view an image of the B-scan using the command:

 .. code-block:: none

-    python -m tools.plot_Bscan examples/cylinder_Bscan_2D_merged.h5 Ez
+    python -m toolboxes.Plotting.plot_Bscan examples/cylinder_Bscan_2D_merged.h5 Ez

 :numref:`cylinder_Bscan_results` shows the B-scan (of the :math:`E_z` field component). Again, the initial part of the signal (~0.5-1.5 ns) represents the direct wave from transmitter to receiver. Then comes the reflected wave (~2-3 ns) from the metal cylinder which creates the hyperbolic shape.

--- a/docs/source/input_hash_cmds.rst
+++ b/docs/source/input_hash_cmds.rst
@@ -965,7 +965,7 @@ For example to save a snapshot of the electromagnetic fields in the model at a s
 PML commands
 ============

-The default behaviour for the absorbing boundary conditions (ABC) is first order Complex Frequency Shifted (CFS) Perfectly Matched Layers (PML), with thicknesses of 10 cells on each of the six sides of the model domain. The thickness of the PML can be altered by using the following command (further customisation of the PML is possible using our Python API):
+The default behaviour for the absorbing boundary conditions (ABC) is first order Complex Frequency Shifted (CFS) Perfectly Matched Layers (PML), with thicknesses of 10 cells on each of the six sides of the model domain.

 #pml_cells:
 ------------
@@ -989,3 +989,47 @@ For example to use a PML with 20 cells (thicker than the default 10 cells) on on
 .. code-block:: none

    #pml_cells: 10 10 20 10 10 20
+
+#pml_formulation:
+-----------------
+
+Allows you to alter the formulation used for the PML. The current options are to use the Higher Order RIPML (HORIPML) - https://doi.org/10.1109/TAP.2011.2180344, or Multipole RIPML (MRIPML) - https://doi.org/10.1109/TAP.2018.2823864. The syntax of the command is:
+
+.. code-block:: none
+
+    #pml_formulation: str
+
+* ``str`` can be either 'HORIPML' or 'MRIPML'
+
+For example to use the Multipole RIPML:
+
+.. code-block:: none
+
+    #pml_formulation: MRIPML
+
+#pml_cfs:
+---------
+
+Allows you (advanced) control of the parameters that are used to build each order of the PML. Up to a second order PML can currently be specified, i.e. by using two ``#pml_cfs`` commands. The syntax of the command is:
+
+.. code-block:: none
+
+    #pml_cfs: str1 str2 f1 f2 str3 str4 f3 f4 str5 str6 f5 f6
+
+* ``str1`` is the type of scaling to use for the CFS :math:`\alpha` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``.
+* ``str2`` is the direction of the scaling to use for the CFS :math:`\alpha` parameter. It can be ``forward`` or ``reverse``.
+* ``f1 f2`` are the minimum and maximum values for the CFS :math:`\alpha` parameter.
+* ``str3`` is the type of scaling to use for the CFS :math:`\kappa` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``.
+* ``str4`` is the direction of the scaling to use for the CFS :math:`\kappa` parameter. It can be ``forward`` or ``reverse``.
+* ``f3 f4`` are the minimum and maximum values for the CFS :math:`\kappa` parameter. The minimum value for the CFS :math:`\kappa` parameter is one.
+* ``str5`` is the type of scaling to use for the CFS :math:`\sigma` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``.
+* ``str6`` is the direction of the scaling to use for the CFS :math:`\sigma` parameter. It can be ``forward`` or ``reverse``.
+* ``f5 f6`` are the minimum and maximum values for the CFS :math:`\sigma` parameter.
+
+The CFS values (which are internally specified) used for the default standard first order PML are: ``#pml_cfs: constant forward 0 0 constant forward 1 1 quartic forward 0 None``. Specifying 'None' for the maximum value of :math:`\sigma` forces gprMax to calculate it internally based on the relative permittivity and permeability of the underlying materials in the model.
+
+The parameters will be applied to all slabs of the PML that are switched on.
+
+.. tip::
+    
+    ``forward`` direction implies minimum parameter value at the inner boundary of the PML and maximum parameter value at the edge of computational domain, ``reverse`` is the opposite.
--- a/examples/antenna_like_GSSI_400_fs.py
+++ b/examples/antenna_like_GSSI_400_fs.py
@@ -30,7 +30,7 @@ scene = gprMax.Scene()
 title = gprMax.Title(name=fn.with_suffix("").name)
 domain = gprMax.Domain(p1=(x, y, z))
 dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
-time_window = gprMax.TimeWindow(time=6e-9)
+time_window = gprMax.TimeWindow(time=15e-9)

 scene.add(title)
 scene.add(domain)
@@ -44,17 +44,20 @@ for obj in gssi_objects:
    scene.add(obj)

 gv1 = gprMax.GeometryView(
-    p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl), filename="antenna_like_GSSI_400", output_type="n"
+    p1=(0, 0, 0), 
+    p2=(x, y, z), 
+    dl=(dl, dl, dl), 
+    filename="antenna_like_GSSI_400", output_type="n"
 )
 gv2 = gprMax.GeometryView(
-    p1=(ant_pos[0] - 0.150 / 2, ant_pos[1] - 0.150 / 2, ant_pos[2] - 0.050),
-    p2=(ant_pos[0] + 0.150 / 2, ant_pos[1] + 0.150 / 2, ant_pos[2] + 0.010),
+    p1=(ant_pos[0] - 0.150, ant_pos[1] - 0.150, ant_pos[2]),
+    p2=(ant_pos[0] + 0.150, ant_pos[1] + 0.150, ant_pos[2] + 0.050),
    dl=(dl, dl, dl),
    filename="antenna_like_GSSI_400_pcb",
    output_type="f",
 )
-scene.add(gv1)
-scene.add(gv2)
+#scene.add(gv1)
+#scene.add(gv2)

 # Run model
 gprMax.run(scenes=[scene], geometry_only=False, outputfile=fn, gpu=None)
--- a/gprMax/init.py
+++ b/gprMax/init.py
@@ -14,9 +14,11 @@ from .cmds_geometry.add_grass import AddGrass
 from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness
 from .cmds_geometry.add_surface_water import AddSurfaceWater
 from .cmds_geometry.box import Box
+from .cmds_geometry.cone import Cone
 from .cmds_geometry.cylinder import Cylinder
 from .cmds_geometry.cylindrical_sector import CylindricalSector
 from .cmds_geometry.edge import Edge
+from .cmds_geometry.ellipsoid import Ellipsoid
 from .cmds_geometry.fractal_box import FractalBox
 from .cmds_geometry.geometry_objects_read import GeometryObjectsRead
 from .cmds_geometry.plate import Plate
@@ -38,6 +40,7 @@ from .cmds_multiuse import (
    RxArray,
    Snapshot,
    SoilPeplinski,
+    Subgrid,
    TransmissionLine,
    VoltageSource,
    Waveform,
@@ -47,6 +50,7 @@ from .cmds_singleuse import (
    Domain,
    ExcitationFile,
    OMPThreads,
+    OutputDir,
    PMLProps,
    RxSteps,
    SrcSteps,
--- a/gprMax/cmds_multiuse.py
+++ b/gprMax/cmds_multiuse.py
@@ -1737,11 +1737,6 @@ class PMLCFS(UserObjectMulti):
        ):
            logger.exception(f"{self.params_str()} minimum and maximum scaling values must be greater than zero.")
            raise ValueError
-        # TODO: Fix handling of kappa for 2nd order PMLs
-        # if float(kappamin) < 1:
-        #     logger.exception(f'{self.params_str()} minimum scaling value for '
-        #                      'kappa must be greater than or equal to one.')
-        #     raise ValueError

        cfsalpha = CFSParameter()
        cfsalpha.ID = "alpha"
@@ -1760,8 +1755,11 @@ class PMLCFS(UserObjectMulti):
        cfssigma.scalingprofile = sigmascalingprofile
        cfssigma.scalingdirection = sigmascalingdirection
        cfssigma.min = float(sigmamin)
+        if sigmamax == "None":
+            sigmamax = None
        if sigmamax is not None:
-            cfssigma.max = float(sigmamax)
+            sigmamax = float(sigmamax)
+        cfssigma.max = sigmamax
        cfs = CFS()
        cfs.alpha = cfsalpha
        cfs.kappa = cfskappa
@@ -1775,7 +1773,7 @@ class PMLCFS(UserObjectMulti):
            f"{cfskappa.scalingdirection}, min: {cfskappa.min:g}, max: "
            f"{cfskappa.max:g}), sigma (scaling: {cfssigma.scalingprofile}, "
            f"scaling direction: {cfssigma.scalingdirection}, min: "
-            f"{cfssigma.min:g}, max: {cfssigma.max:g}) created."
+            f"{cfssigma.min:g}, max: {cfssigma.max}) created."
        )

        grid.pmls["cfs"].append(cfs)
@@ -1800,4 +1798,4 @@ class Subgrid(UserObjectMulti):
            self.children_geometry.append(node)
        else:
            logger.exception("This object is unknown to gprMax.")
-            raise ValueError
+            raise ValueError
--- a/gprMax/cmds_singleuse.py
+++ b/gprMax/cmds_singleuse.py
@@ -433,7 +433,7 @@ class ExcitationFile(UserObjectSingle):
            waveformtime = np.arange(0, G.timewindow + G.dt, G.dt)
            timestr = "simulation time array"

-        for waveformID in waveformIDs:
+        for i, waveformID in enumerate(waveformIDs):
            if any(x.ID == waveformID for x in G.waveforms):
                logger.exception(f"Waveform with ID {waveformID} already exists")
                raise ValueError
@@ -442,7 +442,7 @@ class ExcitationFile(UserObjectSingle):
            w.type = "user"

            # Select correct column of waveform values depending on array shape
-            singlewaveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, waveform]
+            singlewaveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, i]

            # Truncate waveform array if it is longer than time array
            if len(singlewaveformvalues) > len(waveformtime):
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -209,33 +209,29 @@ class SimulationConfig:

        self.args = args

-        if args.mpi and args.geometry_fixed:
+        if self.args.mpi and self.args.geometry_fixed:
            logger.exception("The geometry fixed option cannot be used with MPI.")
            raise ValueError

-        if args.gpu and args.opencl:
+        if self.args.gpu and self.args.opencl:
            logger.exception("You cannot use both CUDA and OpenCl simultaneously.")
            raise ValueError

        # General settings for the simulation
-        #   inputfilepath: path to inputfile location.
-        #   outputfilepath: path to outputfile location.
-        #   progressbars: whether to show progress bars on stdoout or not.
        #   solver: cpu, cuda, opencl.
-        #   subgrid: whether the simulation uses sub-grids.
        #   precision: data type for electromagnetic field output (single/double).
        #   progressbars: progress bars on stdoout or not - switch off
        #                   progressbars when logging level is greater than
        #                   info (20)

-        self.general = {"solver": "cpu", "subgrid": False, "precision": "single", "progressbars": args.log_level <= 20}
+        self.general = {"solver": "cpu", "precision": "single", "progressbars": args.log_level <= 20}

        self.em_consts = {
            "c": c,  # Speed of light in free space (m/s)
            "e0": e0,  # Permittivity of free space (F/m)
            "m0": m0,  # Permeability of free space (H/m)
-            "z0": np.sqrt(m0 / e0),
-        }  # Impedance of free space (Ohms)
+            "z0": np.sqrt(m0 / e0), # Impedance of free space (Ohms)
+        }  

        # Store information about host machine
        self.hostinfo = get_host_info()
@@ -246,7 +242,7 @@ class SimulationConfig:
            # Both single and double precision are possible on GPUs, but single
            # provides best performance.
            self.general["precision"] = "single"
-            self.devices = {"devs": [], "nvcc_opts": None}  # pycuda device objects  # nvcc compiler options
+            self.devices = {"devs": [], "nvcc_opts": None}  # pycuda device objects; nvcc compiler options
            # Suppress nvcc warnings on Microsoft Windows
            if sys.platform == "win32":
                self.devices["nvcc_opts"] = ["-w"]
@@ -258,7 +254,7 @@ class SimulationConfig:
        if self.args.opencl is not None:
            self.general["solver"] = "opencl"
            self.general["precision"] = "single"
-            self.devices = {"devs": [], "compiler_opts": None}  # pyopencl available device(s)
+            self.devices = {"devs": [], "compiler_opts": None}  # pyopencl device device(s); compiler options

            # Suppress unused variable warnings on gcc
            # if sys.platform != 'win32': self.devices['compiler_opts'] = ['-w']
@@ -266,8 +262,8 @@ class SimulationConfig:
            # Add pyopencl available device(s)
            self.devices["devs"] = detect_opencl()

-        # Subgrid parameter may not exist if user enters via CLI
-        try:
+        # Subgrids
+        if hasattr(self.args, "subgrid") and self.args.subgrid:
            self.general["subgrid"] = self.args.subgrid
            # Double precision should be used with subgrid for best accuracy
            self.general["precision"] = "double"
@@ -278,7 +274,7 @@ class SimulationConfig:
                    "You cannot currently use CUDA or OpenCL-based " "solvers with models that contain sub-grids."
                )
                raise ValueError
-        except AttributeError:
+        else:
            self.general["subgrid"] = False

        # Scenes parameter may not exist if user enters via CLI
--- a/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py
+++ b/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py
@@ -35,11 +35,19 @@ x_args = {
                                        int NZ_PHI2,
                                        int NY_R,
                                        const unsigned int* __restrict__ ID,
-                                        const $REAL* __restrict__ Ex, $REAL *Ey,
+                                        const $REAL* __restrict__ Ex, 
+                                        $REAL *Ey,
                                        $REAL *Ez,
-                                        const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1,
+                                        const $REAL* __restrict__ Hx, 
+                                        const $REAL* __restrict__ Hy, 
+                                        const $REAL* __restrict__ Hz, 
+                                        $REAL *PHI1,
                                        $REAL *PHI2,
-                                        const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                        const $REAL* __restrict__ RA, 
+                                        const $REAL* __restrict__ RB, 
+                                        const $REAL* __restrict__ RE, 
+                                        const $REAL* __restrict__ RF, 
+                                        $REAL d)
                    """
    ),
    "opencl": Template(
@@ -95,9 +103,16 @@ y_args = {
                                        $REAL *Ex,
                                        const $REAL* __restrict__ Ey,
                                        $REAL *Ez,
-                                        const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1,
+                                        const $REAL* __restrict__ Hx, 
+                                        const $REAL* __restrict__ Hy, 
+                                        const $REAL* __restrict__ Hz, 
+                                        $REAL *PHI1,
                                        $REAL *PHI2,
-                                        const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                        const $REAL* __restrict__ RA, 
+                                        const $REAL* __restrict__ RB, 
+                                        const $REAL* __restrict__ RE, 
+                                        const $REAL* __restrict__ RF, 
+                                        $REAL d)
                    """
    ),
    "opencl": Template(
@@ -153,9 +168,16 @@ z_args = {
                                        $REAL *Ex,
                                        $REAL *Ey,
                                        const $REAL* __restrict__ Ez,
-                                        const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1,
+                                        const $REAL* __restrict__ Hx,
+                                        const $REAL* __restrict__ Hy, 
+                                        const $REAL* __restrict__ Hz, 
+                                        $REAL *PHI1,
                                        $REAL *PHI2,
-                                        const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                        const $REAL* __restrict__ RA, 
+                                        const $REAL* __restrict__ RB, 
+                                        const $REAL* __restrict__ RE, 
+                                        const $REAL* __restrict__ RF, 
+                                        $REAL d)
                    """
    ),
    "opencl": Template(
--- a/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py
+++ b/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py
@@ -35,11 +35,19 @@ x_args = {
                                        int NZ_PHI2,
                                        int NY_R,
                                        const unsigned int* __restrict__ ID,
-                                        const $REAL* __restrict__ Ex, $REAL *Ey,
+                                        const $REAL* __restrict__ Ex, 
+                                        $REAL *Ey,
                                        $REAL *Ez,
-                                        const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1,
+                                        const $REAL* __restrict__ Hx, 
+                                        const $REAL* __restrict__ Hy, 
+                                        const $REAL* __restrict__ Hz, 
+                                        $REAL *PHI1,
                                        $REAL *PHI2,
-                                        const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                        const $REAL* __restrict__ RA, 
+                                        const $REAL* __restrict__ RB, 
+                                        const $REAL* __restrict__ RE, 
+                                        const $REAL* __restrict__ RF, 
+                                        $REAL d)
                    """
    ),
    "opencl": Template(
@@ -95,9 +103,16 @@ y_args = {
                                        $REAL *Ex,
                                        const $REAL* __restrict__ Ey,
                                        $REAL *Ez,
-                                        const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1,
+                                        const $REAL* __restrict__ Hx, 
+                                        const $REAL* __restrict__ Hy, 
+                                        const $REAL* __restrict__ Hz, 
+                                        $REAL *PHI1,
                                        $REAL *PHI2,
-                                        const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                        const $REAL* __restrict__ RA, 
+                                        const $REAL* __restrict__ RB, 
+                                        const $REAL* __restrict__ RE, 
+                                        const $REAL* __restrict__ RF, 
+                                        $REAL d)
                    """
    ),
    "opencl": Template(
@@ -153,9 +168,16 @@ z_args = {
                                        $REAL *Ex,
                                        $REAL *Ey,
                                        const $REAL* __restrict__ Ez,
-                                        const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1,
+                                        const $REAL* __restrict__ Hx, 
+                                        const $REAL* __restrict__ Hy, 
+                                        const $REAL* __restrict__ Hz, 
+                                        $REAL *PHI1,
                                        $REAL *PHI2,
-                                        const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                        const $REAL* __restrict__ RA, 
+                                        const $REAL* __restrict__ RB, 
+                                        const $REAL* __restrict__ RE, 
+                                        const $REAL* __restrict__ RF, 
+                                        $REAL d)
                    """
    ),
    "opencl": Template(
--- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py
+++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py
@@ -35,11 +35,19 @@ x_args = {
                                          int NZ_PHI2,
                                          int NY_R,
                                          const unsigned int* __restrict__ ID,
-                                          const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, const $REAL* __restrict__ Hx, $REAL *Hy,
+                                          const $REAL* __restrict__ Ex, 
+                                          const $REAL* __restrict__ Ey, 
+                                          const $REAL* __restrict__ Ez, 
+                                          const $REAL* __restrict__ Hx, 
+                                          $REAL *Hy,
                                          $REAL *Hz,
                                          $REAL *PHI1,
                                          $REAL *PHI2,
-                                          const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                          const $REAL* __restrict__ RA, 
+                                          const $REAL* __restrict__ RB, 
+                                          const $REAL* __restrict__ RE, 
+                                          const $REAL* __restrict__ RF, 
+                                          $REAL d)
                    """
    ),
    "opencl": Template(
@@ -92,13 +100,19 @@ y_args = {
                                          int NZ_PHI2,
                                          int NY_R,
                                          const unsigned int* __restrict__ ID,
-                                          const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez,
+                                          const $REAL* __restrict__ Ex, 
+                                          const $REAL* __restrict__ Ey, 
+                                          const $REAL* __restrict__ Ez,
                                          $REAL *Hx,
                                          const $REAL* __restrict__ Hy,
                                          $REAL *Hz,
                                          $REAL *PHI1,
                                          $REAL *PHI2,
-                                          const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                          const $REAL* __restrict__ RA, 
+                                          const $REAL* __restrict__ RB, 
+                                          const $REAL* __restrict__ RE, 
+                                          const $REAL* __restrict__ RF, 
+                                          $REAL d)
                    """
    ),
    "opencl": Template(
@@ -151,13 +165,19 @@ z_args = {
                                          int NZ_PHI2,
                                          int NY_R,
                                          const unsigned int* __restrict__ ID,
-                                          const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez,
+                                          const $REAL* __restrict__ Ex, 
+                                          const $REAL* __restrict__ Ey, 
+                                          const $REAL* __restrict__ Ez,
                                          $REAL *Hx,
                                          $REAL *Hy,
                                          const $REAL* __restrict__ Hz,
                                          $REAL *PHI1,
                                          $REAL *PHI2,
-                                          const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                          const $REAL* __restrict__ RA, 
+                                          const $REAL* __restrict__ RB, 
+                                          const $REAL* __restrict__ RE, 
+                                          const $REAL* __restrict__ RF, 
+                                          $REAL d)
                    """
    ),
    "opencl": Template(
--- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py
+++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py
@@ -35,11 +35,19 @@ x_args = {
                                          int NZ_PHI2,
                                          int NY_R,
                                          const unsigned int* __restrict__ ID,
-                                          const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, const $REAL* __restrict__ Hx, $REAL *Hy,
+                                          const $REAL* __restrict__ Ex, 
+                                          const $REAL* __restrict__ Ey, 
+                                          const $REAL* __restrict__ Ez, 
+                                          const $REAL* __restrict__ Hx, 
+                                          $REAL *Hy,
                                          $REAL *Hz,
                                          $REAL *PHI1,
                                          $REAL *PHI2,
-                                          const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                          const $REAL* __restrict__ RA, 
+                                          const $REAL* __restrict__ RB, 
+                                          const $REAL* __restrict__ RE, 
+                                          const $REAL* __restrict__ RF, 
+                                          $REAL d)
                    """
    ),
    "opencl": Template(
@@ -92,13 +100,19 @@ y_args = {
                                          int NZ_PHI2,
                                          int NY_R,
                                          const unsigned int* __restrict__ ID,
-                                          const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez,
+                                          const $REAL* __restrict__ Ex, 
+                                          const $REAL* __restrict__ Ey, 
+                                          const $REAL* __restrict__ Ez,
                                          $REAL *Hx,
                                          const $REAL* __restrict__ Hy,
                                          $REAL *Hz,
                                          $REAL *PHI1,
                                          $REAL *PHI2,
-                                          const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                          const $REAL* __restrict__ RA, 
+                                          const $REAL* __restrict__ RB, 
+                                          const $REAL* __restrict__ RE, 
+                                          const $REAL* __restrict__ RF, 
+                                          $REAL d)
                    """
    ),
    "opencl": Template(
@@ -151,13 +165,19 @@ z_args = {
                                          int NZ_PHI2,
                                          int NY_R,
                                          const unsigned int* __restrict__ ID,
-                                          const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez,
+                                          const $REAL* __restrict__ Ex, 
+                                          const $REAL* __restrict__ Ey, 
+                                          const $REAL* __restrict__ Ez,
                                          $REAL *Hx,
                                          $REAL *Hy,
                                          const $REAL* __restrict__ Hz,
                                          $REAL *PHI1,
                                          $REAL *PHI2,
-                                          const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d)
+                                          const $REAL* __restrict__ RA, 
+                                          const $REAL* __restrict__ RB, 
+                                          const $REAL* __restrict__ RE, 
+                                          const $REAL* __restrict__ RF, 
+                                          $REAL d)
                    """
    ),
    "opencl": Template(
--- a/gprMax/cuda_opencl/knl_source_updates.py
+++ b/gprMax/cuda_opencl/knl_source_updates.py
@@ -126,10 +126,10 @@ update_magnetic_dipole = {
                                            $REAL dx,
                                            $REAL dy,
                                            $REAL dz,
-                                            __global const int* restrict rcinfo1,
-                                            __global const $REAL* restrict rcinfo2,
-                                            __global const $REAL* restrict rcwaveforms,
-                                            __global const unsigned int* estrict ID,
+                                            __global const int* restrict srcinfo1,
+                                            __global const $REAL* restrict srcinfo2,
+                                            __global const $REAL* restrict srcwaveforms,
+                                            __global const unsigned int* restrict ID,
                                            __global $REAL *Hx,
                                            __global $REAL *Hy,
                                            __global $REAL *Hz
@@ -208,10 +208,10 @@ update_voltage_source = {
                                            $REAL dx,
                                            $REAL dy,
                                            $REAL dz,
-                                            __global const int* restrict rcinfo1,
-                                            __global const $REAL* restrict rcinfo2,
-                                            __global const $REAL* restrict rcwaveforms,
-                                            __global const unsigned int* estrict ID,
+                                            __global const int* restrict srcinfo1,
+                                            __global const $REAL* restrict srcinfo2,
+                                            __global const $REAL* restrict srcwaveforms,
+                                            __global const unsigned int* restrict ID,
                                            __global $REAL *Ex,
                                            __global $REAL *Ey,
                                            __global $REAL *Ez
--- a/gprMax/gprMax.py
+++ b/gprMax/gprMax.py
@@ -25,7 +25,7 @@ from .utilities.logging import logging_config

 # Arguments (used for API) and their default values (used for API and CLI)
 args_defaults = {
-    "scenes": None,
+    "scenes": [],
    "inputfile": None,
    "outputfile": None,
    "n": 1,
@@ -44,37 +44,39 @@ args_defaults = {

 # Argument help messages (used for CLI argparse)
 help_msg = {
-    "scenes": "(list, req): Scenes to run the model. "
-    "Multiple scene objects can given in order to run multiple "
-    "simulation runs. Each scene must contain the essential "
-    "simulation objects",
-    "inputfile": "(str, opt): Input file path. Can also run simulation " "by providing an input file.",
+    "scenes": "(list, req): Scenes to run the model. Multiple scene objects "
+              "can given in order to run multiple simulation runs. Each scene "
+              "must contain the essential simulation objects",
+    "inputfile": "(str, opt): Input file path. Can also run simulation by "
+                 "providing an input file.",
    "outputfile": "(str, req): File path to the output data file.",
    "n": "(int, req): Number of required simulation runs.",
-    "i": "(int, opt): Model number to start/restart simulation "
-    "from. It would typically be used to restart a series of "
-    "models from a specific model number, with the n argument, "
-    "e.g. to restart from A-scan 45 when creating a B-scan "
-    "with 60 traces.",
-    "mpi": "(bool, opt): Flag to use Message Passing Interface (MPI) "
-    "task farm. This option is most usefully combined with n to "
-    "allow individual models to be farmed out using a MPI task "
-    "farm, e.g. to create a B-scan with 60 traces and use MPI to "
-    "farm out each trace. For further details see the parallel "
-    "performance section of the User Guide.",
-    "gpu": "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA " "GPU device ID(s) for specific GPU card(s).",
-    "opencl": "(list/bool, opt): Flag to use OpenCL or list of OpenCL " "device ID(s) for specific compute device(s).",
+    "i": "(int, opt): Model number to start/restart simulation from. It would "
+         "typically be used to restart a series of models from a specific "
+         "model number, with the n argument, e.g. to restart from A-scan 45 "
+         "when creating a B-scan with 60 traces.",
+    "mpi": "(bool, opt): Flag to use Message Passing Interface (MPI) task farm. "
+           "This option is most usefully combined with n to allow individual "
+           "models to be farmed out using a MPI task farm, e.g. to create a "
+           "B-scan with 60 traces and use MPI to farm out each trace. For "
+           "further details see the performance section of the User Guide.",
+    "gpu": "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA GPU "
+           "device ID(s) for specific GPU card(s).",
+    "opencl": "(list/bool, opt): Flag to use OpenCL or list of OpenCL device "
+              "ID(s) for specific compute device(s).",
    "subgrid": "(bool, opt): Flag to use sub-gridding.",
-    "autotranslate": "(bool, opt): For sub-gridding - 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.",
-    "geometry_only": "(bool, opt): Build a model and produce any " "geometry views but do not run the simulation.",
-    "geometry_fixed": "(bool, opt): Run a series of models where the " "geometry does not change between models.",
-    "write_processed": "(bool, opt): Writes another input file after "
-    "any Python code (#python blocks) and in the "
-    "original input file has been processed.",
+    "autotranslate": "(bool, opt): For sub-gridding - 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.",
+    "geometry_only": "(bool, opt): Build a model and produce any geometry "
+                     "views but do not run the simulation.",
+    "geometry_fixed": "(bool, opt): Run a series of models where the geometry "
+                      "does not change between models.",
+    "write_processed": "(bool, opt): Writes another input file after any "
+                       "Python code (#python blocks) and in the original input "
+                       "file has been processed.",
    "log_level": "(int, opt): Level of logging to use.",
    "log_file": "(bool, opt): Write logging information to file.",
 }
@@ -209,9 +211,7 @@ def run_main(args):
    """

    results = {}
-
    logging_config(level=args.log_level, log_file=args.log_file)
-
    config.sim_config = config.SimulationConfig(args)

    # MPI running with (OpenMP/CUDA/OpenCL)
@@ -222,4 +222,5 @@ def run_main(args):
        context = Context()

    results = context.run()
+    
    return results
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -56,10 +56,11 @@ class FDTDGrid:
        self.pmls["formulation"] = "HORIPML"
        self.pmls["cfs"] = []
        self.pmls["slabs"] = []
-        # Ordered dictionary required so that PMLs are always updated in the
-        # same order. The order itself does not matter, however, if must be the
-        # same from model to model otherwise the numerical precision from adding
-        # the PML corrections will be different.
+        # Ordered dictionary required so *updating* the PMLs always follows the
+        # same order (the order for *building* PMLs does not matter). The order 
+        # itself does not matter, however, if must be the same from model to 
+        # model otherwise the numerical precision from adding the PML 
+        # corrections will be different.
        self.pmls["thickness"] = OrderedDict((key, 10) for key in PML.boundaryIDs)

        self.materials = []
@@ -353,6 +354,7 @@ class CUDAGrid(FDTDGrid):
            self.Hx_dev = gpuarray.to_gpu(self.Hx)
            self.Hy_dev = gpuarray.to_gpu(self.Hy)
            self.Hz_dev = gpuarray.to_gpu(self.Hz)
+            
        elif config.sim_config.general["solver"] == "opencl":
            import pyopencl.array as clarray

--- a/gprMax/hash_cmds_file.py
+++ b/gprMax/hash_cmds_file.py
@@ -217,6 +217,7 @@ def check_cmd_names(processedlines, checkessential=True):
            "#title",
            "#omp_threads",
            "#time_step_stability_factor",
+            "#pml_formulation",
            "#pml_cells",
            "#excitation_file",
            "#src_steps",
@@ -248,6 +249,7 @@ def check_cmd_names(processedlines, checkessential=True):
            "#rx",
            "#rx_array",
            "#snapshot",
+            "#pml_cfs",
            "#include_file",
        ]
    }
--- a/gprMax/hash_cmds_multiuse.py
+++ b/gprMax/hash_cmds_multiuse.py
@@ -29,6 +29,7 @@ from .cmds_multiuse import (
    Material,
    MaterialList,
    MaterialRange,
+    PMLCFS,
    Rx,
    RxArray,
    Snapshot,
@@ -375,7 +376,7 @@ def process_multicmds(multicmds):
            tmp = cmdinstance.split()

            if len(tmp) < 2:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least 2 parameters")
+                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least two parameters")
                raise ValueError

            tokens = len(tmp)
@@ -386,4 +387,28 @@ def process_multicmds(multicmds):
            material_list = MaterialList(list_of_materials=lmats, id=tmp[tokens - 1])
            scene_objects.append(material_list)

+    cmdname = '#pml_cfs'
+    if multicmds[cmdname] is not None:        
+        for cmdinstance in multicmds[cmdname]:
+            tmp = cmdinstance.split()
+            
+            if len(tmp) != 12:
+                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly twelve parameters")
+                raise ValueError
+
+            pml_cfs = PMLCFS(alphascalingprofile=tmp[0],
+                             alphascalingdirection=tmp[1],
+                             alphamin=tmp[2],
+                             alphamax=tmp[3],
+                             kappascalingprofile=tmp[4],
+                             kappascalingdirection=tmp[5],
+                             kappamin=tmp[6],
+                             kappamax=tmp[7],
+                             sigmascalingprofile=tmp[8],
+                             sigmascalingdirection=tmp[9],
+                             sigmamin=tmp[10],
+                             sigmamax=tmp[11])
+
+            scene_objects.append(pml_cfs)
+
    return scene_objects
--- a/gprMax/hash_cmds_singleuse.py
+++ b/gprMax/hash_cmds_singleuse.py
@@ -120,20 +120,40 @@ def process_singlecmds(singlecmds):

        scene_objects.append(tw)

+    cmd = "#pml_formulation"
+    if singlecmds[cmd] is not None:
+        tmp = singlecmds[cmd].split()
+        if len(tmp) != 1:
+            logger.exception(f"{cmd} requires one parameter")
+            raise ValueError
+        else:
+            pml_formulation = tmp[0]
+
    cmd = "#pml_cells"
    if singlecmds[cmd] is not None:
        tmp = singlecmds[cmd].split()
        if len(tmp) not in [1, 6]:
            logger.exception(f"{cmd} requires either one or six parameter(s)")
            raise ValueError
-        if len(tmp) == 1:
-            pml_cells = PMLProps(thickness=int(tmp[0]))
+        
+        if 'pml_formulation' in locals():
+            if len(tmp) == 1:
+                pml_props = PMLProps(formulation=pml_formulation, thickness=int(tmp[0]))
+            else:
+                pml_props = PMLProps(
+                    formulation=pml_formulation, 
+                    x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), 
+                    xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5])
+                )
        else:
-            pml_cells = PMLProps(
-                x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5])
-            )
-
-        scene_objects.append(pml_cells)
+            if len(tmp) == 1:
+                pml_props = PMLProps(thickness=int(tmp[0]))
+            else:
+                pml_props = PMLProps(
+                    x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5])
+                )
+            
+        scene_objects.append(pml_props)

    cmd = "#src_steps"
    if singlecmds[cmd] is not None:
--- a/gprMax/pml.py
+++ b/gprMax/pml.py
@@ -433,19 +433,6 @@ class CUDAPML(PML):
            1,
        )

-    def get_update_funcs(self, kernelselectric, kernelsmagnetic):
-        """Gets update functions from PML kernels.
-
-        Args:
-            kernelselectric: pycuda SourceModule containing PML kernels for
-                                electric updates.
-            kernelsmagnetic: pycuda SourceModule containing PML kernels for
-                                magnetic updates.
-        """
-
-        self.update_electric_dev = kernelselectric.get_function("order" + str(len(self.CFS)) + "_" + self.direction)
-        self.update_magnetic_dev = kernelsmagnetic.get_function("order" + str(len(self.CFS)) + "_" + self.direction)
-
    def update_electric(self):
        """Updates electric field components with the PML correction on the GPU."""
        self.update_electric_dev(
@@ -549,12 +536,6 @@ class OpenCLPML(PML):
        self.HPhi1_dev = clarray.to_device(self.queue, self.HPhi1)
        self.HPhi2_dev = clarray.to_device(self.queue, self.HPhi2)

-    def set_blocks_per_grid():
-        pass
-
-    def get_update_funcs():
-        pass
-
    def update_electric(self):
        """Updates electric field components with the PML correction on the
        compute device.
--- a/gprMax/receivers.py
+++ b/gprMax/receivers.py
@@ -26,9 +26,7 @@ class Rx:

    allowableoutputs = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz", "Ix", "Iy", "Iz"]
    defaultoutputs = allowableoutputs[:-3]
-
    allowableoutputs_dev = allowableoutputs[:-3]
-    maxnumoutputs_dev = 0

    def __init__(self):
        self.ID = None
@@ -61,9 +59,6 @@ def htod_rx_arrays(G, queue=None):
        rxcoords[i, 0] = rx.xcoord
        rxcoords[i, 1] = rx.ycoord
        rxcoords[i, 2] = rx.zcoord
-        # Store maximum number of output components
-        if len(rx.outputs) > Rx.maxnumoutputs_dev:
-            Rx.maxnumoutputs_dev = len(rx.outputs)

    # Array to store field components for receivers on compute device -
    #   rows are field components; columns are iterations; pages are receivers
--- a/gprMax/updates.py
+++ b/gprMax/updates.py
@@ -407,7 +407,9 @@ class CUDAUpdates:
            knlH = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
            pml.update_magnetic_dev = knlH.get_function(knl_name)

-        self._copy_mat_coeffs(knlE, knlH)
+            # Copy material coefficient arrays to constant memory of GPU - must
+            # be done for each kernel
+            self._copy_mat_coeffs(knlE, knlH)

    def _set_rx_knl(self):
        """Receivers - initialises arrays on GPU, prepares kernel and gets kernel
--- a/setup.py
+++ b/setup.py
@@ -18,6 +18,7 @@

 import glob
 import os
+import platform
 import re
 import shutil
 import subprocess
@@ -187,35 +188,25 @@ else:
        )
        cpuID = " ".join(cpuID.split())
        if "Apple" in cpuID:
-            gccpath = glob.glob("/opt/homebrew/bin/gcc-[4-9]*")
-            gccpath += glob.glob("/opt/homebrew/bin/gcc-[10-11]*")
-            if gccpath:
-                # Use newest gcc found
-                os.environ["CC"] = gccpath[-1].split(os.sep)[-1]
-                rpath = "/opt/homebrew/opt/gcc/lib/gcc/" + gccpath[-1].split(os.sep)[-1][-1] + "/"
-            else:
-                raise (
-                    "Cannot find gcc in /opt/homebrew/bin. gprMax requires gcc "
-                    + "to be installed - easily done through the Homebrew package "
-                    + "manager (http://brew.sh). Note: gcc with OpenMP support "
-                    + "is required."
-                )
+            gccbasepath = "/opt/homebrew/bin/"
        else:
-            gccpath = glob.glob("/usr/local/bin/gcc-[4-9]*")
-            gccpath += glob.glob("/usr/local/bin/gcc-[10-11]*")
-            if gccpath:
-                # Use newest gcc found
-                os.environ["CC"] = gccpath[-1].split(os.sep)[-1]
-            else:
-                raise (
-                    "Cannot find gcc in /usr/local/bin. gprMax requires gcc "
-                    + "to be installed - easily done through the Homebrew package "
-                    + "manager (http://brew.sh). Note: gcc with OpenMP support "
-                    + "is required."
-                )
+            gccbasepath = "/usr/local/bin/"
+        gccpath = glob.glob(gccbasepath + "gcc-[0-9][0-9]")
+        if gccpath:
+            # Use newest gcc found
+            os.environ["CC"] = gccpath[-1].split(os.sep)[-1]
+            if "Apple" in cpuID:
+                rpath = "/opt/homebrew/opt/gcc/lib/gcc/" + gccpath[-1].split(os.sep)[-1][-1] + "/"
+        else:
+            raise (
+                f"Cannot find gcc in {gccbasepath}. gprMax requires gcc "
+                + "to be installed - easily done through the Homebrew package "
+                + "manager (http://brew.sh). Note: gcc with OpenMP support "
+                + "is required."
+            )

-        # Minimum supported macOS deployment target
-        MIN_MACOS_VERSION = "10.13"
+        # Set minimum supported macOS deployment target to installed macOS version
+        MIN_MACOS_VERSION = platform.mac_ver()[0]
        try:
            os.environ["MACOSX_DEPLOYMENT_TARGET"]
            del os.environ["MACOSX_DEPLOYMENT_TARGET"]
--- a/testing/analytical_solutions.py
+++ b/testing/analytical_solutions.py
@@ -17,8 +17,9 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.

 import numpy as np
+from scipy.constants import c
+from scipy.constants import epsilon_0 as e0

-import gprMax.config as config
 from gprMax.waveforms import Waveform


@@ -78,35 +79,35 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
    Ex_y = y
    Ex_z = z - 0.5 * dz
    Er_x = np.sqrt((Ex_x**2 + Ex_y**2 + Ex_z**2))
-    tau_Ex = Er_x / config.sim_config.em_consts["c"]
+    tau_Ex = Er_x / c

    # Coordinates of Rx for Ey FDTD component
    Ey_x = x
    Ey_y = y + 0.5 * dy
    Ey_z = z - 0.5 * dz
    Er_y = np.sqrt((Ey_x**2 + Ey_y**2 + Ey_z**2))
-    tau_Ey = Er_y / config.sim_config.em_consts["c"]
+    tau_Ey = Er_y / c

    # Coordinates of Rx for Ez FDTD component
    Ez_x = x
    Ez_y = y
    Ez_z = z
    Er_z = np.sqrt((Ez_x**2 + Ez_y**2 + Ez_z**2))
-    tau_Ez = Er_z / config.sim_config.em_consts["c"]
+    tau_Ez = Er_z / c

    # Coordinates of Rx for Hx FDTD component
    Hx_x = x
    Hx_y = y + 0.5 * dy
    Hx_z = z
    Hr_x = np.sqrt((Hx_x**2 + Hx_y**2 + Hx_z**2))
-    tau_Hx = Hr_x / config.sim_config.em_consts["c"]
+    tau_Hx = Hr_x / c

    # Coordinates of Rx for Hy FDTD component
    Hy_x = x + 0.5 * dx
    Hy_y = y
    Hy_z = z
    Hr_y = np.sqrt((Hy_x**2 + Hy_y**2 + Hy_z**2))
-    tau_Hy = Hr_y / config.sim_config.em_consts["c"]
+    tau_Hy = Hr_y / c

    # Initialise fields
    fields = np.zeros((iterations, 6))
@@ -134,7 +135,7 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
        fdot_Hy = wdot.calculate_value((timestep * dt) - tau_Hy, dt) * dl

        # Ex
-        fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_x**5)) * (
+        fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * e0 * Er_x**5)) * (
            3 * (fint_Ex + (tau_Ex * f_Ex)) + (tau_Ex**2 * fdot_Ex)
        )

@@ -145,12 +146,12 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
            tmp = 0
        fields[timestep, 1] = (
            tmp
-            * ((Ey_x * Ey_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_y**5))
+            * ((Ey_x * Ey_z) / (4 * np.pi * e0 * Er_y**5))
            * (3 * (fint_Ey + (tau_Ey * f_Ey)) + (tau_Ey**2 * fdot_Ey))
        )

        # Ez
-        fields[timestep, 2] = (1 / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_z**5)) * (
+        fields[timestep, 2] = (1 / (4 * np.pi * e0 * Er_z**5)) * (
            (2 * Ez_z**2 - (Ez_x**2 + Ez_y**2)) * (fint_Ez + (tau_Ez * f_Ez))
            - (Ez_x**2 + Ez_y**2) * tau_Ez**2 * fdot_Ez
        )
--- a/testing/diff_output_files.py
+++ b/testing/diff_output_files.py
@@ -0,0 +1,93 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
+import logging
+from pathlib import Path
+
+import h5py
+import numpy as np
+
+logger = logging.getLogger(__name__)
+
+def diff_output_files(filename1, filename2):
+    """Calculates differences between two output files.
+
+    Args:
+        filename1: string of filename (including path) of output file 1.
+        filename2: string of filename (including path) of output file 2.
+
+    Returns:
+        time: numpy array containing time.
+        datadiffs: numpy array containing power (dB) of differences.
+    """
+
+    file1 = h5py.File(Path(filename1), "r")
+    file2 = h5py.File(Path(filename2), "r")
+    # Path to receivers in files
+    path = "rxs/rx1/"
+
+    # Get available field output component names
+    outputs1 = list(file1[path].keys())
+    outputs2 = list(file2[path].keys())
+    if outputs1 != outputs2:
+        logger.exception("Field output components are not the same in each file")
+        raise ValueError
+
+    # Check that type of float used to store fields matches
+    floattype1 = file1[path + outputs1[0]].dtype
+    floattype2 = file2[path + outputs2[0]].dtype
+    if floattype1 != floattype2:
+        logger.warning(
+            f"Type of floating point number in test model ({file1[path + outputs1[0]].dtype}) "
+            f"does not match type in reference solution ({file2[path + outputs2[0]].dtype})\n"
+        )
+    
+    # Arrays for storing time
+    time1 = np.zeros((file1.attrs["Iterations"]), dtype=floattype1)
+    time1 = np.linspace(0, (file1.attrs["Iterations"] - 1), num=file1.attrs["Iterations"])
+    time2 = np.zeros((file2.attrs["Iterations"]), dtype=floattype2)
+    time2 = np.linspace(0, (file2.attrs["Iterations"] - 1), num=file2.attrs["Iterations"])
+
+    # Arrays for storing field data
+    data1 = np.zeros((file1.attrs["Iterations"], len(outputs1)), dtype=floattype1)
+    data2 = np.zeros((file2.attrs["Iterations"], len(outputs2)), dtype=floattype2)
+    for ID, name in enumerate(outputs1):
+        data1[:, ID] = file1[path + str(name)][:]
+        data2[:, ID] = file2[path + str(name)][:]
+        if np.any(np.isnan(data1[:, ID])) or np.any(np.isnan(data2[:, ID])):
+            logger.exception("Data contains NaNs")
+            raise ValueError
+
+    file1.close()
+    file2.close()
+
+    # Diffs
+    datadiffs = np.zeros(data1.shape, dtype=np.float64)
+    for i in range(len(outputs2)):
+        maxi = np.amax(np.abs(data1[:, i]))
+        datadiffs[:, i] = np.divide(
+            np.abs(data2[:, i] - data1[:, i]), maxi, out=np.zeros_like(data1[:, i]), where=maxi != 0
+        )  # Replace any division by zero with zero
+
+        # Calculate power (ignore warning from taking a log of any zero values)
+        with np.errstate(divide="ignore"):
+            datadiffs[:, i] = 20 * np.log10(datadiffs[:, i])
+        # Replace any NaNs or Infs from zero division
+        datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0
+
+    return time1, datadiffs
--- a/testing/models_pmls/plot_pml_comparison.py
+++ b/testing/models_pmls/plot_pml_comparison.py
@@ -1,135 +0,0 @@
-# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
-#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
-#
-# This file is part of gprMax.
-#
-# gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# gprMax is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
-
-import itertools
-import logging
-from operator import add
-from pathlib import Path
-
-import h5py
-import matplotlib.pyplot as plt
-import numpy as np
-
-logger = logging.getLogger(__name__)
-
-# Create/setup plot figure
-# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php
-# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"]
-colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"]
-# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"]
-colors = itertools.cycle(colorIDs)
-# for i in range(2):
-#     next(colors)
-lines = itertools.cycle(("--", ":", "-.", "-"))
-markers = ["o", "d", "^", "s", "*"]
-
-parts = Path(__file__).parts
-path = "rxs/rx1/"
-basename = "pml_3D_pec_plate"
-PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"]
-maxerrors = []
-testmodels = ["pml_3D_pec_plate_" + s for s in PMLIDs]
-
-fig, ax = plt.subplots(
-    subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w"
-)
-
-for x, model in enumerate(testmodels):
-    # Open output file and read iterations
-    fileref = h5py.File(Path(*parts[:-1], basename, basename + "_ref.h5"), "r")
-    filetest = h5py.File(Path(*parts[:-1], basename, basename + str(x + 1) + ".h5"), "r")
-
-    # Get available field output component names
-    outputsref = list(fileref[path].keys())
-    outputstest = list(filetest[path].keys())
-    if outputsref != outputstest:
-        logger.exception("Field output components do not match reference solution")
-        raise ValueError
-
-    # Check that type of float used to store fields matches
-    if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype:
-        logger.warning(
-            f"Type of floating point number in test model ({filetest[path + outputstest[0]].dtype}) "
-            f"does not match type in reference solution ({fileref[path + outputsref[0]].dtype})\n"
-        )
-    floattyperef = fileref[path + outputsref[0]].dtype
-    floattypetest = filetest[path + outputstest[0]].dtype
-    # logger.info(f'Data type: {floattypetest}')
-
-    # Arrays for storing time
-    # timeref = np.zeros((fileref.attrs['Iterations']), dtype=floattyperef)
-    # timeref = np.linspace(0, (fileref.attrs['Iterations'] - 1) * fileref.attrs['dt'], num=fileref.attrs['Iterations']) / 1e-9
-    # timetest = np.zeros((filetest.attrs['Iterations']), dtype=floattypetest)
-    # timetest = np.linspace(0, (filetest.attrs['Iterations'] - 1) * filetest.attrs['dt'], num=filetest.attrs['Iterations']) / 1e-9
-    timeref = np.zeros((fileref.attrs["Iterations"]), dtype=floattyperef)
-    timeref = np.linspace(0, (fileref.attrs["Iterations"] - 1), num=fileref.attrs["Iterations"])
-    timetest = np.zeros((filetest.attrs["Iterations"]), dtype=floattypetest)
-    timetest = np.linspace(0, (filetest.attrs["Iterations"] - 1), num=filetest.attrs["Iterations"])
-
-    # Arrays for storing field data
-    dataref = np.zeros((fileref.attrs["Iterations"], len(outputsref)), dtype=floattyperef)
-    datatest = np.zeros((filetest.attrs["Iterations"], len(outputstest)), dtype=floattypetest)
-    for ID, name in enumerate(outputsref):
-        dataref[:, ID] = fileref[path + str(name)][:]
-        datatest[:, ID] = filetest[path + str(name)][:]
-        if np.any(np.isnan(datatest[:, ID])):
-            logger.exception("Test data contains NaNs")
-            raise ValueError
-
-    fileref.close()
-    filetest.close()
-
-    # Diffs
-    datadiffs = np.zeros(datatest.shape, dtype=np.float64)
-    for i in range(len(outputstest)):
-        maxi = np.amax(np.abs(dataref[:, i]))
-        datadiffs[:, i] = np.divide(
-            np.abs(datatest[:, i] - dataref[:, i]), maxi, out=np.zeros_like(dataref[:, i]), where=maxi != 0
-        )  # Replace any division by zero with zero
-
-        # Calculate power (ignore warning from taking a log of any zero values)
-        with np.errstate(divide="ignore"):
-            datadiffs[:, i] = 20 * np.log10(datadiffs[:, i])
-        # Replace any NaNs or Infs from zero division
-        datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0
-
-    # Print maximum error value
-    start = 210
-    maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]")
-    logger.info(f"{model}: Max. error {maxerrors[x]}")
-
-    # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..)
-    ax.plot(timeref[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=model)
-    ax.set_xticks(np.arange(0, 2200, step=100))
-    ax.set_xlim([0, 2100])
-    ax.set_yticks(np.arange(-160, 0, step=20))
-    ax.set_ylim([-160, -20])
-    ax.set_axisbelow(True)
-    ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed")
-
-mylegend = list(map(add, PMLIDs, maxerrors))
-legend = ax.legend(mylegend, loc=1, fontsize=14)
-frame = legend.get_frame()
-frame.set_edgecolor("white")
-frame.set_alpha(0)
-
-plt.show()
-
-# Save a PDF/PNG of the figure
-# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
-# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
--- a/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py
+++ b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py
@@ -0,0 +1,80 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
+import itertools
+import logging
+from operator import add
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from testing.diff_output_files import diff_output_files
+
+logger = logging.getLogger(__name__)
+
+# Create/setup plot figure
+# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php
+# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"]
+colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"]
+# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"]
+colors = itertools.cycle(colorIDs)
+# for i in range(2):
+#     next(colors)
+lines = itertools.cycle(("--", ":", "-.", "-"))
+markers = ["o", "d", "^", "s", "*"]
+
+fn = Path(__file__)
+basename = "pml_3D_pec_plate"
+PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"]
+maxerrors = []
+testmodels = [basename + "_" + s for s in PMLIDs]
+
+fig, ax = plt.subplots(
+    subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w"
+)
+
+for x, model in enumerate(testmodels):
+    time, datadiffs = diff_output_files(fn.parent.joinpath(basename + "_ref.h5"), 
+                                        fn.parent.joinpath(basename + str(x + 1) + ".h5"))
+
+    # Print maximum error value
+    start = 210
+    maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]")
+    logger.info(f"{model}: Max. error {maxerrors[x]}")
+
+    # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..)
+    ax.plot(time[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=model)
+    ax.set_xticks(np.arange(0, 2200, step=100))
+    ax.set_xlim([0, 2100])
+    ax.set_yticks(np.arange(-160, 0, step=20))
+    ax.set_ylim([-160, -20])
+    ax.set_axisbelow(True)
+    ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed")
+
+mylegend = list(map(add, PMLIDs, maxerrors))
+legend = ax.legend(mylegend, loc=1, fontsize=14)
+frame = legend.get_frame()
+frame.set_edgecolor("white")
+frame.set_alpha(0)
+
+plt.show()
+
+# Save a PDF/PNG of the figure
+# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
+# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
--- a/testing/models_pmls/pml_basic/plot_pml_comparison.py
+++ b/testing/models_pmls/pml_basic/plot_pml_comparison.py
@@ -0,0 +1,82 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
+import itertools
+import logging
+from operator import add
+from pathlib import Path
+
+import matplotlib.gridspec as gridspec
+import matplotlib.pyplot as plt
+import numpy as np
+
+from testing.diff_output_files import diff_output_files
+
+logger = logging.getLogger(__name__)
+
+# Create/setup plot figure
+# Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php
+colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545"]
+colors = itertools.cycle(colorIDs)
+lines = itertools.cycle(("--", ":", "-.", "-"))
+markers = ["o", "d", "^", "s", "*"]
+
+fn = Path(__file__)
+basename = "pml_basic"
+PMLIDs = ["off", "x0", "y0", "z0", "xmax", "ymax", "zmax"]
+maxerrors = []
+
+for x, PMLID in enumerate(PMLIDs):
+    file1 = fn.parent.joinpath(basename + str(x + 1) + "_CPU.h5")
+    file2 = fn.parent.joinpath(basename + str(x + 1) + "_GPU.h5")
+    time, datadiffs = diff_output_files(file1, file2)
+
+    # Print maximum error value
+    start = 150
+    maxerrors.append(f": {np.amax(np.amax(datadiffs[start::, :])):.1f} [dB]")
+    print(f"{PMLID}: Max. error {maxerrors[x]}")
+
+    # Plot diffs
+    gs = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3)
+    fig, ax = plt.subplots(figsize=(20, 10), facecolor="w", edgecolor="w")
+    ax.remove()
+    fig.suptitle(f"{PMLID}")
+    
+    outputs = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"]
+    for i, output in enumerate(outputs):
+        if i < 3:
+            ax = plt.subplot(gs[i, 0])
+        else:
+            ax = plt.subplot(gs[i - 3, 1])
+        ax.plot(time[start::], datadiffs[start::, i], color=next(colors), lw=2)
+        ax.set_xticks(np.arange(0, 1800, step=200))
+        ax.set_xlim([0, 1600])
+        ax.set_yticks(np.arange(-400, 80, step=40))
+        ax.set_ylim([-400, 40])
+        ax.set_axisbelow(True)
+        ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed")
+        ax.set_xlabel("Time [iterations]")
+        ax.set_ylabel(f"{output} error [dB]")
+
+    # Save a PDF/PNG of the figure
+    fig.savefig(basename + "_diffs_" + PMLID + ".pdf", dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
+    # fig.savefig(basename + "_diffs_" + PMLID + ".png", dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
+
+plt.show()
+
+
--- a/testing/models_pmls/pml_basic/pml_basic.py
+++ b/testing/models_pmls/pml_basic/pml_basic.py
--- a/testing/test_models.py
+++ b/testing/test_models.py
@@ -40,14 +40,12 @@ if sys.platform == "linux":
        python -m testing.test_models
 """

-# Specify directoty with set of models to test
+# Specify directory with set of models to test
 modelset = "models_basic"
 # modelset += 'models_advanced'
-# modelset += 'models_pmls'

 basepath = Path(__file__).parents[0] / modelset

-
 # List of available basic test models
 testmodels = [
    "hertzian_dipole_fs_analytical",
@@ -64,9 +62,6 @@ testmodels = [
 # List of available advanced test models
 # testmodels = ['antenna_GSSI_1500_fs', 'antenna_MALA_1200_fs']

-# List of available PML models
-# testmodels = ['pml_x0', 'pml_y0', 'pml_z0', 'pml_xmax', 'pml_ymax', 'pml_zmax', 'pml_3D_pec_plate']
-
 # Select a specific model if desired
 # testmodels = [testmodels[0]]
 testresults = dict.fromkeys(testmodels)
@@ -80,7 +75,7 @@ for i, model in enumerate(testmodels):

    # Run model
    file = basepath / model / model
-    gprMax.run(inputfile=file.with_suffix(".in"), gpu=None)
+    gprMax.run(inputfile=file.with_suffix(".in"), gpu=None, opencl=None)

    # Special case for analytical comparison
    if model == "hertzian_dipole_fs_analytical":
@@ -116,6 +111,7 @@ for i, model in enumerate(testmodels):
        dataref = hertzian_dipole_fs(
            filetest.attrs["Iterations"], filetest.attrs["dt"], filetest.attrs["dx_dy_dz"], rxposrelative
        )
+        filetest.close()

    else:
        # Get output for model and reference files
@@ -134,15 +130,15 @@ for i, model in enumerate(testmodels):
            raise ValueError

        # Check that type of float used to store fields matches
-        if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype:
-            logger.warning(
-                f"Type of floating point number in test model "
-                f"({filetest[path + outputstest[0]].dtype}) does not "
-                f"match type in reference solution ({fileref[path + outputsref[0]].dtype})\n"
-            )
        float_or_doubleref = fileref[path + outputsref[0]].dtype
        float_or_doubletest = filetest[path + outputstest[0]].dtype
-
+        if float_or_doubleref != float_or_doubletest:
+            logger.warning(
+                f"Type of floating point number in test model "
+                f"({float_or_doubletest}) does not "
+                f"match type in reference solution ({float_or_doubleref})\n"
+            )
+        
        # Arrays for storing time
        timeref = np.zeros((fileref.attrs["Iterations"]), dtype=float_or_doubleref)
        timeref = (
@@ -166,7 +162,7 @@ for i, model in enumerate(testmodels):
                raise ValueError

        fileref.close()
-    filetest.close()
+        filetest.close()

    # Diffs
    datadiffs = np.zeros(datatest.shape, dtype=np.float64)
--- a/toolboxes/GPRAntennaModels/GSSI.py
+++ b/toolboxes/GPRAntennaModels/GSSI.py
@@ -420,12 +420,12 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
    return scene_objects


-def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs):
+def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
    """Inserts a description of an antenna similar to the GSSI 400MHz antenna.
        This model represents an update to the previous model of the GSSI 400MHz 
        antenna and was created and optimised by Stadler et al. (2022) 
        in: https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9686638. 
-        Can be used with 0.5mm, 1mm (default) or 2mm spatial resolution.
+        Can be used with 2mm spatial resolution.
        The external dimensions of the antenna are 300x300x178mm.
        One output point is defined between the arms of the receiver bowtie.
        The bowties are aligned with the y axis so the output is the y component
@@ -451,7 +451,6 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs):
    casesize = (0.3, 0.3, 0.178)  # original
    casethickness = 0.002
    shieldthickness = 0.002
-    foamsurroundthickness = 0.003
    pcbthickness = 0.002
    bowtiebase = 0.06
    bowtieheight = 0.06  # original 0.056
@@ -488,25 +487,12 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs):
    # Coordinates of source excitation point in antenna
    tx = x + 0.01 + 0.005 + 0.056, y + casethickness + 0.005 + 0.143, z + skidthickness

-    if resolution == 0.0005:
-        dx = 0.0005
-        dy = 0.0005
-        dz = 0.0005
-        tx = x + 0.01 + 0.005 + 0.056, y + casethickness + 0.005 + 0.1435, z + skidthickness
-    elif resolution == 0.001:
-        dx = 0.001
-        dy = 0.001
-        dz = 0.001
-    elif resolution == 0.002:
-        dx = 0.002
-        dy = 0.002
-        dz = 0.002
-        foamsurroundthickness = 0.002
-        metalboxheight = 0.088
-        tx = x + 0.01 + 0.004 + 0.056, y + casethickness + 0.005 + 0.143 - 0.002, z + skidthickness - 0.002
-    else:
-        logger.exception("This antenna module can only be used with a spatial discretisation of 0.5mm, 1mm, 2mm")
-        raise ValueError
+    dx = 0.002
+    dy = 0.002
+    dz = 0.002
+    foamsurroundthickness = 0.002
+    metalboxheight = 0.088
+    tx = x + 0.01 + 0.004 + 0.056, y + casethickness + 0.005 + 0.143 - 0.002, z + skidthickness - 0.002

    # Material definitions
    absorber = gprMax.Material(er=absorberEr, se=absorbersig, mr=1, sm=0, id="absorber")
@@ -565,68 +551,25 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs):
        scene_objects.extend((b1, b2, b3, b4, b5))

        # PCB
-        if resolution == 0.0005:
-            b6 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.0235, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.034 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-            )
-            b7 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.0235, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.194 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-            )
-            scene_objects.extend((b6, b7))
-
-        elif resolution == 0.001:
-            b6 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.023, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.034 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-            )
-            b7 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.023, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.194 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-            )
-            scene_objects.extend((b6, b7))
-
-        elif resolution == 0.002:
-            b6 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
-                p2=(
-                    x + 0.01 + 0.005 + 0.033 + bowtiebase,
-                    y + casethickness + 0.006 + 0.202 + patchheight,
-                    z + skidthickness + pcbthickness - 0.002,
-                ),
-                material_id="pcb",
-            )
-            b7 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
-                p2=(
-                    x + 0.01 + 0.005 + 0.195 + bowtiebase,
-                    y + casethickness + 0.006 + 0.202 + patchheight,
-                    z + skidthickness + pcbthickness - 0.002,
-                ),
-                material_id="pcb",
-            )
-            scene_objects.extend((b6, b7))
+        b6 = gprMax.Box(
+            p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+            p2=(
+                x + 0.01 + 0.005 + 0.033 + bowtiebase,
+                y + casethickness + 0.006 + 0.202 + patchheight,
+                z + skidthickness + pcbthickness - 0.002,
+            ),
+            material_id="pcb",
+        )
+        b7 = gprMax.Box(
+            p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+            p2=(
+                x + 0.01 + 0.005 + 0.195 + bowtiebase,
+                y + casethickness + 0.006 + 0.202 + patchheight,
+                z + skidthickness + pcbthickness - 0.002,
+            ),
+            material_id="pcb",
+        )
+        scene_objects.extend((b6, b7))

    elif smooth_dec == "no":
        # Plastic case
@@ -682,366 +625,124 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs):
        scene_objects.extend((b8, b9, b10, b11, b12))

        # PCB
-        if resolution == 0.0005:
-            b13 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.0235, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.034 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-                averaging="n",
-            )
-            b14 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.0235, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.194 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-                averaging="n",
-            )
-            scene_objects.extend((b13, b14))
-
-        elif resolution == 0.001:
-            b13 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.023, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.034 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-                averaging="n",
-            )
-            b14 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.023, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.194 + bowtiebase,
-                    y + casethickness + 0.005 + 0.204 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-                averaging="n",
-            )
-            scene_objects.extend((b13, b14))
-
-        elif resolution == 0.002:
-            b13 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
-                p2=(
-                    x + 0.01 + 0.005 + 0.033 + bowtiebase,
-                    y + casethickness + 0.006 + 0.202 + patchheight,
-                    z + skidthickness + pcbthickness - 0.002,
-                ),
-                material_id="pcb",
-                averaging="n",
-            )
-            b14 = gprMax.Box(
-                p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness),
-                p2=(
-                    x + 0.01 + 0.005 + 0.195 + bowtiebase,
-                    y + casethickness + 0.006 + 0.202 + patchheight,
-                    z + skidthickness + pcbthickness,
-                ),
-                material_id="pcb",
-                averaging="n",
-            )
-            scene_objects.extend((b13, b14))
+        b13 = gprMax.Box(
+            p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+            p2=(
+                x + 0.01 + 0.005 + 0.033 + bowtiebase,
+                y + casethickness + 0.006 + 0.202 + patchheight,
+                z + skidthickness + pcbthickness - 0.002,
+            ),
+            material_id="pcb",
+            averaging="n",
+        )
+        b14 = gprMax.Box(
+            p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness),
+            p2=(
+                x + 0.01 + 0.005 + 0.195 + bowtiebase,
+                y + casethickness + 0.006 + 0.202 + patchheight,
+                z + skidthickness + pcbthickness,
+            ),
+            material_id="pcb",
+            averaging="n",
+        )
+        scene_objects.extend((b13, b14))

    # PCB components
    # My own bowties with triangle commands
-    if resolution == 0.0005:
-        # "left" side
-        # extension plates
-        p1 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.0235, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.026 + bowtiebase,
-                y + casethickness + 0.005 + 0.0235 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        p2 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.026 + bowtiebase,
-                y + casethickness + 0.005 + 0.204 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        # triangles
-        t1 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.0835, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.0835, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.0835 + bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-        t2 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.204 - bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
+    # "left" side
+    # extension plates
+    p1 = gprMax.Plate(
+        p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+        p2=(
+            x + 0.01 + 0.005 + 0.025 + bowtiebase,
+            y + casethickness + 0.005 + 0.021 + patchheight,
+            z + skidthickness - 0.002,
+        ),
+        material_id="pec",
+    )
+    p2 = gprMax.Plate(
+        p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
+        p2=(
+            x + 0.01 + 0.005 + 0.025 + bowtiebase,
+            y + casethickness + 0.005 + 0.203 + patchheight,
+            z + skidthickness - 0.002,
+        ),
+        material_id="pec",
+    )
+    # triangles
+    t1 = gprMax.Triangle(
+        p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
+        p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
+        p3=(
+            x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2),
+            y + casethickness + 0.005 + 0.081 + bowtieheight,
+            z + skidthickness - 0.002,
+        ),
+        thickness=0,
+        material_id="pec",
+    )
+    t2 = gprMax.Triangle(
+        p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
+        p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
+        p3=(
+            x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2),
+            y + casethickness + 0.005 + 0.203 - bowtieheight,
+            z + skidthickness - 0.002,
+        ),
+        thickness=0,
+        material_id="pec",
+    )
+    # "right" side
+    p3 = gprMax.Plate(
+        p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+        p2=(
+            x + 0.01 + 0.005 + 0.187 + bowtiebase,
+            y + casethickness + 0.005 + 0.021 + patchheight,
+            z + skidthickness - 0.002,
+        ),
+        material_id="pec",
+    )
+    p4 = gprMax.Plate(
+        p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
+        p2=(
+            x + 0.01 + 0.005 + 0.187 + bowtiebase,
+            y + casethickness + 0.005 + 0.203 + patchheight,
+            z + skidthickness - 0.002,
+        ),
+        material_id="pec",
+    )
+    # triangles
+    t3 = gprMax.Triangle(
+        p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
+        p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
+        p3=(
+            x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2),
+            y + casethickness + 0.005 + 0.081 + bowtieheight,
+            z + skidthickness - 0.002,
+        ),
+        thickness=0,
+        material_id="pec",
+    )
+    t4 = gprMax.Triangle(
+        p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
+        p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
+        p3=(
+            x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2),
+            y + casethickness + 0.005 + 0.203 - bowtieheight,
+            z + skidthickness - 0.002,
+        ),
+        thickness=0,
+        material_id="pec",
+    )

-        # "right" side
-        p3 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.0235, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.186 + bowtiebase,
-                y + casethickness + 0.005 + 0.0235 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        p4 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.186 + bowtiebase,
-                y + casethickness + 0.005 + 0.204 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        # triangles
-        t3 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.0835, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.0835, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.0835 + bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_ID="pec",
-        )
-        t4 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.204 - bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-
-        # Edges that represent wire between bowtie halves in 1mm model
-        e1 = gprMax.Edge(p1=(tx[0] + 0.16, tx[1] - dy, tx[2]), p2=(tx[0] + 0.16, tx[1], tx[2]), material_id="pec")
-        e2 = gprMax.Edge(
-            p1=(tx[0] + 0.16, tx[1] + dy, tx[2]), p2=(tx[0] + 0.16, tx[1] + 2 * dy, tx[2]), material_id="pec"
-        )
-        e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec")
-        e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec")
-        scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4))
-
-    elif resolution == 0.001:
-        # "left" side
-        # extension plates
-        p1 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.023, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.026 + bowtiebase,
-                y + casethickness + 0.005 + 0.023 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        p2 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.026 + bowtiebase,
-                y + casethickness + 0.005 + 0.204 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        # triangles
-        t1 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.083, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.083, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.083 + bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-        t2 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.204 - bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-
-        # "right" side
-        p3 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.023, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.186 + bowtiebase,
-                y + casethickness + 0.005 + 0.023 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        p4 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(
-                x + 0.01 + 0.005 + 0.186 + bowtiebase,
-                y + casethickness + 0.005 + 0.204 + patchheight,
-                z + skidthickness,
-            ),
-            material_id="pec",
-        )
-        # triangles
-        t3 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.083, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.083, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.083 + bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-        t4 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness),
-            p3=(
-                x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.204 - bowtieheight,
-                z + skidthickness,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-
-        # Edges that represent wire between bowtie halves in 1mm model
-        e1 = gprMax.Edge(p1=(tx[0] + 0.16, tx[1] - dy, tx[2]), p2=(tx[0] + 0.16, tx[1], tx[2]), material_id="pec")
-        e2 = gprMax.Edge(
-            p1=(tx[0] + 0.16, tx[1] + dy, tx[2]), p2=(tx[0] + 0.16, tx[1] + 2 * dy, tx[2]), material_id="pec"
-        )
-        e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec")
-        e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec")
-        scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4))
-
-    elif resolution == 0.002:
-        # "left" side
-        # extension plates
-        p1 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
-            p2=(
-                x + 0.01 + 0.005 + 0.025 + bowtiebase,
-                y + casethickness + 0.005 + 0.021 + patchheight,
-                z + skidthickness - 0.002,
-            ),
-            material_id="pec",
-        )
-        p2 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
-            p2=(
-                x + 0.01 + 0.005 + 0.025 + bowtiebase,
-                y + casethickness + 0.005 + 0.203 + patchheight,
-                z + skidthickness - 0.002,
-            ),
-            material_id="pec",
-        )
-        # triangles
-        t1 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
-            p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
-            p3=(
-                x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.081 + bowtieheight,
-                z + skidthickness - 0.002,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-        t2 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
-            p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
-            p3=(
-                x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.203 - bowtieheight,
-                z + skidthickness - 0.002,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-        # "right" side
-        p3 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
-            p2=(
-                x + 0.01 + 0.005 + 0.187 + bowtiebase,
-                y + casethickness + 0.005 + 0.021 + patchheight,
-                z + skidthickness - 0.002,
-            ),
-            material_id="pec",
-        )
-        p4 = gprMax.Plate(
-            p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
-            p2=(
-                x + 0.01 + 0.005 + 0.187 + bowtiebase,
-                y + casethickness + 0.005 + 0.203 + patchheight,
-                z + skidthickness - 0.002,
-            ),
-            material_id="pec",
-        )
-        # triangles
-        t3 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
-            p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
-            p3=(
-                x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.081 + bowtieheight,
-                z + skidthickness - 0.002,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-        t4 = gprMax.Triangle(
-            p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
-            p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
-            p3=(
-                x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2),
-                y + casethickness + 0.005 + 0.203 - bowtieheight,
-                z + skidthickness - 0.002,
-            ),
-            thickness=0,
-            material_id="pec",
-        )
-
-        # Edges that represent wire between bowtie halves in 2mm model
-        e1 = gprMax.Edge(p1=(tx[0] + 0.162, tx[1] - dy, tx[2]), p2=(tx[0] + 0.162, tx[1], tx[2]), material_id="pec")
-        e2 = gprMax.Edge(
-            p1=(tx[0] + 0.162, tx[1] + dy, tx[2]), p2=(tx[0] + 0.162, tx[1] + 2 * dy, tx[2]), material_id="pec"
-        )
-        e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec")
-        e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec")
-        scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4))
+    # Edges that represent wire between bowtie halves in 2mm model
+    e1 = gprMax.Edge(p1=(tx[0] + 0.162, tx[1] - dy, tx[2]), p2=(tx[0] + 0.162, tx[1], tx[2]), material_id="pec")
+    e2 = gprMax.Edge(
+        p1=(tx[0] + 0.162, tx[1] + dy, tx[2]), p2=(tx[0] + 0.162, tx[1] + 2 * dy, tx[2]), material_id="pec"
+    )
+    e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec")
+    e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec")
+    scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4))

    # Metallic plate extension
    b15 = gprMax.Box(
@@ -1069,20 +770,18 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs):
    scene_objects.extend((b15, b16))

    # Source
-    # Excitation - Gaussian pulse
-    w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian")
-    scene_objects.append(w1)
-
    if src_type == "voltage_source":
+        w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian")
        vs1 = gprMax.VoltageSource(
            polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_gaussian"
        )
-        scene_objects.append(vs1)
+        scene_objects.extend((w1, vs1))
    elif src_type == "transmission_line":
+        w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian")
        tl1 = gprMax.TransmissionLine(
            polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_gaussian"
        )
-        scene_objects.append(tl1)
+        scene_objects.extend((w1, tl1))
    else:
        # Optimised custom pulse
        exc1 = gprMax.ExcitationFile(
@@ -1094,35 +793,19 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs):
        scene_objects.extend((exc1, vs1))

    # Receiver
-    # Zero waveform to use with transmission line at receiver output
-    w2 = gprMax.Waveform(wave_type="gaussian", amp=0, freq=excitationfreq, id="my_zero_wave")
-    scene_objects.append(w2)
-
-    if resolution == 0.001 or resolution == 0.0005:
-        if src_type == "transmission_line":
-            tl2 = gprMax.TransmissionLine(
-                polarisation="y",
-                p1=(tx[0] + 0.16, tx[1], tx[2]),
-                resistance=receiverresistance,
-                waveform_id="my_zero_wave",
-            )
-            scene_objects.append(tl2)
-        elif src_type == "voltage_source":
-            r1 = gprMax.Rx(p1=(tx[0] + 0.16, tx[1], tx[2]), id="rxbowtie", outputs="Ey")
-            scene_objects.append(r1)
-
-    elif resolution == 0.002:
-        if src_type == "transmission_line":
-            tl2 = gprMax.TransmissionLine(
-                polarisation="y",
-                p1=(tx[0] + 0.162, tx[1], tx[2]),
-                resistance=receiverresistance,
-                waveform_id="my_zero_wave",
-            )
-            scene_objects.append(tl2)
-        elif src_type == "voltage_source":
-            r1 = gprMax.Rx(p1=(tx[0] + 0.162, tx[1], tx[2]), id="rxbowtie", outputs="Ey")
-            scene_objects.append(r1)
+    if src_type == "transmission_line":
+        # Zero waveform to use with transmission line at receiver output
+        w2 = gprMax.Waveform(wave_type="gaussian", amp=0, freq=excitationfreq, id="my_zero_wave")
+        tl2 = gprMax.TransmissionLine(
+            polarisation="y",
+            p1=(tx[0] + 0.162, tx[1], tx[2]),
+            resistance=receiverresistance,
+            waveform_id="my_zero_wave",
+        )
+        scene_objects.extend((w2, tl2))
+    else:
+        r1 = gprMax.Rx(p1=(tx[0] + 0.162, tx[1], tx[2]), id="rxbowtie", outputs="Ey")
+        scene_objects.append(r1)

    # Geometry views
    gv1 = gprMax.GeometryView(
--- a/toolboxes/GPRAntennaModels/README.rst
+++ b/toolboxes/GPRAntennaModels/README.rst
@@ -14,7 +14,7 @@ Manufacturer/Model       Dimensions    Resolution(s) Author/Contact
 ======================== ============= ============= ========================================================================================================================================================================================================================= ================
 GSSI 1.5GHz (Model 5100) 170x108x45mm  1, 2mm        Craig Warren (craig.warren@northumbria.ac.uk), Northumbria University, UK                                                                                                                                                 1,2
 MALA 1.2GHz              184x109x46mm  1, 2mm        Craig Warren (craig.warren@northumbria.ac.uk), Northumbria University, UK                                                                                                                                                 1
-GSSI 400MHz              300x300x170mm 0.5, 1, 2mm   Sam Stadler (Sam.Stadler@liag-hannover.de), `Leibniz Institute for Applied Geophysics <https://www.leibniz-liag.de/en/research/methods/electromagnetic-methods/ground-penetrating-radar/guided-gpr-waves.html>`_, Germany 3
+GSSI 400MHz              300x300x170mm 2mm           Sam Stadler (sam.stadler@leibniz-liag.de), `Leibniz Institute for Applied Geophysics <https://www.leibniz-liag.de/en/research/methods/electromagnetic-methods/ground-penetrating-radar.html>`_, Germany 3
 ======================== ============= ============= ========================================================================================================================================================================================================================= ================

 **License**: `Creative Commons Attribution-ShareAlike 4.0 International License <http://creativecommons.org/licenses/by-sa/4.0/>`_
@@ -36,7 +36,7 @@ Descriptions of how the models were created can be found in the aforementioned a
 How to use the package
 ======================

-The antenna models can be accessed from within a block of Python code in an input file. The models are inserted at location x,y,z. The coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction. The models must be used with cubic spatial resolutions of either 0.5mm (GSSI 400MHz antenna only), 1mm (default), or 2mm by setting the keyword argument, e.g. ``resolution=0.002``.
+The antenna models can be accessed from within a block of Python code in an input file. The models are inserted at location x,y,z. The coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction. The models must be used with cubic spatial resolutions of either 1mm or 2mm by setting the keyword argument, e.g. ``resolution=0.002``.

 .. note::