Updates to get antenna models working with revised antenna library.

2025-08-07 15:10:13 +08:00 · 2019-10-31 14:27:22 +00:00
--- a/.gitignore
+++ b/.gitignore
@@ -21,3 +21,8 @@ dist/

 # Jupyter notebook checkpoints
 .ipynb_checkpoints/
+
+# Craig's files
+user_models/GPU-paper-landmines/
+user_models/IWAGPR2017-challenge/
+user_models/to-update/
--- a/gprMax/init.py
+++ b/gprMax/init.py
@@ -49,6 +49,7 @@ from .cmds_geometry.fractal_box import FractalBox
 from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness
 from .cmds_geometry.add_surface_water import AddSurfaceWater
 from .cmds_geometry.add_grass import AddGrass
+from .cmds_geometry.geometry_objects_read import GeometryObjectsRead

 from .hash_cmds_file import user_libs_fn_to_scene_obj

--- a/gprMax/cmds_geometry/geometry_objects_read.py
+++ b/gprMax/cmds_geometry/geometry_objects_read.py
@@ -17,11 +17,17 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.

 import logging
-import os
+from pathlib import Path
+
+import h5py

 import gprMax.config as config
 from .cmds_geometry import UserObjectGeometry
+from ..cython.geometry_primitives import build_voxels_from_array
 from ..exceptions import CmdInputError
+from ..hash_cmds_file import check_cmd_names
+from ..hash_cmds_multiuse import process_multicmds
+from ..utilities import round_value

 log = logging.getLogger(__name__)

@@ -31,6 +37,7 @@ class GeometryObjectsRead(UserObjectGeometry):
    log.debug('More work required here.')

    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
        self.order = 1
        self.hash = '#geometry_objects_read'

@@ -49,10 +56,12 @@ class GeometryObjectsRead(UserObjectGeometry):
        xs, ys, zs = uip.discretise_point(p1)

        # See if material file exists at specified path and if not try input file directory
-        if not os.path.isfile(matfile):
-            matfile = os.path.abspath(os.path.join(G.inputdirectory, matfile))
+        matfile = Path(matfile)

-        matstr = os.path.splitext(os.path.split(matfile)[1])[0]
+        if not matfile.exists():
+            matfile = Path(config.sim_config.input_file_path.parent, matfile)
+
+        matstr = matfile.with_suffix('').name
        numexistmaterials = len(G.materials)

        # Read materials from file
@@ -64,7 +73,7 @@ class GeometryObjectsRead(UserObjectGeometry):
        singlecmdsimport, multicmdsimport, geometryimport = check_cmd_names(materials, checkessential=False)

        # Process parameters for commands that can occur multiple times in the model
-        process_multicmds(multicmdsimport, G)
+        process_multicmds(multicmdsimport)

        # Update material type
        for material in G.materials:
@@ -75,8 +84,9 @@ class GeometryObjectsRead(UserObjectGeometry):
                    material.type = 'imported'

        # See if geometry object file exists at specified path and if not try input file directory
-        if not os.path.isfile(geofile):
-            geofile = os.path.abspath(os.path.join(G.inputdirectory, geofile))
+        geofile = Path(geofile)
+        if not geofile.exists():
+            geofile = Path(config.sim_config.input_file_path.parent, geofile)

        # Open geometry object file and read/check spatial resolution attribute
        f = h5py.File(geofile, 'r')
--- a/gprMax/cmds_multiple.py
+++ b/gprMax/cmds_multiple.py
@@ -17,7 +17,6 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.

 import logging
-import sys

 import numpy as np

@@ -415,7 +414,7 @@ class TransmissionLine(UserObjectMulti):
            raise CmdInputError(f"'{self.params_str()}' requires at least six parameters")

        # Warn about using a transmission line on GPU
-        if grid.gpu is not None:
+        if config.sim_config.general['cuda']:
            raise CmdInputError(f"'{self.params_str()}' A #transmission_line cannot currently be used with GPU solving. Consider using a #voltage_source instead.")

        # Check polarity & position parameters
@@ -430,12 +429,12 @@ class TransmissionLine(UserObjectMulti):

        xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str())

-        if resistance <= 0 or resistance >= z0:
+        if resistance <= 0 or resistance >= config.sim_config.em_consts['z0']:
            raise CmdInputError(f"'{self.params_str()}' requires a resistance greater than zero and less than the impedance of free space, i.e. 376.73 Ohms")

        # Check if there is a waveformID in the waveforms list
        if not any(x.ID == waveform_id for x in grid.waveforms):
-            raise CmdInputError(f"'{self.params_str()}' there is no waveform with the identifier {tmp[5]}")
+            raise CmdInputError(f"'{self.params_str()}' there is no waveform with the identifier {waveform_id}")

        t = TransmissionLineUser(grid)
        t.polarisation = polarisation
--- a/gprMax/cmds_single_use.py
+++ b/gprMax/cmds_single_use.py
@@ -18,7 +18,7 @@

 import inspect
 import logging
-import os
+from pathlib import Path

 from colorama import init
 from colorama import Fore
@@ -427,6 +427,10 @@ class ExcitationFile(UserObjectSingle):
    :type fill_value: float, optional
    """

+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.order = 11
+
    def create(self, G, uip):
        try:
            kwargs = dict()
@@ -442,55 +446,56 @@ class ExcitationFile(UserObjectSingle):
            except KeyError:
                raise CmdInputError('#excitation_file: requires either one or three parameter(s)')

-            # See if file exists at specified path and if not try input file directory
-            if not os.path.isfile(excitationfile):
-                excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
+        # See if file exists at specified path and if not try input file directory
+        excitationfile = Path(excitationfile)
+        if not excitationfile.exists():
+            excitationfile = Path(config.sim_config.input_file_path.parent, excitationfile)

-            log.info(f'\nExcitation file: {excitationfile}')
+        log.info(f'Excitation file: {excitationfile}')

-            # Get waveform names
-            with open(excitationfile, 'r') as f:
-                waveformIDs = f.readline().split()
+        # Get waveform names
+        with open(excitationfile, 'r') as f:
+            waveformIDs = f.readline().split()

-            # Read all waveform values into an array
-            waveformvalues = np.loadtxt(excitationfile, skiprows=1, dtype=floattype)
+        # Read all waveform values into an array
+        waveformvalues = np.loadtxt(excitationfile, skiprows=1, dtype=config.sim_config.dtypes['float_or_double'])

-            # Time array (if specified) for interpolation, otherwise use simulation time
-            if waveformIDs[0].lower() == 'time':
-                waveformIDs = waveformIDs[1:]
-                waveformtime = waveformvalues[:, 0]
-                waveformvalues = waveformvalues[:, 1:]
-                timestr = 'user-defined time array'
-            else:
-                waveformtime = np.arange(0, G.timewindow + G.dt, G.dt)
-                timestr = 'simulation time array'
+        # Time array (if specified) for interpolation, otherwise use simulation time
+        if waveformIDs[0].lower() == 'time':
+            waveformIDs = waveformIDs[1:]
+            waveformtime = waveformvalues[:, 0]
+            waveformvalues = waveformvalues[:, 1:]
+            timestr = 'user-defined time array'
+        else:
+            waveformtime = np.arange(0, G.timewindow + G.dt, G.dt)
+            timestr = 'simulation time array'

-            for waveform in range(len(waveformIDs)):
-                if any(x.ID == waveformIDs[waveform] for x in G.waveforms):
-                    raise CmdInputError(f'Waveform with ID {waveformIDs[waveform]} already exists')
-                w = Waveform()
-                w.ID = waveformIDs[waveform]
-                w.type = 'user'
+        for waveform in range(len(waveformIDs)):
+            if any(x.ID == waveformIDs[waveform] for x in G.waveforms):
+                raise CmdInputError(f'Waveform with ID {waveformIDs[waveform]} already exists')
+            w = Waveform()
+            w.ID = waveformIDs[waveform]
+            w.type = 'user'

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

-                # Truncate waveform array if it is longer than time array
-                if len(singlewaveformvalues) > len(waveformtime):
-                    singlewaveformvalues = singlewaveformvalues[:len(waveformtime)]
-                # Zero-pad end of waveform array if it is shorter than time array
-                elif len(singlewaveformvalues) < len(waveformtime):
-                    singlewaveformvalues = np.lib.pad(singlewaveformvalues,
-                                                      (0, len(singlewaveformvalues) -
-                                                      len(waveformvalues)),
-                                                      'constant', constant_values=0)
+            # Truncate waveform array if it is longer than time array
+            if len(singlewaveformvalues) > len(waveformtime):
+                singlewaveformvalues = singlewaveformvalues[:len(waveformtime)]
+            # Zero-pad end of waveform array if it is shorter than time array
+            elif len(singlewaveformvalues) < len(waveformtime):
+                singlewaveformvalues = np.lib.pad(singlewaveformvalues,
+                                                  (0, len(singlewaveformvalues) -
+                                                  len(waveformvalues)),
+                                                  'constant', constant_values=0)

-                # Interpolate waveform values
-                w.userfunc = interpolate.interp1d(waveformtime, singlewaveformvalues, **kwargs)
+            # Interpolate waveform values
+            w.userfunc = interpolate.interp1d(waveformtime, singlewaveformvalues, **kwargs)

-                log.info(f"User waveform {w.ID} created using {timestr} and, if required, interpolation parameters (kind: {kwargs['kind']}, fill value: {kwargs['fill_value']}).")
+            log.info(f"User waveform {w.ID} created using {timestr} and, if required, interpolation parameters (kind: {kwargs['kind']}, fill value: {kwargs['fill_value']}).")

-                G.waveforms.append(w)
+            G.waveforms.append(w)


 class OutputDir(UserObjectSingle):
@@ -501,7 +506,7 @@ class OutputDir(UserObjectSingle):
    """
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
-        self.order = 11
+        self.order = 12

    def create(self, grid, uip):
        config.model_configs[grid.model_num].set_output_file_path(self.kwargs['dir'])
@@ -516,7 +521,7 @@ class NumberOfModelRuns(UserObjectSingle):
    """
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
-        self.order = 12
+        self.order = 13

    def create(self, grid, uip):
        try:
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -210,7 +210,7 @@ class SimulationConfig:
            # Suppress nvcc warnings on Microsoft Windows
            if sys.platform == 'win32': self.cuda['nvcc_opts'] = '-w'
            self.get_gpus()
-            self.set_gpus()
+            self.set_single_gpu()

        # Subgrid parameter may not exist if user enters via CLI
        try:
@@ -242,7 +242,7 @@ class SimulationConfig:

        self.cuda['gpus'], self.cuda['gpus_str'] = detect_check_gpus(self.cuda['gpus'])

-    def set_gpus(self):
+    def set_single_gpu(self):
        """Adjust list of GPU object(s) to specify single GPU."""
        self.cuda['gpus'] = self.cuda['gpus'][0]
        self.cuda['gpus_str'] = self.cuda['gpus_str'][0]
@@ -304,8 +304,10 @@ class SimulationConfig:

    def set_input_file_path(self):
        """Set input file path for CLI or API."""
+        # API
        if self.args.inputfile is None:
            self.input_file_path = Path(self.args.outputfile)
+        # CLI
        else:
            self.input_file_path = Path(self.args.inputfile)

--- a/gprMax/fields_outputs.py
+++ b/gprMax/fields_outputs.py
@@ -124,7 +124,6 @@ def write_hdf5_sub_grid_outputfile(outputfile, G):
    parent = outputfile.parent

    for sg in G.subgrids:
-
        # Create an outputfile for each subgrid
        fp = stem + '_' + sg.name + suffix
        fp = parent / Path(fp)
@@ -146,6 +145,9 @@ def write_data(outputfile, G):
    Args:
        outputfile (str): Name of the output file.
        G (FDTDGrid): Parameters describing a grid in a model.
+
+    Returns:
+        f (file object): File object.
    """

    f = h5py.File(outputfile, 'w')
--- a/gprMax/geometry_outputs.py
+++ b/gprMax/geometry_outputs.py
@@ -405,9 +405,9 @@ class GeometryViewFineMultiGrid:
    """Geometry view for all grids in the simulation.

        Slicing is not supported by this class :( - only the full extent of the grids
-        are output. The subgrids are output without the non-working regions If you
-        require domainslicing GeometryView seperately for each grid you require and
-        view them at once in Paraview.
+        are output. The subgrids are output without the non-working regions.
+        If you require domain slicing, GeometryView seperately for each grid you
+        require and view them at once in Paraview.
    """

    def __init__(self, xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, G):
@@ -418,7 +418,9 @@ class GeometryViewFineMultiGrid:
            filename (str): Filename to save to.
            fileext (str): File extension of VTK file - either '.vti' for a per cell
                    geometry view, or '.vtp' for a per cell edge geometry view.
+            G (FDTDGrid): Parameters describing a grid in a model.
        """
+
        self.G = G
        self.nx = G.nx
        self.ny = G.ny
@@ -460,41 +462,36 @@ class GeometryViewFineMultiGrid:
            PolygonalData (.vtp) for a per-cell-edge geometry view.

            N.B. No Python 3 support for VTK at time of writing (03/2015)
-
-        Args:
-            G (FDTDGrid): Parameters describing a grid in a model.
        """

        G = self.G

        with open(self.filename, 'wb') as f:
-
            # refine parameters for subgrid
-
            f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
-            f.write('<VTKFile type="PolyData" version="1.0" byte_order="{}">\n'.format(config.sim_config.vtk_byteorder).encode('utf-8'))
-            f.write('<PolyData>\n<Piece NumberOfPoints="{}" NumberOfVerts="0" NumberOfLines="{}" NumberOfStrips="0" NumberOfPolys="0">\n'.format(self.vtk_numpoints, self.vtk_numlines).encode('utf-8'))
+            f.write(f"""<VTKFile type="PolyData" version="1.0" byte_order="{config.sim_config.vtk_byteorder}">\n""".encode('utf-8'))
+            f.write(f"""<PolyData>\n<Piece NumberOfPoints="{self.vtk_numpoints}" NumberOfVerts="0" NumberOfLines="{self.vtk_numlines}" NumberOfStrips="0" NumberOfPolys="0">\n""".encode('utf-8'))

            f.write('<Points>\n<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="0" />\n</Points>\n'.encode('utf-8'))
-            f.write('<Lines>\n<DataArray type="UInt32" Name="connectivity" format="appended" offset="{}" />\n'.format(self.vtk_connectivity_offset).encode('utf-8'))
-            f.write('<DataArray type="UInt32" Name="offsets" format="appended" offset="{}" />\n</Lines>\n'.format(self.vtk_offsets_offset).encode('utf-8'))
+            f.write(f"""<Lines>\n<DataArray type="UInt32" Name="connectivity" format="appended" offset="{self.vtk_connectivity_offset}" />\n""".encode('utf-8'))
+            f.write(f"""<DataArray type="UInt32" Name="offsets" format="appended" offset="{self.vtk_offsets_offset}" />\n</Lines>\n""".encode('utf-8'))
            f.write('<CellData Scalars="Material">\n'.encode('utf-8'))
-            f.write('<DataArray type="UInt32" Name="Material" format="appended" offset="{}" />\n'.format(self.vtk_materials_offset).encode('utf-8'))
+            f.write(f"""<DataArray type="UInt32" Name="Material" format="appended" offset="{self.vtk_materials_offset}" />\n""".encode('utf-8'))
            f.write('</CellData>\n'.encode('utf-8'))

            f.write('</Piece>\n</PolyData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))

            # Write points
-            log.info('writing points main grid')
+            log.info('\nWriting points main grid')
            datasize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components
            f.write(pack('I', datasize))
            for i in range(0, G.nx + 1):
                for j in range(0, G.ny + 1):
-                    for k in range(0, self.G.nz + 1):
+                    for k in range(0, G.nz + 1):
                        f.write(pack('fff', i * G.dx, j * G.dy, k * G.dz))

            for sg_v in self.sg_views:
-                log.info('writing points subgrid')
+                log.info('Writing points subgrid')
                sg_v.write_points(f, G)

            n_x_lines = self.nx * (self.ny + 1) * (self.nz + 1)
@@ -509,7 +506,7 @@ class GeometryViewFineMultiGrid:
            z_lines = np.zeros((n_z_lines, 2), dtype=np.uint32)
            z_materials = np.zeros((n_z_lines), dtype=np.uint32)

-            log.info('calculate connectivity main grid')
+            log.info('Calculate connectivity main grid')
            label = 0
            counter_x = 0
            counter_y = 0
@@ -541,7 +538,7 @@ class GeometryViewFineMultiGrid:

                        label = label + 1

-            log.info('calculate connectivity subgrids')
+            log.info('Calculate connectivity subgrids')
            for sg_v in self.sg_views:
                sg_v.populate_connectivity_and_materials(label)
                # use the last subgrids label for the next view
@@ -596,15 +593,14 @@ class GeometryViewFineMultiGrid:

        f.write('\n\n<gprMax>\n'.encode('utf-8'))
        for material in G.materials:
-            f.write('<Material name="{}">{}</Material>\n'.format(material.ID, material.numID).encode('utf-8'))
+            f.write(f"""<Material name="{material.ID}">{material.numID}</Material>\n""".encode('utf-8'))
        if not materialsonly:
            f.write('<PML name="PML boundary region">1</PML>\n'.encode('utf-8'))
            for index, src in enumerate(G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines):
-                f.write('<Sources name="{}">{}</Sources>\n'.format(src.ID, index + 2).encode('utf-8'))
+                f.write(f"""<Sources name="{src.ID}">{index + 2}</Sources>\n""".encode('utf-8'))
            for index, rx in enumerate(G.rxs):
-                f.write('<Receivers name="{}">{}</Receivers>\n'.format(rx.ID, index + 1).encode('utf-8'))
+                f.write(f"""<Receivers name="{rx.ID}">{index + 1}</Receivers>\n""".encode('utf-8'))
        f.write('</gprMax>\n'.encode('utf-8'))
-        return None


 class SubgridGeometryView:
@@ -649,7 +645,7 @@ class SubgridGeometryView:

    def populate_connectivity_and_materials(self, label):
        """Label is the starting label. 0 if no other grids are present but
-        +1 the last label used for a multigrid view.
+            +1 the last label used for a multigrid view.
        """

        sg = self.sg
--- a/gprMax/hash_cmds_file.py
+++ b/gprMax/hash_cmds_file.py
@@ -24,7 +24,7 @@ import sys

 import gprMax.config as config
 from .exceptions import CmdInputError
-from .hash_cmds_geometry import process_geometrycmds
+# from .hash_cmds_geometry import process_geometrycmds
 from .hash_cmds_multiuse import process_multicmds
 from .hash_cmds_singleuse import process_singlecmds

--- a/gprMax/sources.py
+++ b/gprMax/sources.py
@@ -317,12 +317,12 @@ class TransmissionLine(Source):
        # Cell position of where line connects to antenna/main grid
        self.antpos = 10

-        self.voltage = np.zeros(self.nl, dtype=config.dtypes['float_or_double'])
-        self.current = np.zeros(self.nl, dtype=config.dtypes['float_or_double'])
-        self.Vinc = np.zeros(G.iterations, dtype=config.dtypes['float_or_double'])
-        self.Iinc = np.zeros(G.iterations, dtype=config.dtypes['float_or_double'])
-        self.Vtotal = np.zeros(G.iterations, dtype=config.dtypes['float_or_double'])
-        self.Itotal = np.zeros(G.iterations, dtype=config.dtypes['float_or_double'])
+        self.voltage = np.zeros(self.nl, dtype=config.sim_config.dtypes['float_or_double'])
+        self.current = np.zeros(self.nl, dtype=config.sim_config.dtypes['float_or_double'])
+        self.Vinc = np.zeros(G.iterations, dtype=config.sim_config.dtypes['float_or_double'])
+        self.Iinc = np.zeros(G.iterations, dtype=config.sim_config.dtypes['float_or_double'])
+        self.Vtotal = np.zeros(G.iterations, dtype=config.sim_config.dtypes['float_or_double'])
+        self.Itotal = np.zeros(G.iterations, dtype=config.sim_config.dtypes['float_or_double'])

    def calculate_incident_V_I(self, G):
        """Calculates the incident voltage and current with a long length
--- a/user_libs/antennas/MALA.py
+++ b/user_libs/antennas/MALA.py
@@ -147,7 +147,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001):
    c8 = gprMax.Cylinder(p1=(x + 0.147, y + 0.008, z + skidthickness),
                         p2=(x + 0.147, y + 0.008,
                         z + skidthickness + casesize[2] - casethickness),
-                         r=0.007, 'free_space')
+                         r=0.007, material_id='free_space')
    b6 = gprMax.Box(p1=(x + 0.054, y + casesize[1] - 0.016, z + skidthickness),
                    p2=(x + 0.056, y + casesize[1] - 0.014,
                    z + skidthickness + casesize[2] - casethickness),
@@ -341,7 +341,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001):
                          p2=(tx[0] + 0.020, tx[1] + bowtieheight + 0.006, tx[2]),
                          material_id='txresupper')
        e16 = gprMax.Edge(p1=(tx[0] + 0.020 + dx, tx[1] + bowtieheight + 0.002, tx[2]),
-                          p2=tx[0] + 0.020 + dx, tx[1] + bowtieheight + 0.006, tx[2]),
+                          p2=(tx[0] + 0.020 + dx, tx[1] + bowtieheight + 0.006, tx[2]),
                          material_id='txresupper')
        scene_objects.extend((e11, e12, e13, e14, e15, e16))

@@ -354,7 +354,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001):
                          p2=(tx[0] - 0.023 + dx + 0.076, tx[1] - bowtieheight - dy, tx[2]),
                          material_id='rxreslower')
        e19 = gprMax.Edge(p1=(tx[0] + 0.076, tx[1] - bowtieheight - 0.004, tx[2]),
-                          p1=(tx[0] + 0.076, tx[1] - bowtieheight - dy, tx[2]),
+                          p2=(tx[0] + 0.076, tx[1] - bowtieheight - dy, tx[2]),
                          material_id='rxreslower')
        e20 = gprMax.Edge(p1=(tx[0] + dx + 0.076, tx[1] - bowtieheight - 0.004, tx[2]),
                          p2=(tx[0] + dx + 0.076, tx[1] - bowtieheight - dy, tx[2]),
--- a/user_models/antenna_like_GSSI_1500_fs.in
+++ b/user_models/antenna_like_GSSI_1500_fs.in
@@ -1,9 +0,0 @@
-#title: GSSI 1.5GHz 'like' antenna in free-space
-#domain: 0.250 0.188 0.220
-#dx_dy_dz: 0.001 0.001 0.001
-#time_window: 6e-9
-
-#python:
-from user_libs.antennas.GSSI import antenna_like_GSSI_1500
-antenna_like_GSSI_1500(0.125, 0.094, 0.100, resolution=0.001)
-#end_python:
--- a/user_models/antenna_like_GSSI_1500_fs.py
+++ b/user_models/antenna_like_GSSI_1500_fs.py
@@ -0,0 +1,35 @@
+from pathlib import Path
+
+import gprMax
+from user_libs.antennas.GSSI import antenna_like_GSSI_1500
+
+# File path for output
+fn = Path(__file__)
+
+# Discretisation
+dl = 0.002
+
+# Domain
+x = 0.250
+y = 0.188
+z = 0.220
+
+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)
+
+scene.add(title)
+scene.add(domain)
+scene.add(dxdydz)
+scene.add(time_window)
+
+# Import antenna model and add to model
+gssi_objects = antenna_like_GSSI_1500(0.125, 0.094, 0.100, resolution=dl)
+for obj in gssi_objects:
+    scene.add(obj)
+
+# Run model
+gprMax.run(scenes=[scene], geometry_only=False, outputfile=fn)
--- a/user_models/antenna_like_GSSI_1500_patterns_E.in
+++ b/user_models/antenna_like_GSSI_1500_patterns_E.in
@@ -1,48 +0,0 @@
-#title: GSSI 1.5GHz antenna field patterns
-#dx_dy_dz: 0.001 0.001 0.001
-#pml_cells: 14
-
-#python:
-import os
-import numpy as np
-
-from gprMax.input_cmd_funcs import *
-from user_libs.antennas.GSSI import antenna_like_GSSI_1500
-
-filename = os.path.splitext(os.path.split(inputfile)[1])[0]
-
-timewindows = np.array([4.5e-9]) # For 0.3m max
-radii = np.linspace(0.1, 0.3, 20)
-theta = np.linspace(3, 357, 60)
-materials = ['5 0 1 0 er5'] # Can add more to list and use selector integer to choose
-selector = 0
-
-fs = np.array([0.040, 0.040, 0.040])
-domain = np.array([2 * fs[0] + 0.170, 2 * fs[1] + 2 * radii[-1], 2 * fs[2] + 2 * radii[-1]])
-antennaposition = np.array([domain[0] / 2, fs[1] + radii[-1], fs[2] + radii[-1]])
-antenna_like_GSSI_1500(antennaposition[0], antennaposition[1], antennaposition[2])
-
-print('#domain: {:.3f} {:.3f} {:.3f}'.format(domain[0], domain[1], domain[2]))
-print('#time_window: {:.3e}'.format(timewindows[selector]))
-
-## Can introduce soil model
-#print('#soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 mySoil')
-#print('#fractal_box: 0 0 0 {} {} {} 1.5 1 1 1 50 mySoil mySoilBox 1'.format(domain[0], domain[1], fs[2] + radii[-1]))
-
-print('#material: {}'.format(materials[selector]))
-print('#box: 0 0 0 {} {} {} {} n'.format(domain[0], domain[1], fs[2] + radii[-1], materials[selector].split()[-1]))
-
-## Save the position of the antenna to file for use when processing results
-np.savetxt(os.path.join(os.path.dirname(inputfile), filename + '_rxsorigin.txt'), antennaposition, fmt="%f")
-
-## Generate receiver points for pattern
-for radius in range(len(radii)):
-    ## E-plane circle (yz plane, x=0, phi=pi/2,3pi/2)
-    x = radii[radius] * np.sin(theta * np.pi /180) * np.cos(90 * np.pi / 180)
-    y = radii[radius] * np.sin(theta * np.pi /180) * np.sin(90 * np.pi / 180)
-    z = radii[radius] * np.cos(theta * np.pi /180)
-    for rxpt in range(len(theta)):
-        print('#rx: {:.3f} {:.3f} {:.3f}'.format(x[rxpt] + antennaposition[0], y[rxpt] + antennaposition[1], z[rxpt] + antennaposition[2]))
-
-geometry_view(0, 0, 0, domain[0], domain[1], domain[2], 0.001, 0.001, 0.001, filename, 'n')
-#end_python:
--- a/user_models/antenna_like_GSSI_1500_patterns_H.in
+++ b/user_models/antenna_like_GSSI_1500_patterns_H.in
@@ -1,48 +0,0 @@
-#title: GSSI 1.5GHz antenna field patterns
-#dx_dy_dz: 0.001 0.001 0.001
-#pml_cells: 14
-
-#python:
-import os
-import numpy as np
-
-from gprMax.input_cmd_funcs import *
-from user_libs.antennas.GSSI import antenna_like_GSSI_1500
-
-filename = os.path.splitext(os.path.split(inputfile)[1])[0]
-
-timewindows = np.array([4.5e-9]) # For 0.3m max
-radii = np.linspace(0.1, 0.3, 20)
-theta = np.linspace(3, 357, 60)
-materials = ['5 0 1 0 er5'] # Can add more to list and use selector integer to choose
-selector = 0
-
-fs = np.array([0.040, 0.040, 0.040])
-domain = np.array([2 * fs[0] + 2 * radii[-1], 2 * fs[1] + 0.107, 2 * fs[2] + 2 * radii[-1]])
-antennaposition = np.array([fs[0] + radii[-1], domain[1] / 2, fs[2] + radii[-1]])
-antenna_like_GSSI_1500(antennaposition[0], antennaposition[1], antennaposition[2])
-
-print('#domain: {:.3f} {:.3f} {:.3f}'.format(domain[0], domain[1], domain[2]))
-print('#time_window: {:.3e}'.format(timewindows[selector]))
-
-## Can introduce soil model
-#print('#soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 mySoil')
-#print('#fractal_box: 0 0 0 {} {} {} 1.5 1 1 1 50 mySoil mySoilBox 1'.format(domain[0], domain[1], fs[2] + radii[-1]))
-
-print('#material: {}'.format(materials[selector]))
-print('#box: 0 0 0 {} {} {} {} n'.format(domain[0], domain[1], fs[2] + radii[-1], materials[selector].split()[-1]))
-
-## Save the position of the antenna to file for use when processing results
-np.savetxt(os.path.join(os.path.dirname(inputfile), filename + '_rxsorigin.txt'), antennaposition, fmt="%f")
-
-## Generate receiver points for pattern
-for radius in range(len(radii)):
-    ## H-plane circle (xz plane, y=0, phi=0,pi)
-    x = radii[radius] * np.sin(theta * np.pi /180) * np.cos(180 * np.pi / 180)
-    y = radii[radius] * np.sin(theta * np.pi /180) * np.sin(180 * np.pi / 180)
-    z = radii[radius] * np.cos(theta * np.pi /180)
-    for rxpt in range(len(theta)):
-        print('#rx: {:.3f} {:.3f} {:.3f}'.format(x[rxpt] + antennaposition[0], y[rxpt] + antennaposition[1], z[rxpt] + antennaposition[2]))
-
-geometry_view(0, 0, 0, domain[0], domain[1], domain[2], 0.001, 0.001, 0.001, filename, 'n')
-#end_python:
--- a/user_models/antenna_like_GSSI_400_fs.in
+++ b/user_models/antenna_like_GSSI_400_fs.in
@@ -1,9 +0,0 @@
-#title: GSSI 400MHz 'like' antenna in free-space
-#domain: 0.380 0.380 0.360
-#dx_dy_dz: 0.001 0.001 0.001
-#time_window: 12e-9
-
-#python:
-from user_libs.antennas.GSSI import antenna_like_GSSI_400
-antenna_like_GSSI_400(0.190, 0.190, 0.140, resolution=0.001)
-#end_python:
--- a/user_models/antenna_like_GSSI_400_fs.py
+++ b/user_models/antenna_like_GSSI_400_fs.py
@@ -0,0 +1,35 @@
+from pathlib import Path
+
+import gprMax
+from user_libs.antennas.GSSI import antenna_like_GSSI_400
+
+# File path for output
+fn = Path(__file__)
+
+# Discretisation
+dl = 0.001
+
+# Domain
+x = 0.380
+y = 0.380
+z = 0.360
+
+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=12e-9)
+
+scene.add(title)
+scene.add(domain)
+scene.add(dxdydz)
+scene.add(time_window)
+
+# Import antenna model and add to model
+gssi_objects = antenna_like_GSSI_400(0.190, 0.190, 0.140, resolution=dl)
+for obj in gssi_objects:
+    scene.add(obj)
+
+# Run model
+gprMax.run(scenes=[scene], geometry_only=False, outputfile=fn)
--- a/user_models/antenna_like_MALA_1200_fs.in
+++ b/user_models/antenna_like_MALA_1200_fs.in
@@ -1,9 +0,0 @@
-#title: MALA 1.2GHz 'like' antenna in free-space
-#domain: 0.264 0.189 0.220
-#dx_dy_dz: 0.001 0.001 0.001
-#time_window: 6e-9
-
-#python:
-from user_libs.antennas.MALA import antenna_like_MALA_1200
-antenna_like_MALA_1200(0.132, 0.095, 0.100, 0.001)
-#end_python:
--- a/user_models/antenna_like_MALA_1200_fs.py
+++ b/user_models/antenna_like_MALA_1200_fs.py
@@ -0,0 +1,35 @@
+from pathlib import Path
+
+import gprMax
+from user_libs.antennas.MALA import antenna_like_MALA_1200
+
+# File path for output
+fn = Path(__file__)
+
+# Discretisation
+dl = 0.002
+
+# Domain
+x = 0.264
+y = 0.189
+z = 0.220
+
+scene = gprMax.Scene()
+
+title = gprMax.Title(name=fn.name)
+domain = gprMax.Domain(p1=(x, y, z))
+dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+time_window = gprMax.TimeWindow(time=6e-9)
+
+scene.add(title)
+scene.add(domain)
+scene.add(dxdydz)
+scene.add(time_window)
+
+# Import antenna model and add to model
+gssi_objects = antenna_like_MALA_1200(0.132, 0.095, 0.100, resolution=dl)
+for obj in gssi_objects:
+    scene.add(obj)
+
+# Run model
+gprMax.run(scenes=[scene], geometry_only=False, outputfile=fn)
--- a/user_models/cylinder_2D_py.in
+++ b/user_models/cylinder_2D_py.in
@@ -1,53 +0,0 @@
-#python:
-from gprMax.input_cmd_funcs import *
-command('title', 'A or B scan from a metal cylinder buried in a dielectric half-space')
-z_dim = 0.002
-resolution = 0.002
-tsim = 3e-9
-B_scan = False
-
-domain = domain(x=64e-2, y=30e-2, z=z_dim)
-dx = dx_dy_dz(resolution, resolution, z_dim)
-time_window(tsim)
-
-material(permittivity=6, conductivity=0,
-         permeability=1, magconductivity=0, name='half_space')
-
-identifier = waveform('ricker', amplitude=1, frequency=1.5e9,
-                      identifier='my_ricker')
-
-if B_scan:
-    x_ant = 8e-2
-else:
-    x_ant = domain.x/2 - 1e-2 # in the middle of the x-axis
-
-tx = hertzian_dipole('z',
-                     x_ant, domain.y - 4e-2, 0, # minus 4 cm in y-direction
-                     identifier)
-rx(tx.x + 2e-2, tx.y, tx.z)     # 2 cm away in x-direction from tx
-
-if B_scan:
-    src_steps(dx=0.8e-2/4)
-    rx_steps(dx=0.8e-2/4)
-
-b0, b1 = box(0, 0, 0,
-             domain.x, domain.y - 4e-2, z_dim, # same as domain, minus 4 cm in y-direction
-             'half_space')
-
-c_x, c_y = (domain.x/2, b1.y - 5e-2) # in the middle of the x-axis and 5 cm down from the half_space
-cylinder(c_x, c_y, 0,
-         c_x, c_y, z_dim,
-         radius=1e-2, material='pec')
-
-# Outputs, geometry and snapshots
-geometry_view(0, 0, 0,
-              domain.x, domain.y, domain.z,
-              dx.x, dx.y, dx.z,
-              'cylinder', 'n')
-
-N = 32
-for i in range(1, N+1):
-    snapshot(0, 0, 0,
-             domain.x, domain.y, domain.z,
-             dx.x, dx.y, dx.z, i*(tsim/N), 'snapshot' + str(i))
-#end_python:
--- a/user_models/cylinder_Bscan_GSSI_1500.in
+++ b/user_models/cylinder_Bscan_GSSI_1500.in
@@ -1,16 +0,0 @@
-#title: B-scan from a metal cylinder buried in a dielectric half-space with a GSSI 1.5GHz 'like' antenna
-#domain: 0.480 0.148 0.235
-#dx_dy_dz: 0.001 0.001 0.001
-#time_window: 6e-9
-
-#material: 6 0 1 0 half_space
-
-#box: 0 0 0 0.480 0.148 0.170 half_space
-#cylinder: 0.240 0 0.080 0.240 0.148 0.080 0.010 pec
-
-#python:
-from user_libs.antennas.GSSI import antenna_like_GSSI_1500
-antenna_like_GSSI_1500(0.105 + current_model_run * 0.005, 0.074, 0.170, 0.001)
-#end_python:
-
-geometry_view: 0 0 0 0.480 0.148 0.235 0.001 0.001 0.001 cylinder_GSSI_1500 n
--- a/user_models/heterogeneous_soil.py
+++ b/user_models/heterogeneous_soil.py
@@ -1,42 +0,0 @@
-import gprMax
-
-# single objects
-dxdydz = gprMax.Discretisation(p1=(1e-3, 1e-3, 1e-3))
-messages = gprMax.Messages(yn='y')
-tw = gprMax.TimeWindow(time=6e-9)
-domain = gprMax.Domain(p1=(0.15, 0.15, 0.1))
-title = gprMax.Title(name='Heterogeneous soil using a stochastic distribution of dielectric properties given by a mixing model from Peplinski')
-waveform1 = gprMax.Waveform(wave_type='ricker', amp=1, freq=1.5e9, id='my_ricker')
-dipole = gprMax.HertzianDipole(p1=(0.045, 0.075, 0.085), polarisation='y', waveform_id='my_ricker')
-
-
-sp = gprMax.SoilPeplinski(sand_fraction=0.5,
-                     clay_fraction=0.5,
-                     bulk_density=2.0,
-                     sand_density=2.66,
-                     water_fraction_lower=0.001,
-                     water_fraction_upper=0.25,
-                     id='my_soil')
-
-fb = gprMax.FractalBox(p1=(0, 0, 0), p2=(0.15, 0.15, 0.070), frac_dim=1.5, weighting=[1, 1, 1], n_materials=50, mixing_model_id='my_soil', id='my_soil_box')
-asf = gprMax.AddSurfaceRoughness(p1=(0, 0, 0.070), p2=(0.15, 0.15, 0.070), frac_dim=1.5, weighting=[1, 1], limits=[0.065, 0.080], fractal_box_id='my_soil_box')
-gv = gprMax.GeometryView(p1=(0, 0, 0), p2=(0.15, 0.15, 0.1), dl=(0.001, 0.001, 0.001), filename='heterogeneous_soil', output_type='n')
-
-
-rx = gprMax.Rx(p1=(0.045, 0.075 + 10e-3, 0.085))
-
-scene = gprMax.Scene()
-scene.add(dxdydz)
-scene.add(messages)
-scene.add(tw)
-scene.add(domain)
-scene.add(title)
-scene.add(waveform1)
-scene.add(dipole)
-scene.add(sp)
-scene.add(fb)
-scene.add(asf)
-scene.add(gv)
-scene.add(rx)
-
-gprMax.run(scenes=[scene], n=1, geometry_only=False, outputfile='mysimulation')