Merge branch 'add_reframe_tests' into mpi

gprMax/hash_cmds_multiuse.py

@@ -61,7 +61,9 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) != 4:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly four parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly four parameters"
+                )
                 raise ValueError
 
             waveform = Waveform(wave_type=tmp[0], amp=float(tmp[1]), freq=float(tmp[2]), id=tmp[3])
@@ -88,7 +90,9 @@ def process_multicmds(multicmds):
                     end=float(tmp[7]),
                 )
             else:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least six parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least six parameters"
+                )
                 raise ValueError
 
             scene_objects.append(voltage_source)
@@ -98,11 +102,20 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) < 5:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least five parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires at least five parameters"
+                )
                 raise ValueError
             if len(tmp) == 5:
                 hertzian_dipole = HertzianDipole(
-                    polarisation=tmp[0], p1=(float(tmp[1]), float(tmp[2]), float(tmp[3])), waveform_id=tmp[4]
+                    polarisation=tmp[0],
+                    p1=(float(tmp[1]), float(tmp[2]), float(tmp[3])),
+                    waveform_id=tmp[4],
                 )
             elif len(tmp) == 7:
                 hertzian_dipole = HertzianDipole(
@@ -113,7 +126,9 @@ def process_multicmds(multicmds):
                     end=float(tmp[6]),
                 )
             else:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters"
+                )
                 raise ValueError
 
             scene_objects.append(hertzian_dipole)
@@ -123,11 +138,20 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) < 5:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least five parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires at least five parameters"
+                )
                 raise ValueError
             if len(tmp) == 5:
                 magnetic_dipole = MagneticDipole(
-                    polarisation=tmp[0], p1=(float(tmp[1]), float(tmp[2]), float(tmp[3])), waveform_id=tmp[4]
+                    polarisation=tmp[0],
+                    p1=(float(tmp[1]), float(tmp[2]), float(tmp[3])),
+                    waveform_id=tmp[4],
                 )
             elif len(tmp) == 7:
                 magnetic_dipole = MagneticDipole(
@@ -138,7 +162,9 @@ def process_multicmds(multicmds):
                     end=float(tmp[6]),
                 )
             else:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters"
+                )
                 raise ValueError
 
             scene_objects.append(magnetic_dipole)
@@ -148,7 +174,9 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) < 6:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least six parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least six parameters"
+                )
                 raise ValueError
 
             if len(tmp) == 6:
@@ -168,17 +196,19 @@ def process_multicmds(multicmds):
                     end=tmp[7],
                 )
             else:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters"
+                )
                 raise ValueError
 
             scene_objects.append(tl)
 
-    cmd = "#excitation_file"
-    if multicmds[cmd] is not None:
+    cmdname = "#excitation_file"
+    if multicmds[cmdname] is not None:
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) not in [1, 3]:
-                logger.exception(f"{cmd} requires either one or three parameter(s)")
+                logger.exception(f"{cmdname} requires either one or three parameter(s)")
                 raise ValueError
 
             if len(tmp) > 1:
@@ -193,12 +223,23 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) != 3 and len(tmp) < 5:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " has an incorrect number of parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " has an incorrect number of parameters"
+                )
                 raise ValueError
             if len(tmp) == 3:
                 rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])))
             else:
-                rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])), id=tmp[3], outputs=" ".join(tmp[4:]))
+                rx = Rx(
+                    p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
+                    id=tmp[3],
+                    outputs=" ".join(tmp[4:]),
+                )
 
             scene_objects.append(rx)
 
@@ -207,7 +248,9 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) != 9:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly nine parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly nine parameters"
+                )
                 raise ValueError
 
             p1 = (float(tmp[0]), float(tmp[1]), float(tmp[2]))
@@ -222,7 +265,14 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) != 11:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly eleven parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires exactly eleven parameters"
+                )
                 raise ValueError
 
             p1 = (float(tmp[0]), float(tmp[1]), float(tmp[2]))
@@ -245,10 +295,14 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) != 5:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly five parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly five parameters"
+                )
                 raise ValueError
 
-            material = Material(er=float(tmp[0]), se=float(tmp[1]), mr=float(tmp[2]), sm=float(tmp[3]), id=tmp[4])
+            material = Material(
+                er=float(tmp[0]), se=float(tmp[1]), mr=float(tmp[2]), sm=float(tmp[3]), id=tmp[4]
+            )
             scene_objects.append(material)
 
     cmdname = "#add_dispersion_debye"
@@ -257,7 +311,14 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
 
             if len(tmp) < 4:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least four parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires at least four parameters"
+                )
                 raise ValueError
 
             poles = int(tmp[0])
@@ -269,7 +330,9 @@ def process_multicmds(multicmds):
                 er_delta.append(float(tmp[pole]))
                 tau.append(float(tmp[pole + 1]))
 
-            debye_dispersion = AddDebyeDispersion(poles=poles, er_delta=er_delta, tau=tau, material_ids=material_ids)
+            debye_dispersion = AddDebyeDispersion(
+                poles=poles, er_delta=er_delta, tau=tau, material_ids=material_ids
+            )
             scene_objects.append(debye_dispersion)
 
     cmdname = "#add_dispersion_lorentz"
@@ -278,7 +341,14 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
 
             if len(tmp) < 5:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least five parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires at least five parameters"
+                )
                 raise ValueError
 
             poles = int(tmp[0])
@@ -303,7 +373,14 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
 
             if len(tmp) < 5:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least five parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires at least five parameters"
+                )
                 raise ValueError
 
             poles = int(tmp[0])
@@ -315,7 +392,9 @@ def process_multicmds(multicmds):
                 tau.append(float(tmp[pole]))
                 alpha.append(float(tmp[pole + 1]))
 
-            drude_dispersion = AddDrudeDispersion(poles=poles, material_ids=material_ids, tau=tau, alpha=alpha)
+            drude_dispersion = AddDrudeDispersion(
+                poles=poles, material_ids=material_ids, tau=tau, alpha=alpha
+            )
             scene_objects.append(drude_dispersion)
 
     cmdname = "#soil_peplinski"
@@ -324,7 +403,14 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
 
             if len(tmp) != 7:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at exactly seven parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires at exactly seven parameters"
+                )
                 raise ValueError
             soil = SoilPeplinski(
                 sand_fraction=float(tmp[0]),
@@ -342,7 +428,14 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) != 11:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly eleven parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires exactly eleven parameters"
+                )
                 raise ValueError
 
             p1 = float(tmp[0]), float(tmp[1]), float(tmp[2])
@@ -357,7 +450,14 @@ def process_multicmds(multicmds):
         for cmdinstance in multicmds[cmdname]:
             tmp = cmdinstance.split()
             if len(tmp) != 7:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly seven parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires exactly seven parameters"
+                )
                 raise ValueError
 
             p1 = float(tmp[0]), float(tmp[1]), float(tmp[2])
@@ -371,7 +471,14 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
 
             if len(tmp) != 9:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at exactly nine parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires at exactly nine parameters"
+                )
                 raise ValueError
             material_range = MaterialRange(
                 er_lower=float(tmp[0]),
@@ -392,7 +499,9 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
 
             if len(tmp) < 2:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least two parameters")
+                logger.exception(
+                    "'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least two parameters"
+                )
                 raise ValueError
 
             tokens = len(tmp)
@@ -409,7 +518,14 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
 
             if len(tmp) != 12:
-                logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly twelve parameters")
+                logger.exception(
+                    "'"
+                    + cmdname
+                    + ": "
+                    + " ".join(tmp)
+                    + "'"
+                    + " requires exactly twelve parameters"
+                )
                 raise ValueError
 
             pml_cfs = PMLCFS(

reframe_tests/base_tests.py

@@ -212,5 +212,9 @@ class GprMaxRegressionTest(GprMaxBaseTest):
             )
 
 
+class GprMaxAPIRegressionTest(GprMaxRegressionTest):
+    executable = "time -p python"
+
+
 class GprMaxMpiTest(GprMaxBaseTest):
     pass

reframe_tests/reframe_tests.py

@@ -1,5 +1,5 @@
 import reframe as rfm
-from base_tests import GprMaxRegressionTest
+from base_tests import GprMaxAPIRegressionTest, GprMaxRegressionTest
 from reframe.core.builtins import parameter, run_after
 
 """ReFrame tests for basic functionality
@@ -82,3 +82,71 @@ class BasicModelsTest(GprMaxRegressionTest):
         self.executable_opts = [self.input_file, "-o", self.output_file]
         self.postrun_cmds = [f"python -m toolboxes.Plotting.plot_Ascan -save {self.output_file}"]
         self.keep_files = [self.input_file, self.output_file, f"{self.model}.pdf"]
+
+
+@rfm.simple_test
+class AntennaModelsTest(GprMaxRegressionTest):
+    tags = {"test", "serial", "regression", "antenna"}
+
+    # List of available antenna test models
+    model = parameter(
+        [
+            "antenna_wire_dipole_fs",
+        ]
+    )
+    num_cpus_per_task = 16
+
+    @run_after("init")
+    def set_filenames(self):
+        self.input_file = f"{self.model}.in"
+        self.output_file = f"{self.model}.h5"
+        self.executable_opts = [self.input_file, "-o", self.output_file]
+        self.postrun_cmds = [f"python -m toolboxes.Plotting.plot_Ascan -save {self.output_file}"]
+        self.postrun_cmds = [
+            f"python -m toolboxes.Plotting.plot_antenna_params -save {self.output_file}"
+        ]
+
+        antenna_t1_params = f"{self.model}_t1_params.pdf"
+        antenna_ant_params = f"{self.model}_ant_params.pdf"
+        plot_ascan_output = f"{self.model}.pdf"
+        geometry_view = f"{self.model}.vtu"
+        self.keep_files = [
+            self.input_file,
+            self.output_file,
+            antenna_t1_params,
+            antenna_ant_params,
+            plot_ascan_output,
+            geometry_view,
+        ]
+
+
+@rfm.simple_test
+class SubgridTest(GprMaxAPIRegressionTest):
+    tags = {"test", "api", "serial", "regression", "subgrid"}
+
+    # List of available subgrid test models
+    model = parameter(
+        [
+            "cylinder_fs",
+            # "gssi_400_over_fractal_subsurface",  # Takes ~1hr 30m on ARCHER2
+        ]
+    )
+    num_cpus_per_task = 16
+
+    @run_after("init")
+    def set_filenames(self):
+        self.input_file = f"{self.model}.py"
+        self.output_file = f"{self.model}.h5"
+        self.executable_opts = [self.input_file, "-o", self.output_file]
+        self.postrun_cmds = [f"python -m toolboxes.Plotting.plot_Ascan -save {self.output_file}"]
+
+        geometry_view = f"{self.model}.vti"
+        subgrid_geometry_view = f"{self.model}_sg.vti"
+        plot_ascan_output = f"{self.model}.pdf"
+        self.keep_files = [
+            self.input_file,
+            self.output_file,
+            geometry_view,
+            subgrid_geometry_view,
+            plot_ascan_output,
+        ]

reframe_tests/src/antenna_wire_dipole_fs.in

@@ -0,0 +1,15 @@
+#title: Wire antenna - half-wavelength dipole in free-space
+#domain: 0.050 0.050 0.200
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 60e-9
+
+#waveform: gaussian 1 1e9 mypulse
+#transmission_line: z 0.025 0.025 0.100 73 mypulse
+
+## 150mm length
+#edge: 0.025 0.025 0.025 0.025 0.025 0.175 pec
+
+## 1mm gap at centre of dipole
+#edge: 0.025 0.025 0.100 0.025 0.025 0.101 free_space
+
+#geometry_view: 0.020 0.020 0.020 0.030 0.030 0.180 0.001 0.001 0.001 antenna_wire_dipole_fs f

reframe_tests/src/cylinder_fs.py

@@ -0,0 +1,115 @@
+"""Cylinder in freespace
+
+This example model demonstrates how to use subgrids at a basic level.
+
+The geometry is 3D (required for any use of subgrids) and is of a water-filled
+cylindrical object in freespace. The subgrid encloses the cylinderical object
+using a fine spatial discretisation (1mm), and a courser spatial discretisation
+(5mm) is used in the rest of the model (main grid). A simple Hertzian dipole
+source is used with a waveform shaped as the first derivative of a gaussian.
+"""
+
+from pathlib import Path
+
+import gprMax
+from gprMax.materials import calculate_water_properties
+
+# File path - used later to specify name of output files
+fn = Path(__file__)
+parts = fn.parts
+
+# Subgrid spatial discretisation in x, y, z directions
+dl_sg = 1e-3
+
+# Subgrid ratio - must always be an odd integer multiple
+ratio = 5
+dl = dl_sg * ratio
+
+# Domain extent
+x = 0.500
+y = 0.500
+z = 0.500
+
+# Time window
+tw = 6e-9
+
+scene = gprMax.Scene()
+
+title = gprMax.Title(name=fn.name)
+dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+domain = gprMax.Domain(p1=(x, y, z))
+time_window = gprMax.TimeWindow(time=tw)
+
+wf = gprMax.Waveform(wave_type="gaussiandot", amp=1, freq=1.5e9, id="mypulse")
+hd = gprMax.HertzianDipole(polarisation="z", p1=(0.205, 0.400, 0.250), waveform_id="mypulse")
+rx = gprMax.Rx(p1=(0.245, 0.400, 0.250))
+
+scene.add(title)
+scene.add(dxdydz)
+scene.add(domain)
+scene.add(time_window)
+scene.add(wf)
+scene.add(hd)
+scene.add(rx)
+
+# Cylinder parameters
+c1 = (0.225, 0.250, 0.100)
+c2 = (0.225, 0.250, 0.400)
+r = 0.010
+sg1 = (c1[0] - r, c1[1] - r, c1[2])
+sg2 = (c2[0] + r, c2[1] + r, c2[2])
+
+# Create subgrid
+subgrid = gprMax.SubGridHSG(p1=sg1, p2=sg2, ratio=ratio, id="sg")
+scene.add(subgrid)
+
+# Create water material
+eri, er, tau, sig = calculate_water_properties()
+water = gprMax.Material(er=eri, se=sig, mr=1, sm=0, id="water")
+subgrid.add(water)
+water = gprMax.AddDebyeDispersion(poles=1, er_delta=[er - eri], tau=[tau], material_ids=["water"])
+subgrid.add(water)
+
+# Add cylinder to subgrid
+cylinder = gprMax.Cylinder(p1=c1, p2=c2, r=r, material_id="water")
+subgrid.add(cylinder)
+
+# Create some geometry views for both subgrid and main grid
+gvsg = gprMax.GeometryView(
+    p1=sg1,
+    p2=sg2,
+    dl=(dl_sg, dl_sg, dl_sg),
+    filename=fn.with_suffix("").parts[-1] + "_sg",
+    output_type="n",
+)
+subgrid.add(gvsg)
+
+gv1 = gprMax.GeometryView(
+    p1=(0, 0, 0),
+    p2=(x, y, z),
+    dl=(dl, dl, dl),
+    filename=fn.with_suffix("").parts[-1],
+    output_type="n",
+)
+scene.add(gv1)
+
+# Create some snapshots of entire domain
+for i in range(5):
+    s = gprMax.Snapshot(
+        p1=(0, 0, 0),
+        p2=(x, y, z),
+        dl=(dl, dl, dl),
+        time=(i + 0.5) * 1e-9,
+        filename=fn.with_suffix("").parts[-1] + "_" + str(i + 1),
+    )
+    scene.add(s)
+
+gprMax.run(
+    scenes=[scene],
+    n=1,
+    geometry_only=False,
+    outputfile=fn,
+    subgrid=True,
+    autotranslate=True,
+    log_level=25,
+)

reframe_tests/src/gssi_400_over_fractal_subsurface.py

@@ -0,0 +1,166 @@
+"""GPR antenna model (like a GSSI 400MHz antenna) over layered media with a
+    rough subsurface interface.
+
+This example model demonstrates how to use subgrids at a more advanced level -
+    combining use of an imported antenna model and rough subsurface interface.
+
+The geometry is 3D (required for any use of subgrids) and is of a 2 layered
+subsurface. The top layer in a sandy soil and the bottom layer a soil with
+higher permittivity (both have some simple conductive loss). There is a rough
+interface between the soil layers. A GPR antenna model (like a GSSI 400MHz
+antenna) is imported and placed on the surface of the layered media. The antenna
+is meshed using a subgrid with a fine spatial discretisation (1mm), and a
+courser spatial discretisation (9mm) is used in the rest of the model (main
+grid).
+"""
+
+from pathlib import Path
+
+import numpy as np
+
+import gprMax
+from toolboxes.GPRAntennaModels.GSSI import antenna_like_GSSI_400
+
+# File path - used later to specify name of output files
+fn = Path(__file__)
+parts = fn.parts
+
+# Subgrid spatial discretisation in x, y, z directions
+dl_sg = 1e-3
+
+# Subgrid ratio - must always be an odd integer multiple
+ratio = 9
+dl = dl_sg * ratio
+
+# Domain extent
+x = 3
+y = 1
+z = 2
+
+# Time window
+# Estimated two way travel time over 1 metre in material with highest
+# permittivity, slowest velocity.
+tw = 2 / 3e8 * (np.sqrt(3.2) + np.sqrt(9))
+
+scene = gprMax.Scene()
+
+title = gprMax.Title(name=fn.name)
+dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+domain = gprMax.Domain(p1=(x, y, z))
+time_window = gprMax.TimeWindow(time=tw)
+
+scene.add(title)
+scene.add(dxdydz)
+scene.add(domain)
+scene.add(time_window)
+
+# Dimensions of antenna case
+antenna_case = (0.3, 0.3, 0.178)
+
+# Position of antenna
+antenna_p = (x / 2, y / 2, 170 * dl)
+
+# Extra distance surrounding antenna for subgrid
+bounding_box = 2 * dl
+
+# Subgrid extent
+sg_x0 = antenna_p[0] - antenna_case[0] / 2 - bounding_box
+sg_y0 = antenna_p[1] - antenna_case[1] / 2 - bounding_box
+sg_z0 = antenna_p[2] - bounding_box
+sg_x1 = antenna_p[0] + antenna_case[0] / 2 + bounding_box
+sg_y1 = antenna_p[1] + antenna_case[1] / 2 + bounding_box
+sg_z1 = antenna_p[2] + antenna_case[2] + bounding_box
+
+# Create subgrid
+sg = gprMax.SubGridHSG(p1=[sg_x0, sg_y0, sg_z0], p2=[sg_x1, sg_y1, sg_z1], ratio=ratio, id="sg")
+scene.add(sg)
+
+# Create and add a box of homogeneous material to main grid - sandy_soil
+sandy_soil = gprMax.Material(er=3.2, se=0.397e-3, mr=1, sm=0, id="sandy_soil")
+scene.add(sandy_soil)
+b1 = gprMax.Box(p1=(0, 0, 0), p2=(x, y, antenna_p[2]), material_id="sandy_soil")
+scene.add(b1)
+
+# Position box of sandy_soil in the subgrid.
+#   It has to be positioned manually because it traverses the main grid/subgrid
+#   interface. Grid traversal is when objects extend beyond the outer surface.
+#   Setting autotranslate to false allows you to place objects beyond the outer
+#   surface.
+
+# PML separation from the outer surface
+ps = ratio // 2 + 2
+# Number of PML cells in the subgrid
+pc = 6
+# Inner surface/outer surface separation
+isos = 3 * ratio
+
+# Calculate maximum z-coordinate (height) for box of sandy_soil in subgrid
+h = antenna_p[2] - sg_z0 + (ps + pc + isos) * dl_sg
+
+# Create and add a box of homogeneous material to subgrid - sandy_soil
+sg.add(sandy_soil)
+b2 = gprMax.Box(p1=(0, 0, 0), p2=(411 * dl_sg, 411 * dl_sg, h), material_id="sandy_soil")
+# Set autotranslate for the box object to false
+b2.autotranslate = False
+sg.add(b2)
+
+# Import antenna model and add components to subgrid
+gssi_objects = antenna_like_GSSI_400(*antenna_p, resolution=dl_sg)
+for obj in gssi_objects:
+    sg.add(obj)
+
+# Create and add a homogeneous material with a rough surface
+soil = gprMax.Material(er=9, se=0.397e-3, mr=1, sm=0, id="soil")
+scene.add(soil)
+
+fb = gprMax.FractalBox(
+    p1=(0, 0, 0),
+    p2=(3, 1, 1),
+    frac_dim=1.5,
+    weighting=(1, 1, 1),
+    n_materials=1,
+    mixing_model_id="soil",
+    id="fbox",
+    seed=1,
+)
+scene.add(fb)
+
+rough_surf = gprMax.AddSurfaceRoughness(
+    p1=(0, 0, 1),
+    p2=(3, 1, 1),
+    frac_dim=1.5,
+    weighting=(1, 1),
+    limits=(0.4, 1.2),
+    fractal_box_id="fbox",
+    seed=1,
+)
+scene.add(rough_surf)
+
+# Create some snapshots and geometry views
+for i in range(1, 51):
+    snap = gprMax.Snapshot(
+        p1=(0, y / 2, 0),
+        p2=(x, y / 2 + dl, z),
+        dl=(dl, dl, dl),
+        filename=Path(*parts[:-1], f"{parts[-1]}_{str(i)}").name,
+        time=i * tw / 50,
+    )
+    scene.add(snap)
+
+gvsg = gprMax.GeometryView(
+    p1=(sg_x0, sg_y0, sg_z0),
+    p2=(sg_x1, sg_y1, sg_z1),
+    dl=(dl_sg, dl_sg, dl_sg),
+    filename=fn.with_suffix("").parts[-1] + "_sg",
+    output_type="n",
+)
+sg.add(gvsg)
+
+gv1 = gprMax.GeometryView(
+    p1=(0, 0, 0), p2=domain.props.p1, dl=dl, filename=fn.with_suffix("").parts[-1], output_type="n"
+)
+scene.add(gv1)
+
+gprMax.run(
+    scenes=[scene], n=1, geometry_only=False, outputfile=fn, subgrid=True, autotranslate=True
+)

toolboxes/GPRAntennaModels/GSSI.py

@@ -6,6 +6,7 @@
 # Please use the attribution at http://dx.doi.org/10.1190/1.3548506
 
 import logging
+from pathlib import Path
 
 import gprMax
 
@@ -64,7 +65,9 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
         patchheight = 0.016
         tx = x + 0.114, y + 0.052, z + skidthickness
     else:
-        logger.exception("This antenna module can only be used with a spatial discretisation of 1mm or 2mm")
+        logger.exception(
+            "This antenna module can only be used with a spatial discretisation of 1mm or 2mm"
+        )
         raise ValueError
 
     # If using parameters from an optimisation
@@ -98,7 +101,9 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
             sourceresistance = 230  # Correction for old (< 123) GprMax3D bug (optimised to 4)
             rxres = 925  # Resistance at Rx bowtie
             absorber1 = gprMax.Material(er=1.58, se=0.428, mr=1, sm=0, id="absorber1")
-            absorber2 = gprMax.Material(er=3, se=0, mr=1, sm=0, id="absorber2")  # Foam modelled as PCB material
+            absorber2 = gprMax.Material(
+                er=3, se=0, mr=1, sm=0, id="absorber2"
+            )  # Foam modelled as PCB material
             pcb = gprMax.Material(er=3, se=0, mr=1, sm=0, id="pcb")
             hdpe = gprMax.Material(er=2.35, se=0, mr=1, sm=0, id="hdpe")
             rxres = gprMax.Material(er=3, se=(1 / rxres) * (dy / (dx * dz)), mr=1, sm=0, id="rxres")
@@ -117,7 +122,9 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
                 tau=[1.00723e-11, 1.55686e-10, 3.44129e-10],
                 material_ids=["absorber1"],
             )
-            absorber2 = gprMax.Material(er=3, se=0, mr=1, sm=0, id="absorber2")  # Foam modelled as PCB material
+            absorber2 = gprMax.Material(
+                er=3, se=0, mr=1, sm=0, id="absorber2"
+            )  # Foam modelled as PCB material
             pcb = gprMax.Material(er=3, se=0, mr=1, sm=0, id="pcb")
             hdpe = gprMax.Material(er=2.35, se=0, mr=1, sm=0, id="hdpe")
             rxres = gprMax.Material(er=3, se=(1 / rxres) * (dy / (dx * dz)), mr=1, sm=0, id="rxres")
@@ -299,12 +306,20 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
         scene_objects.extend((p9, p10))
 
         # Edges that represent wire between bowtie halves in 1mm model
-        e1 = gprMax.Edge(p1=(tx[0] - 0.059, tx[1] - dy, tx[2]), p2=(tx[0] - 0.059, tx[1], tx[2]), material_id="pec")
+        e1 = gprMax.Edge(
+            p1=(tx[0] - 0.059, tx[1] - dy, tx[2]),
+            p2=(tx[0] - 0.059, tx[1], tx[2]),
+            material_id="pec",
+        )
         e2 = gprMax.Edge(
-            p1=(tx[0] - 0.059, tx[1] + dy, tx[2]), p2=(tx[0] - 0.059, tx[1] + 0.002, tx[2]), material_id="pec"
+            p1=(tx[0] - 0.059, tx[1] + dy, tx[2]),
+            p2=(tx[0] - 0.059, tx[1] + 0.002, 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] + dz, tx[2]), p2=(tx[0], tx[1] + 0.002, tx[2]), material_id="pec")
+        e4 = gprMax.Edge(
+            p1=(tx[0], tx[1] + dz, tx[2]), p2=(tx[0], tx[1] + 0.002, tx[2]), material_id="pec"
+        )
         scene_objects.extend((e1, e2, e3, e4))
 
     elif resolution == 0.002:
@@ -359,7 +374,9 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
     scene_objects.extend((p11, p12))
 
     # Skid
-    b10 = gprMax.Box(p1=(x, y, z), p2=(x + casesize[0], y + casesize[1], z + skidthickness), material_id="hdpe")
+    b10 = gprMax.Box(
+        p1=(x, y, z), p2=(x + casesize[0], y + casesize[1], z + skidthickness), material_id="hdpe"
+    )
     scene_objects.append(b10)
 
     # Geometry views
@@ -384,17 +401,25 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
         # Gaussian pulse
         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"
+            polarisation="y",
+            p1=(tx[0], tx[1], tx[2]),
+            resistance=sourceresistance,
+            waveform_id="my_gaussian",
         )
         scene_objects.extend((w1, vs1))
 
     elif optstate == "GiannakisPaper":
         # Optimised custom pulse
         exc1 = gprMax.ExcitationFile(
-            filepath="toolboxes/GPRAntennaModels/GSSI_1500MHz_pulse.txt", kind="linear", fill_value="extrapolate"
+            filepath=Path(__file__).with_name("GSSI_1500MHz_pulse.txt"),
+            kind="linear",
+            fill_value="extrapolate",
         )
         vs1 = gprMax.VoltageSource(
-            polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_pulse"
+            polarisation="y",
+            p1=(tx[0], tx[1], tx[2]),
+            resistance=sourceresistance,
+            waveform_id="my_pulse",
         )
         scene_objects.extend((exc1, vs1))
 
@@ -402,7 +427,9 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
     if resolution == 0.001:
         if optstate == "WarrenThesis" or optstate == "DebyeAbsorber":
             e1 = gprMax.Edge(
-                p1=(tx[0] - 0.059, tx[1], tx[2]), p2=(tx[0] - 0.059, tx[1] + dy, tx[2]), material_id="rxres"
+                p1=(tx[0] - 0.059, tx[1], tx[2]),
+                p2=(tx[0] - 0.059, tx[1] + dy, tx[2]),
+                material_id="rxres",
             )
             scene_objects.append(e1)
         r1 = gprMax.Rx(p1=(tx[0] - 0.059, tx[1], tx[2]), id="rxbowtie", outputs="Ey")
@@ -411,7 +438,9 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
     elif resolution == 0.002:
         if optstate == "WarrenThesis" or optstate == "DebyeAbsorber":
             e1 = gprMax.Edge(
-                p1=(tx[0] - 0.060, tx[1], tx[2]), p2=(tx[0] - 0.060, tx[1] + dy, tx[2]), material_id="rxres"
+                p1=(tx[0] - 0.060, tx[1], tx[2]),
+                p2=(tx[0] - 0.060, tx[1] + dy, tx[2]),
+                material_id="rxres",
             )
             scene_objects.append(e1)
         r1 = gprMax.Rx(p1=(tx[0] - 0.060, tx[1], tx[2]), id="rxbowtie", outputs="Ey")
@@ -492,7 +521,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     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
+    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")
@@ -510,13 +543,21 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
         )
         b2 = gprMax.Box(
             p1=(x + casethickness, y + casethickness, z + skidthickness - 0.002),
-            p2=(x + casesize[0] - casethickness, y + casesize[1] - casethickness, z + casesize[2] - casethickness),
+            p2=(
+                x + casesize[0] - casethickness,
+                y + casesize[1] - casethickness,
+                z + casesize[2] - casethickness,
+            ),
             material_id="free_space",
         )
 
         # Metallic enclosure
         b3 = gprMax.Box(
-            p1=(x + casethickness, y + casethickness, z + skidthickness + (metalmiddleplateheight - metalboxheight)),
+            p1=(
+                x + casethickness,
+                y + casethickness,
+                z + skidthickness + (metalmiddleplateheight - metalboxheight),
+            ),
             p2=(
                 x + casesize[0] - casethickness,
                 y + casesize[1] - casethickness,
@@ -552,7 +593,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
 
         # PCB
         b6 = gprMax.Box(
-            p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+            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,
@@ -561,7 +606,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
             material_id="pcb",
         )
         b7 = gprMax.Box(
-            p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+            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,
@@ -581,14 +630,22 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
         )
         b9 = gprMax.Box(
             p1=(x + casethickness, y + casethickness, z + skidthickness - 0.002),
-            p2=(x + casesize[0] - casethickness, y + casesize[1] - casethickness, z + casesize[2] - casethickness),
+            p2=(
+                x + casesize[0] - casethickness,
+                y + casesize[1] - casethickness,
+                z + casesize[2] - casethickness,
+            ),
             material_id="free_space",
             averaging="n",
         )
 
         # Metallic enclosure
         b10 = gprMax.Box(
-            p1=(x + casethickness, y + casethickness, z + skidthickness + (metalmiddleplateheight - metalboxheight)),
+            p1=(
+                x + casethickness,
+                y + casethickness,
+                z + skidthickness + (metalmiddleplateheight - metalboxheight),
+            ),
             p2=(
                 x + casesize[0] - casethickness,
                 y + casesize[1] - casethickness,
@@ -626,7 +683,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
 
         # PCB
         b13 = gprMax.Box(
-            p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
+            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,
@@ -672,7 +733,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     # 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),
+        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,
@@ -683,7 +748,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     )
     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),
+        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,
@@ -714,7 +783,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     # 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),
+        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,
@@ -725,7 +798,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     )
     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),
+        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,
@@ -736,12 +813,18 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     )
 
     # 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")
+    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"
+        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")
+    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
@@ -758,7 +841,9 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     # Skid
     if smooth_dec == "yes":
         b16 = gprMax.Box(
-            p1=(x, y, z), p2=(x + casesize[0], y + casesize[1], z + skidthickness - 0.002), material_id="hdpe"
+            p1=(x, y, z),
+            p2=(x + casesize[0], y + casesize[1], z + skidthickness - 0.002),
+            material_id="hdpe",
         )
     elif smooth_dec == "no":
         b16 = gprMax.Box(
@@ -773,22 +858,33 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
     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"
+            polarisation="y",
+            p1=(tx[0], tx[1], tx[2]),
+            resistance=sourceresistance,
+            waveform_id="my_gaussian",
         )
         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"
+            polarisation="y",
+            p1=(tx[0], tx[1], tx[2]),
+            resistance=sourceresistance,
+            waveform_id="my_gaussian",
         )
         scene_objects.extend((w1, tl1))
     else:
         # Optimised custom pulse
         exc1 = gprMax.ExcitationFile(
-            filepath="toolboxes/GPRAntennaModels/GSSI_400MHz_pulse.txt", kind="linear", fill_value="extrapolate"
+            filepath=Path(__file__).with_name("GSSI_400MHz_pulse.txt"),
+            kind="linear",
+            fill_value="extrapolate",
         )
         vs1 = gprMax.VoltageSource(
-            polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_pulse"
+            polarisation="y",
+            p1=(tx[0], tx[1], tx[2]),
+            resistance=sourceresistance,
+            waveform_id="my_pulse",
         )
         scene_objects.extend((exc1, vs1))
 

toolboxes/Plotting/plot_antenna_params.py

@@ -28,7 +28,9 @@ import numpy as np
 logger = logging.getLogger(__name__)
 
 
-def calculate_antenna_params(filename, tltxnumber=1, tlrxnumber=None, rxnumber=None, rxcomponent=None):
+def calculate_antenna_params(
+    filename, tltxnumber=1, tlrxnumber=None, rxnumber=None, rxcomponent=None
+):
     """Calculates antenna parameters - incident, reflected and total volatges
         and currents; s11, (s21) and input impedance.
 
@@ -224,7 +226,10 @@ def mpl_plot(
 
     # Print some useful values from s11, and input impedance
     s11minfreq = np.where(s11[pltrange] == np.amin(s11[pltrange]))[0][0]
-    logger.info(f"s11 minimum: {np.amin(s11[pltrange]):g} dB at " + f"{freqs[s11minfreq + pltrangemin]:g} Hz")
+    logger.info(
+        f"s11 minimum: {np.amin(s11[pltrange]):g} dB at "
+        + f"{freqs[s11minfreq + pltrangemin]:g} Hz"
+    )
     logger.info(f"At {freqs[s11minfreq + pltrangemin]:g} Hz...")
     logger.info(
         f"Input impedance: {np.abs(zin[s11minfreq + pltrangemin]):.1f}"
@@ -236,7 +241,10 @@ def mpl_plot(
     # Figure 1
     # Plot incident voltage
     fig1, ax = plt.subplots(
-        num="Transmitter transmission line parameters", figsize=(20, 12), facecolor="w", edgecolor="w"
+        num="Transmitter transmission line parameters",
+        figsize=(20, 12),
+        facecolor="w",
+        edgecolor="w",
     )
     gs1 = gridspec.GridSpec(4, 2, hspace=0.7)
     ax = plt.subplot(gs1[0, 0])
@@ -368,7 +376,9 @@ def mpl_plot(
 
     # Figure 2
     # Plot frequency spectra of s11
-    fig2, ax = plt.subplots(num="Antenna parameters", figsize=(20, 12), facecolor="w", edgecolor="w")
+    fig2, ax = plt.subplots(
+        num="Antenna parameters", figsize=(20, 12), facecolor="w", edgecolor="w"
+    )
     gs2 = gridspec.GridSpec(2, 2, hspace=0.3)
     ax = plt.subplot(gs2[0, 0])
     markerline, stemlines, baseline = ax.stem(freqs[pltrange], s11[pltrange], "-.")
@@ -458,13 +468,26 @@ def mpl_plot(
     # ax.grid(which='both', axis='both', linestyle='-.')
 
     if save:
-        savename1 = filename.stem + "_tl_params"
-        savename1 = filename.parent / savename1
-        savename2 = filename.stem + "_ant_params"
-        savename2 = filename.parent / savename2
+        filepath = Path(filename)
+        savename1 = filepath.stem + "_tl_params"
+        savename1 = filepath.parent / savename1
+        savename2 = filepath.stem + "_ant_params"
+        savename2 = filepath.parent / savename2
         # Save a PDF of the figure
-        fig1.savefig(savename1.with_suffix(".pdf"), dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1)
-        fig2.savefig(savename2.with_suffix(".pdf"), dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1)
+        fig1.savefig(
+            savename1.with_suffix(".pdf"),
+            dpi=None,
+            format="pdf",
+            bbox_inches="tight",
+            pad_inches=0.1,
+        )
+        fig2.savefig(
+            savename2.with_suffix(".pdf"),
+            dpi=None,
+            format="pdf",
+            bbox_inches="tight",
+            pad_inches=0.1,
+        )
         # Save a PNG of the figure
         # fig1.savefig(savename1.with_suffix('.png'), dpi=150, format='png',
         #             bbox_inches='tight', pad_inches=0.1)
@@ -485,7 +508,9 @@ if __name__ == "__main__":
         usage="cd gprMax; python -m toolboxes.Plotting.plot_antenna_params outputfile",
     )
     parser.add_argument("outputfile", help="name of output file including path")
-    parser.add_argument("--tltx-num", default=1, type=int, help="transmitter antenna - transmission line number")
+    parser.add_argument(
+        "--tltx-num", default=1, type=int, help="transmitter antenna - transmission line number"
+    )
     parser.add_argument("--tlrx-num", type=int, help="receiver antenna - transmission line number")
     parser.add_argument("--rx-num", type=int, help="receiver antenna - output number")
     parser.add_argument(
@@ -495,7 +520,10 @@ if __name__ == "__main__":
         choices=["Ex", "Ey", "Ez"],
     )
     parser.add_argument(
-        "-save", action="store_true", default=False, help="save plot directly to file, i.e. do not display"
+        "-save",
+        action="store_true",
+        default=False,
+        help="save plot directly to file, i.e. do not display",
     )
     args = parser.parse_args()