class GeneralError(ValueError):
    """Handles general errors. Subclasses the ValueError class."""

    def __init__(self, message, *args):

        self.message = message
        super(GeneralError, self).__init__(message, *args)
        print(Fore.RED) def store_outputs(iteration, Ex, Ey, Ez, Hx, Hy, Hz, G):
    """Stores field component values for every receiver and transmission line.

    Args:
        iteration (int): Current iteration number.
        Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
        G (class): Grid class instance - holds essential parameters describing the model.
    """

    for rx in G.rxs:
        for output in rx.outputs:
            # Store electric or magnetic field components
            if 'I' not in output:
                field = locals()[output]
                rx.outputs[output][iteration] = field[rx.xcoord, rx.ycoord, rx.zcoord]
            # Store current component
            else:
                func = globals()[output]
                rx.outputs[output][iteration] = func(rx.xcoord, rx.ycoord, rx.zcoord, Hx, Hy, Hz, G)

    for tl in G.transmissionlines:
        tl.Vtotal[iteration] = tl.voltage[tl.antpos]
        tl.Itotal[iteration] = tl.current[tl.antpos] def store_outputs(iteration, Ex, Ey, Ez, Hx, Hy, Hz, G):
def write_hdf5_outputfile(outputfile, G): def write_hdf5_outputfile(outputfile, G):
    """Write an output file in HDF5 format.

    Args:
        outputfile (str): Name of the output file.
        G (class): Grid class instance - holds essential parameters describing the model.
    """

    f = h5py.File(outputfile, 'w')
    f.attrs['gprMax'] = __version__
    f.attrs['Title'] = G.title
    f.attrs['Iterations'] = G.iterations
    f.attrs['nx_ny_nz'] = (G.nx, G.ny, G.nz)
    f.attrs['dx_dy_dz'] = (G.dx, G.dy, G.dz)
    f.attrs['dt'] = G.dt
    nsrc = len(G.voltagesources + G.hertziandipoles + G.magneticdipoles + G.transmissionlines)
    f.attrs['nsrc'] = nsrc
    f.attrs['nrx'] = len(G.rxs)
    f.attrs['srcsteps'] = G.srcsteps
    f.attrs['rxsteps'] = G.rxsteps

    # Create group for sources (except transmission lines); add type and positional data attributes
    srclist = G.voltagesources + G.hertziandipoles + G.magneticdipoles
    for srcindex, src in enumerate(srclist):
        grp = f.create_group('/srcs/src' + str(srcindex + 1))
        grp.attrs['Type'] = type(src).__name__
        grp.attrs['Position'] = (src.xcoord * G.dx, src.ycoord * G.dy, src.zcoord * G.dz)

    # Create group for transmission lines; add positional data, line resistance and
    # line discretisation attributes; write arrays for line voltages and currents
    for tlindex, tl in enumerate(G.transmissionlines):
        grp = f.create_group('/tls/tl' + str(tlindex + 1))
        grp.attrs['Position'] = (tl.xcoord * G.dx, tl.ycoord * G.dy, tl.zcoord * G.dz)
        grp.attrs['Resistance'] = tl.resistance
        grp.attrs['dl'] = tl.dl
        # Save incident voltage and current
        grp['Vinc'] = tl.Vinc
        grp['Iinc'] = tl.Iinc
        # Save total voltage and current
        f['/tls/tl' + str(tlindex + 1) + '/Vtotal'] = tl.Vtotal
        f['/tls/tl' + str(tlindex + 1) + '/Itotal'] = tl.Itotal

    # Create group, add positional data and write field component arrays for receivers
    for rxindex, rx in enumerate(G.rxs):
        grp = f.create_group('/rxs/rx' + str(rxindex + 1))
        if rx.ID:
            grp.attrs['Name'] = rx.ID
        grp.attrs['Position'] = (rx.xcoord * G.dx, rx.ycoord * G.dy, rx.zcoord * G.dz)

        for output in rx.outputs:
            f['/rxs/rx' + str(rxindex + 1) + '/' + output] = rx.outputs[output] class FDTDGrid(Grid):
    """
    Holds attributes associated with the entire grid. A convenient
    way for accessing regularly used parameters.
    """

    def __init__(self):
        self.inputfilename = ''
        self.inputdirectory = ''
        self.outputdirectory = ''
        self.title = ''
        self.messages = True
        self.progressbars = self.messages
        self.memoryusage = 0

        # Get information about host machine
        self.hostinfo = None

        # CPU - OpenMP threads
        self.nthreads = 0

        # GPU
        # Threads per block - electric and magnetic field updates
        self.tpb = (256, 1, 1)

        # GPU object
        self.gpu = None

        # Copy snapshot data from GPU to CPU during simulation
        # N.B. This will happen if the requested snapshots are too large to fit
        # on the memory of the GPU. If True this will slow performance significantly
        self.snapsgpu2cpu = False

        # Threshold (dB) down from maximum power (0dB) of main frequency used
        # to calculate highest frequency for numerical dispersion analysis
        self.highestfreqthres = 40
        # Maximum allowable percentage physical phase-velocity phase error
        self.maxnumericaldisp = 2
        # Minimum grid sampling of smallest wavelength for physical wave propagation
        self.mingridsampling = 3

        self.nx = 0
        self.ny = 0
        self.nz = 0
        self.dx = 0
        self.dy = 0
        self.dz = 0
        self.dt = 0
        self.mode = None
        self.iterations = 0
        self.timewindow = 0

        # 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.
        self.pmlthickness = OrderedDict((key, 10) for key in PML.boundaryIDs)
        self.cfs = []
        self.pmls = []
        self.pmlformulation = 'HORIPML'

        self.materials = []
        self.mixingmodels = []
        self.averagevolumeobjects = True
        self.fractalvolumes = []
        self.geometryviews = []
        self.geometryobjectswrite = []
        self.waveforms = []
        self.voltagesources = []
        self.hertziandipoles = []
        self.magneticdipoles = []
        self.transmissionlines = []
        self.rxs = []
        self.srcsteps = [0, 0, 0]
        self.rxsteps = [0, 0, 0]
        self.snapshots = []

    def initialise_geometry_arrays(self):
        """
        Initialise an array for volumetric material IDs (solid);
            boolean arrays for specifying whether materials can have dielectric smoothing (rigid);
            and an array for cell edge IDs (ID).
        Solid and ID arrays are initialised to free_space (one);
            rigid arrays to allow dielectric smoothing (zero).
        """
        self.solid = np.ones((self.nx, self.ny, self.nz), dtype=np.uint32)
        self.rigidE = np.zeros((12, self.nx, self.ny, self.nz), dtype=np.int8)
        self.rigidH = np.zeros((6, self.nx, self.ny, self.nz), dtype=np.int8)
        self.ID = np.ones((6, self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.uint32)
        self.IDlookup = {'Ex': 0, 'Ey': 1, 'Ez': 2, 'Hx': 3, 'Hy': 4, 'Hz': 5}

    def initialise_field_arrays(self):
        """Initialise arrays for the electric and magnetic field components."""
        self.Ex = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype)
        self.Ey = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype)
        self.Ez = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype)
        self.Hx = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype)
        self.Hy = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype)
        self.Hz = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype)

    def initialise_std_update_coeff_arrays(self):
        """Initialise arrays for storing update coefficients."""
        self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=floattype)
        self.updatecoeffsH = np.zeros((len(self.materials), 5), dtype=floattype)

    def initialise_dispersive_arrays(self):
        """Initialise arrays for storing coefficients when there are dispersive materials present."""
        self.Tx = np.zeros((Material.maxpoles, self.nx + 1, self.ny + 1, self.nz + 1), dtype=complextype)
        self.Ty = np.zeros((Material.maxpoles, self.nx + 1, self.ny + 1, self.nz + 1), dtype=complextype)
        self.Tz = np.zeros((Material.maxpoles, self.nx + 1, self.ny + 1, self.nz + 1), dtype=complextype)
        self.updatecoeffsdispersive = np.zeros((len(self.materials), 3 * Material.maxpoles), dtype=complextype)

    def memory_estimate_basic(self):
        """Estimate the amount of memory (RAM) required to run a model."""

        stdoverhead = 50e6

        solidarray = self.nx * self.ny * self.nz * np.dtype(np.uint32).itemsize

        # 12 x rigidE array components + 6 x rigidH array components
        rigidarrays = (12 + 6) * self.nx * self.ny * self.nz * np.dtype(np.int8).itemsize

        # 6 x field arrays + 6 x ID arrays
        fieldarrays = (6 + 6) * (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(floattype).itemsize

        # PML arrays
        pmlarrays = 0
        for (k, v) in self.pmlthickness.items():
            if v > 0:
                if 'x' in k:
                    pmlarrays += ((v + 1) * self.ny * (self.nz + 1))
                    pmlarrays += ((v + 1) * (self.ny + 1) * self.nz)
                    pmlarrays += (v * self.ny * (self.nz + 1))
                    pmlarrays += (v * (self.ny + 1) * self.nz)
                elif 'y' in k:
                    pmlarrays += (self.nx * (v + 1) * (self.nz + 1))
                    pmlarrays += ((self.nx + 1) * (v + 1) * self.nz)
                    pmlarrays += ((self.nx + 1) * v * self.nz)
                    pmlarrays += (self.nx * v * (self.nz + 1))
                elif 'z' in k:
                    pmlarrays += (self.nx * (self.ny + 1) * (v + 1))
                    pmlarrays += ((self.nx + 1) * self.ny * (v + 1))
                    pmlarrays += ((self.nx + 1) * self.ny * v)
                    pmlarrays += (self.nx * (self.ny + 1) * v)

        self.memoryusage = int(stdoverhead + fieldarrays + solidarray + rigidarrays + pmlarrays)

    def memory_check(self, snapsmemsize=0):
        """Check if the required amount of memory (RAM) is available on the host and GPU if specified.

        Args:
            snapsmemsize (int): amount of memory (bytes) required to store all requested snapshots
        """

        # Check if model can be built and/or run on host
        if self.memoryusage > self.hostinfo['ram']:
            raise GeneralError('Memory (RAM) required ~{} exceeds {} detected!\n'.format(human_size(self.memoryusage), human_size(self.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)))

        # Check if model can be run on specified GPU if required
        if self.gpu is not None:
            if self.memoryusage - snapsmemsize > self.gpu.totalmem:
                raise GeneralError('Memory (RAM) required ~{} exceeds {} detected on specified {} - {} GPU!\n'.format(human_size(self.memoryusage), human_size(self.gpu.totalmem, a_kilobyte_is_1024_bytes=True), self.gpu.deviceID, self.gpu.name))

            # If the required memory without the snapshots will fit on the GPU then transfer and store snaphots on host
            if snapsmemsize != 0 and self.memoryusage - snapsmemsize < self.gpu.totalmem:
                self.snapsgpu2cpu = True

    def gpu_set_blocks_per_grid(self):
        """Set the blocks per grid size used for updating the electric and magnetic field arrays on a GPU."""
        self.bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / self.tpb[0])), 1, 1)

    def gpu_initialise_arrays(self):
        """Initialise standard field arrays on GPU."""

        import pycuda.gpuarray as gpuarray

        self.ID_gpu = gpuarray.to_gpu(self.ID)
        self.Ex_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
        self.Ey_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
        self.Ez_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
        self.Hx_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
        self.Hy_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
        self.Hz_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))

    def gpu_initialise_dispersive_arrays(self):
        """Initialise dispersive material coefficient arrays on GPU."""

        import pycuda.gpuarray as gpuarray

        self.Tx_gpu = gpuarray.to_gpu(self.Tx)
        self.Ty_gpu = gpuarray.to_gpu(self.Ty)
        self.Tz_gpu = gpuarray.to_gpu(self.Tz)
        self.updatecoeffsdispersive_gpu = gpuarray.to_gpu(self.updatecoeffsdispersive) def dispersion_analysis(G):
    """
    Analysis of numerical dispersion (Taflove et al, 2005, p112) -
        worse case of maximum frequency and minimum wavelength

    Args:
        G (class): Grid class instance - holds essential parameters describing the model.

    Returns:
        results (dict): Results from dispersion analysis
    """

    # Physical phase velocity error (percentage); grid sampling density;
    # material with maximum permittivity; maximum significant frequency; error message
    results = {'deltavp': False, 'N': False, 'material': False, 'maxfreq': [], 'error': ''}

    # Find maximum significant frequency
    if G.waveforms:
        for waveform in G.waveforms:
            if waveform.type == 'sine' or waveform.type == 'contsine':
                results['maxfreq'].append(4 * waveform.freq)

            elif waveform.type == 'impulse':
                results['error'] = 'impulse waveform used.'

            else:
                # User-defined waveform
                if waveform.type == 'user':
                    iterations = G.iterations

                # Built-in waveform
                else:
                    # Time to analyse waveform - 4*pulse_width as using entire
                    # time window can result in demanding FFT
                    waveform.calculate_coefficients()
                    iterations = round_value(4 * waveform.chi / G.dt)
                    if iterations > G.iterations:
                        iterations = G.iterations

                waveformvalues = np.zeros(G.iterations)
                for iteration in range(G.iterations):
                    waveformvalues[iteration] = waveform.calculate_value(iteration * G.dt, G.dt)

                # Ensure source waveform is not being overly truncated before attempting any FFT
                if np.abs(waveformvalues[-1]) < np.abs(np.amax(waveformvalues)) / 100:
                    # FFT
                    freqs, power = fft_power(waveformvalues, G.dt)
                    # Get frequency for max power
                    freqmaxpower = np.where(np.isclose(power, 0))[0][0]

                    # Set maximum frequency to a threshold drop from maximum power, ignoring DC value
                    try:
                        freqthres = np.where(power[freqmaxpower:] < -G.highestfreqthres)[0][0] + freqmaxpower
                        results['maxfreq'].append(freqs[freqthres])
                    except ValueError:
                        results['error'] = 'unable to calculate maximum power from waveform, most likely due to undersampling.'

                # Ignore case where someone is using a waveform with zero amplitude, i.e. on a receiver
                elif waveform.amp == 0:
                    pass

                # If waveform is truncated don't do any further analysis
                else:
                    results['error'] = 'waveform does not fit within specified time window and is therefore being truncated.'
    else:
        results['error'] = 'no waveform detected.'

    if results['maxfreq']:
        results['maxfreq'] = max(results['maxfreq'])

        # Find minimum wavelength (material with maximum permittivity)
        maxer = 0
        matmaxer = ''
        for x in G.materials:
            if x.se != float('inf'):
                er = x.er
                # If there are dispersive materials calculate the complex relative permittivity
                # at maximum frequency and take the real part
                if x.poles > 0:
                    er = x.calculate_er(results['maxfreq'])
                    er = er.real
                if er > maxer:
                    maxer = er
                    matmaxer = x.ID
        results['material'] = next(x for x in G.materials if x.ID == matmaxer)

        # Minimum velocity
        minvelocity = c / np.sqrt(maxer)

        # Minimum wavelength
        minwavelength = minvelocity / results['maxfreq']

        # Maximum spatial step
        if '3D' in G.mode:
            delta = max(G.dx, G.dy, G.dz)
        elif '2D' in G.mode:
            if G.nx == 1:
                delta = max(G.dy, G.dz)
            elif G.ny == 1:
                delta = max(G.dx, G.dz)
            elif G.nz == 1:
                delta = max(G.dx, G.dy)

        # Courant stability factor
        S = (c * G.dt) / delta

        # Grid sampling density
        results['N'] = minwavelength / delta

        # Check grid sampling will result in physical wave propagation
        if int(np.floor(results['N'])) >= G.mingridsampling:
            # Numerical phase velocity
            vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N'])))

            # Physical phase velocity error (percentage)
            results['deltavp'] = (((vp * c) - c) / c) * 100

        # Store rounded down value of grid sampling density
        results['N'] = int(np.floor(results['N']))

    return results def process_geometrycmds(geometry, G):
    """
    This function checks the validity of command parameters, creates instances
    of classes of parameters, and calls functions to directly set arrays
    solid, rigid and ID.

    Args:
        geometry (list): Geometry commands in the model
    """

    # Disable progress bar if on Windows as it does not update properly
    # when messages are printed for geometry
    if sys.platform == 'win32':
        progressbars = False
    else:
        progressbars = not G.progressbars

    for object in tqdm(geometry, desc='Processing geometry related cmds', unit='cmds', ncols=get_terminal_width() - 1, file=sys.stdout, disable=progressbars):
        tmp = object.split()

        if tmp[0] == '#geometry_objects_read:':
            if len(tmp) != 6:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly five parameters')

            xs = round_value(float(tmp[1]) / G.dx)
            ys = round_value(float(tmp[2]) / G.dy)
            zs = round_value(float(tmp[3]) / G.dz)
            geofile = tmp[4]
            matfile = tmp[5]

            # 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))

            matstr = os.path.splitext(os.path.split(matfile)[1])[0]
            numexistmaterials = len(G.materials)

            # Read materials from file
            with open(matfile, 'r') as f:
                # Read any lines that begin with a hash. Strip out any newline
                # characters and comments that must begin with double hashes.
                materials = [line.rstrip() + '{' + matstr + '}\n' for line in f if(line.startswith('#') and not line.startswith('##') and line.rstrip('\n'))]

            # Check validity of command names
            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)

            # Update material type
            for material in G.materials:
                if material.numID >= numexistmaterials:
                    if material.type:
                        material.type += ',\nimported'
                    else:
                        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))

            # Open geometry object file and read/check spatial resolution attribute
            f = h5py.File(geofile, 'r')
            dx_dy_dz = f.attrs['dx_dy_dz']
            if round_value(dx_dy_dz[0] / G.dx) != 1 or round_value(dx_dy_dz[1] / G.dy) != 1 or round_value(dx_dy_dz[2] / G.dz) != 1:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the spatial resolution of the geometry objects file to match the spatial resolution of the model')

            data = f['/data'][:]

            # Should be int16 to allow for -1 which indicates background, i.e.
            # don't build anything, but AustinMan/Woman maybe uint16
            if data.dtype != 'int16':
                data = data.astype('int16')

            # Check that there are no values in the data greater than the maximum index for the specified materials
            if np.amax(data) > len(materials) - 1:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' found data value(s) ({}) in the geometry objects file greater than the maximum index for the specified materials ({})'.format(np.amax(data), len(materials) - 1))

            # Look to see if rigid and ID arrays are present (these should be
            # present if the original geometry objects were written from gprMax)
            try:
                rigidE = f['/rigidE'][:]
                rigidH = f['/rigidH'][:]
                ID = f['/ID'][:]
                G.solid[xs:xs + data.shape[0], ys:ys + data.shape[1], zs:zs + data.shape[2]] = data + numexistmaterials
                G.rigidE[:, xs:xs + rigidE.shape[1], ys:ys + rigidE.shape[2], zs:zs + rigidE.shape[3]] = rigidE
                G.rigidH[:, xs:xs + rigidH.shape[1], ys:ys + rigidH.shape[2], zs:zs + rigidH.shape[3]] = rigidH
                G.ID[:, xs:xs + ID.shape[1], ys:ys + ID.shape[2], zs:zs + ID.shape[3]] = ID + numexistmaterials
                if G.messages:
                    tqdm.write('Geometry objects from file {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile))
            except KeyError:
                averaging = False
                build_voxels_from_array(xs, ys, zs, numexistmaterials, averaging, data, G.solid, G.rigidE, G.rigidH, G.ID)
                if G.messages:
                    tqdm.write('Geometry objects from file (voxels only) {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile))

        elif tmp[0] == '#edge:':
            if len(tmp) != 8:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly seven parameters')

            xs = round_value(float(tmp[1]) / G.dx)
            xf = round_value(float(tmp[4]) / G.dx)
            ys = round_value(float(tmp[2]) / G.dy)
            yf = round_value(float(tmp[5]) / G.dy)
            zs = round_value(float(tmp[3]) / G.dz)
            zf = round_value(float(tmp[6]) / G.dz)

            if xs < 0 or xs > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
            if xf < 0 or xf > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper x-coordinate {:g}m is not within the model domain'.format(xf * G.dx))
            if ys < 0 or ys > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower y-coordinate {:g}m is not within the model domain'.format(ys * G.dy))
            if yf < 0 or yf > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper y-coordinate {:g}m is not within the model domain'.format(yf * G.dy))
            if zs < 0 or zs > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower z-coordinate {:g}m is not within the model domain'.format(zs * G.dz))
            if zf < 0 or zf > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper z-coordinate {:g}m is not within the model domain'.format(zf * G.dz))
            if xs > xf or ys > yf or zs > zf:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')

            material = next((x for x in G.materials if x.ID == tmp[7]), None)

            if not material:
                raise CmdInputError('Material with ID {} does not exist'.format(tmp[7]))

            # Check for valid orientations
            # x-orientated wire
            if xs != xf:
                if ys != yf or zs != zf:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the edge is not specified correctly')
                else:
                    for i in range(xs, xf):
                        build_edge_x(i, ys, zs, material.numID, G.rigidE, G.rigidH, G.ID)

            # y-orientated wire
            elif ys != yf:
                if xs != xf or zs != zf:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the edge is not specified correctly')
                else:
                    for j in range(ys, yf):
                        build_edge_y(xs, j, zs, material.numID, G.rigidE, G.rigidH, G.ID)

            # z-orientated wire
            elif zs != zf:
                if xs != xf or ys != yf:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the edge is not specified correctly')
                else:
                    for k in range(zs, zf):
                        build_edge_z(xs, ys, k, material.numID, G.rigidE, G.rigidH, G.ID)

            if G.messages:
                tqdm.write('Edge from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, tmp[7]))

        elif tmp[0] == '#plate:':
            if len(tmp) < 8:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least seven parameters')

            # Isotropic case
            elif len(tmp) == 8:
                materialsrequested = [tmp[7]]

            # Anisotropic case
            elif len(tmp) == 9:
                materialsrequested = [tmp[7:]]

            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

            xs = round_value(float(tmp[1]) / G.dx)
            xf = round_value(float(tmp[4]) / G.dx)
            ys = round_value(float(tmp[2]) / G.dy)
            yf = round_value(float(tmp[5]) / G.dy)
            zs = round_value(float(tmp[3]) / G.dz)
            zf = round_value(float(tmp[6]) / G.dz)

            if xs < 0 or xs > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
            if xf < 0 or xf > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper x-coordinate {:g}m is not within the model domain'.format(xf * G.dx))
            if ys < 0 or ys > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower y-coordinate {:g}m is not within the model domain'.format(ys * G.dy))
            if yf < 0 or yf > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper y-coordinate {:g}m is not within the model domain'.format(yf * G.dy))
            if zs < 0 or zs > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower z-coordinate {:g}m is not within the model domain'.format(zs * G.dz))
            if zf < 0 or zf > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper z-coordinate {:g}m is not within the model domain'.format(zf * G.dz))
            if xs > xf or ys > yf or zs > zf:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')

            # Check for valid orientations
            if xs == xf:
                if ys == yf or zs == zf:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the plate is not specified correctly')

            elif ys == yf:
                if xs == xf or zs == zf:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the plate is not specified correctly')

            elif zs == zf:
                if xs == xf or ys == yf:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the plate is not specified correctly')

            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the plate is not specified correctly')

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            # yz-plane plate
            if xs == xf:
                # Isotropic case
                if len(materials) == 1:
                    numIDx = numIDy = numIDz = materials[0].numID

                # Uniaxial anisotropic case
                elif len(materials) == 2:
                    numIDy = materials[0].numID
                    numIDz = materials[1].numID

                for j in range(ys, yf):
                    for k in range(zs, zf):
                        build_face_yz(xs, j, k, numIDy, numIDz, G.rigidE, G.rigidH, G.ID)

            # xz-plane plate
            elif ys == yf:
                # Isotropic case
                if len(materials) == 1:
                    numIDx = numIDy = numIDz = materials[0].numID

                # Uniaxial anisotropic case
                elif len(materials) == 2:
                    numIDx = materials[0].numID
                    numIDz = materials[1].numID

                for i in range(xs, xf):
                    for k in range(zs, zf):
                        build_face_xz(i, ys, k, numIDx, numIDz, G.rigidE, G.rigidH, G.ID)

            # xy-plane plate
            elif zs == zf:
                # Isotropic case
                if len(materials) == 1:
                    numIDx = numIDy = numIDz = materials[0].numID

                # Uniaxial anisotropic case
                elif len(materials) == 2:
                    numIDx = materials[0].numID
                    numIDy = materials[1].numID

                for i in range(xs, xf):
                    for j in range(ys, yf):
                        build_face_xy(i, j, zs, numIDx, numIDy, G.rigidE, G.rigidH, G.ID)

            if G.messages:
                tqdm.write('Plate from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, ', '.join(materialsrequested)))

        elif tmp[0] == '#triangle:':
            if len(tmp) < 12:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eleven parameters')

            # Isotropic case with no user specified averaging
            elif len(tmp) == 12:
                materialsrequested = [tmp[11]]
                averagetriangularprism = G.averagevolumeobjects

            # Isotropic case with user specified averaging
            elif len(tmp) == 13:
                materialsrequested = [tmp[11]]
                if tmp[12].lower() == 'y':
                    averagetriangularprism = True
                elif tmp[12].lower() == 'n':
                    averagetriangularprism = False
                else:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n')

            # Uniaxial anisotropic case
            elif len(tmp) == 14:
                materialsrequested = tmp[11:]

            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

            x1 = round_value(float(tmp[1]) / G.dx) * G.dx
            y1 = round_value(float(tmp[2]) / G.dy) * G.dy
            z1 = round_value(float(tmp[3]) / G.dz) * G.dz
            x2 = round_value(float(tmp[4]) / G.dx) * G.dx
            y2 = round_value(float(tmp[5]) / G.dy) * G.dy
            z2 = round_value(float(tmp[6]) / G.dz) * G.dz
            x3 = round_value(float(tmp[7]) / G.dx) * G.dx
            y3 = round_value(float(tmp[8]) / G.dy) * G.dy
            z3 = round_value(float(tmp[9]) / G.dz) * G.dz
            thickness = float(tmp[10])

            if x1 < 0 or x2 < 0 or x3 < 0 or x1 > G.nx or x2 > G.nx or x3 > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the one of the x-coordinates is not within the model domain')
            if y1 < 0 or y2 < 0 or y3 < 0 or y1 > G.ny or y2 > G.ny or y3 > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the one of the y-coordinates is not within the model domain')
            if z1 < 0 or z2 < 0 or z3 < 0 or z1 > G.nz or z2 > G.nz or z3 > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the one of the z-coordinates is not within the model domain')
            if thickness < 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for thickness')

            # Check for valid orientations
            # yz-plane triangle
            if x1 == x2 and x2 == x3:
                normal = 'x'
            # xz-plane triangle
            elif y1 == y2 and y2 == y3:
                normal = 'y'
            # xy-plane triangle
            elif z1 == z2 and z2 == z3:
                normal = 'z'
            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the triangle is not specified correctly')

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            if thickness > 0:
                # Isotropic case
                if len(materials) == 1:
                    averaging = materials[0].averagable and averagetriangularprism
                    numID = numIDx = numIDy = numIDz = materials[0].numID

                # Uniaxial anisotropic case
                elif len(materials) == 3:
                    averaging = False
                    numIDx = materials[0].numID
                    numIDy = materials[1].numID
                    numIDz = materials[2].numID
                    requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
                    averagedmaterial = [x for x in G.materials if x.ID == requiredID]
                    if averagedmaterial:
                        numID = averagedmaterial.numID
                    else:
                        numID = len(G.materials)
                        m = Material(numID, requiredID)
                        m.type = 'dielectric-smoothed'
                        # Create dielectric-smoothed constituents for material
                        m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
                        m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
                        m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
                        m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)

                        # Append the new material object to the materials list
                        G.materials.append(m)
            else:
                averaging = False
                # Isotropic case
                if len(materials) == 1:
                    numID = numIDx = numIDy = numIDz = materials[0].numID

                # Uniaxial anisotropic case
                elif len(materials) == 3:
                    # numID requires a value but it will not be used
                    numID = None
                    numIDx = materials[0].numID
                    numIDy = materials[1].numID
                    numIDz = materials[2].numID

            build_triangle(x1, y1, z1, x2, y2, z2, x3, y3, z3, normal, thickness, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)

            if G.messages:
                if thickness > 0:
                    if averaging:
                        dielectricsmoothing = 'on'
                    else:
                        dielectricsmoothing = 'off'
                    tqdm.write('Triangle with coordinates {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m and thickness {:g}m of material(s) {} created, dielectric smoothing is {}.'.format(x1, y1, z1, x2, y2, z2, x3, y3, z3, thickness, ', '.join(materialsrequested), dielectricsmoothing))
                else:
                    tqdm.write('Triangle with coordinates {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m of material(s) {} created.'.format(x1, y1, z1, x2, y2, z2, x3, y3, z3, ', '.join(materialsrequested)))

        elif tmp[0] == '#box:':
            if len(tmp) < 8:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least seven parameters')

            # Isotropic case with no user specified averaging
            elif len(tmp) == 8:
                materialsrequested = [tmp[7]]
                averagebox = G.averagevolumeobjects

            # Isotropic case with user specified averaging
            elif len(tmp) == 9:
                materialsrequested = [tmp[7]]
                if tmp[8].lower() == 'y':
                    averagebox = True
                elif tmp[8].lower() == 'n':
                    averagebox = False
                else:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n')

            # Uniaxial anisotropic case
            elif len(tmp) == 10:
                materialsrequested = tmp[7:]

            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

            xs = round_value(float(tmp[1]) / G.dx)
            xf = round_value(float(tmp[4]) / G.dx)
            ys = round_value(float(tmp[2]) / G.dy)
            yf = round_value(float(tmp[5]) / G.dy)
            zs = round_value(float(tmp[3]) / G.dz)
            zf = round_value(float(tmp[6]) / G.dz)

            if xs < 0 or xs > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
            if xf < 0 or xf > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper x-coordinate {:g}m is not within the model domain'.format(xf * G.dx))
            if ys < 0 or ys > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower y-coordinate {:g}m is not within the model domain'.format(ys * G.dy))
            if yf < 0 or yf > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper y-coordinate {:g}m is not within the model domain'.format(yf * G.dy))
            if zs < 0 or zs > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower z-coordinate {:g}m is not within the model domain'.format(zs * G.dz))
            if zf < 0 or zf > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper z-coordinate {:g}m is not within the model domain'.format(zf * G.dz))
            if xs >= xf or ys >= yf or zs >= zf:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            # Isotropic case
            if len(materials) == 1:
                averaging = materials[0].averagable and averagebox
                numID = numIDx = numIDy = numIDz = materials[0].numID

            # Uniaxial anisotropic case
            elif len(materials) == 3:
                averaging = False
                numIDx = materials[0].numID
                numIDy = materials[1].numID
                numIDz = materials[2].numID
                requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
                averagedmaterial = [x for x in G.materials if x.ID == requiredID]
                if averagedmaterial:
                    numID = averagedmaterial.numID
                else:
                    numID = len(G.materials)
                    m = Material(numID, requiredID)
                    m.type = 'dielectric-smoothed'
                    # Create dielectric-smoothed constituents for material
                    m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
                    m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
                    m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
                    m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)

                    # Append the new material object to the materials list
                    G.materials.append(m)

            build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)

            if G.messages:
                if averaging:
                    dielectricsmoothing = 'on'
                else:
                    dielectricsmoothing = 'off'
                tqdm.write('Box from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created, dielectric smoothing is {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, ', '.join(materialsrequested), dielectricsmoothing))

        elif tmp[0] == '#cylinder:':
            if len(tmp) < 9:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eight parameters')

            # Isotropic case with no user specified averaging
            elif len(tmp) == 9:
                materialsrequested = [tmp[8]]
                averagecylinder = G.averagevolumeobjects

            # Isotropic case with user specified averaging
            elif len(tmp) == 10:
                materialsrequested = [tmp[8]]
                if tmp[9].lower() == 'y':
                    averagecylinder = True
                elif tmp[9].lower() == 'n':
                    averagecylinder = False
                else:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n')

            # Uniaxial anisotropic case
            elif len(tmp) == 11:
                materialsrequested = tmp[8:]

            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

            x1 = round_value(float(tmp[1]) / G.dx) * G.dx
            y1 = round_value(float(tmp[2]) / G.dy) * G.dy
            z1 = round_value(float(tmp[3]) / G.dz) * G.dz
            x2 = round_value(float(tmp[4]) / G.dx) * G.dx
            y2 = round_value(float(tmp[5]) / G.dy) * G.dy
            z2 = round_value(float(tmp[6]) / G.dz) * G.dz
            r = float(tmp[7])

            if r <= 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the radius {:g} should be a positive value.'.format(r))

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            # Isotropic case
            if len(materials) == 1:
                averaging = materials[0].averagable and averagecylinder
                numID = numIDx = numIDy = numIDz = materials[0].numID

            # Uniaxial anisotropic case
            elif len(materials) == 3:
                averaging = False
                numIDx = materials[0].numID
                numIDy = materials[1].numID
                numIDz = materials[2].numID
                requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
                averagedmaterial = [x for x in G.materials if x.ID == requiredID]
                if averagedmaterial:
                    numID = averagedmaterial.numID
                else:
                    numID = len(G.materials)
                    m = Material(numID, requiredID)
                    m.type = 'dielectric-smoothed'
                    # Create dielectric-smoothed constituents for material
                    m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
                    m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
                    m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
                    m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)

                    # Append the new material object to the materials list
                    G.materials.append(m)

            build_cylinder(x1, y1, z1, x2, y2, z2, r, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)

            if G.messages:
                if averaging:
                    dielectricsmoothing = 'on'
                else:
                    dielectricsmoothing = 'off'
                tqdm.write('Cylinder with face centres {:g}m, {:g}m, {:g}m and {:g}m, {:g}m, {:g}m, with radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(x1, y1, z1, x2, y2, z2, r, ', '.join(materialsrequested), dielectricsmoothing))

        elif tmp[0] == '#cylindrical_sector:':
            if len(tmp) < 10:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least nine parameters')

            # Isotropic case with no user specified averaging
            elif len(tmp) == 10:
                materialsrequested = [tmp[9]]
                averagecylindricalsector = G.averagevolumeobjects

            # Isotropic case with user specified averaging
            elif len(tmp) == 11:
                materialsrequested = [tmp[9]]
                if tmp[10].lower() == 'y':
                    averagecylindricalsector = True
                elif tmp[10].lower() == 'n':
                    averagecylindricalsector = False
                else:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n')

            # Uniaxial anisotropic case
            elif len(tmp) == 12:
                materialsrequested = tmp[9:]

            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

            normal = tmp[1].lower()
            ctr1 = float(tmp[2])
            ctr2 = float(tmp[3])
            extent1 = float(tmp[4])
            extent2 = float(tmp[5])
            thickness = extent2 - extent1
            r = float(tmp[6])
            sectorstartangle = 2 * np.pi * (float(tmp[7]) / 360)
            sectorangle = 2 * np.pi * (float(tmp[8]) / 360)

            if normal != 'x' and normal != 'y' and normal != 'z':
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the normal direction must be either x, y or z.')
            if r <= 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the radius {:g} should be a positive value.'.format(r))
            if sectorstartangle < 0 or sectorangle <= 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the starting angle and sector angle should be a positive values.')
            if sectorstartangle >= 2 * np.pi or sectorangle >= 2 * np.pi:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the starting angle and sector angle must be less than 360 degrees.')

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            if thickness > 0:
                # Isotropic case
                if len(materials) == 1:
                    averaging = materials[0].averagable and averagecylindricalsector
                    numID = numIDx = numIDy = numIDz = materials[0].numID

                # Uniaxial anisotropic case
                elif len(materials) == 3:
                    averaging = False
                    numIDx = materials[0].numID
                    numIDy = materials[1].numID
                    numIDz = materials[2].numID
                    requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
                    averagedmaterial = [x for x in G.materials if x.ID == requiredID]
                    if averagedmaterial:
                        numID = averagedmaterial.numID
                    else:
                        numID = len(G.materials)
                        m = Material(numID, requiredID)
                        m.type = 'dielectric-smoothed'
                        # Create dielectric-smoothed constituents for material
                        m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
                        m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
                        m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
                        m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)

                        # Append the new material object to the materials list
                        G.materials.append(m)
            else:
                averaging = False
                # Isotropic case
                if len(materials) == 1:
                    numID = numIDx = numIDy = numIDz = materials[0].numID

                # Uniaxial anisotropic case
                elif len(materials) == 3:
                    # numID requires a value but it will not be used
                    numID = None
                    numIDx = materials[0].numID
                    numIDy = materials[1].numID
                    numIDz = materials[2].numID

            # yz-plane cylindrical sector
            if normal == 'x':
                ctr1 = round_value(ctr1 / G.dy) * G.dy
                ctr2 = round_value(ctr2 / G.dz) * G.dz
                level = round_value(extent1 / G.dx)

            # xz-plane cylindrical sector
            elif normal == 'y':
                ctr1 = round_value(ctr1 / G.dx) * G.dx
                ctr2 = round_value(ctr2 / G.dz) * G.dz
                level = round_value(extent1 / G.dy)

            # xy-plane cylindrical sector
            elif normal == 'z':
                ctr1 = round_value(ctr1 / G.dx) * G.dx
                ctr2 = round_value(ctr2 / G.dy) * G.dy
                level = round_value(extent1 / G.dz)

            build_cylindrical_sector(ctr1, ctr2, level, sectorstartangle, sectorangle, r, normal, thickness, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)

            if G.messages:
                if thickness > 0:
                    if averaging:
                        dielectricsmoothing = 'on'
                    else:
                        dielectricsmoothing = 'off'
                    tqdm.write('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, thickness {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, thickness, ', '.join(materialsrequested), dielectricsmoothing))
                else:
                    tqdm.write('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, of material(s) {} created.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, ', '.join(materialsrequested)))

        elif tmp[0] == '#sphere:':
            if len(tmp) < 6:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least five parameters')

            # Isotropic case with no user specified averaging
            elif len(tmp) == 6:
                materialsrequested = [tmp[5]]
                averagesphere = G.averagevolumeobjects

            # Isotropic case with user specified averaging
            elif len(tmp) == 7:
                materialsrequested = [tmp[5]]
                if tmp[6].lower() == 'y':
                    averagesphere = True
                elif tmp[6].lower() == 'n':
                    averagesphere = False
                else:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n')

            # Uniaxial anisotropic case
            elif len(tmp) == 8:
                materialsrequested = tmp[5:]

            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

            # Centre of sphere
            xc = round_value(float(tmp[1]) / G.dx)
            yc = round_value(float(tmp[2]) / G.dy)
            zc = round_value(float(tmp[3]) / G.dz)
            r = float(tmp[4])

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            # Isotropic case
            if len(materials) == 1:
                averaging = materials[0].averagable and averagesphere
                numID = numIDx = numIDy = numIDz = materials[0].numID

            # Uniaxial anisotropic case
            elif len(materials) == 3:
                averaging = False
                numIDx = materials[0].numID
                numIDy = materials[1].numID
                numIDz = materials[2].numID
                requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
                averagedmaterial = [x for x in G.materials if x.ID == requiredID]
                if averagedmaterial:
                    numID = averagedmaterial.numID
                else:
                    numID = len(G.materials)
                    m = Material(numID, requiredID)
                    m.type = 'dielectric-smoothed'
                    # Create dielectric-smoothed constituents for material
                    m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
                    m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
                    m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
                    m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)

                    # Append the new material object to the materials list
                    G.materials.append(m)

            build_sphere(xc, yc, zc, r, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)

            if G.messages:
                if averaging:
                    dielectricsmoothing = 'on'
                else:
                    dielectricsmoothing = 'off'
                tqdm.write('Sphere with centre {:g}m, {:g}m, {:g}m, radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(xc * G.dx, yc * G.dy, zc * G.dz, r, ', '.join(materialsrequested), dielectricsmoothing))

        elif tmp[0] == '#fractal_box:':
            # Default is no dielectric smoothing for a fractal box
            averagefractalbox = False

            if len(tmp) < 14:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least thirteen parameters')
            elif len(tmp) == 14:
                seed = None
            elif len(tmp) == 15:
                seed = int(tmp[14])
            elif len(tmp) == 16:
                seed = int(tmp[14])
                if tmp[15].lower() == 'y':
                    averagefractalbox = True
                elif tmp[15].lower() == 'n':
                    averagefractalbox = False
                else:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n')
            else:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

            xs = round_value(float(tmp[1]) / G.dx)
            xf = round_value(float(tmp[4]) / G.dx)
            ys = round_value(float(tmp[2]) / G.dy)
            yf = round_value(float(tmp[5]) / G.dy)
            zs = round_value(float(tmp[3]) / G.dz)
            zf = round_value(float(tmp[6]) / G.dz)

            if xs < 0 or xs > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
            if xf < 0 or xf > G.nx:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper x-coordinate {:g}m is not within the model domain'.format(xf * G.dx))
            if ys < 0 or ys > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower y-coordinate {:g}m is not within the model domain'.format(ys * G.dy))
            if yf < 0 or yf > G.ny:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper y-coordinate {:g}m is not within the model domain'.format(yf * G.dy))
            if zs < 0 or zs > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower z-coordinate {:g}m is not within the model domain'.format(zs * G.dz))
            if zf < 0 or zf > G.nz:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper z-coordinate {:g}m is not within the model domain'.format(zf * G.dz))
            if xs >= xf or ys >= yf or zs >= zf:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')
            if float(tmp[7]) < 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal dimension')
            if float(tmp[8]) < 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal weighting in the x direction')
            if float(tmp[9]) < 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal weighting in the y direction')
            if float(tmp[10]) < 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal weighting in the z direction')
            if round_value(tmp[11]) < 0:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the number of bins')

            # Find materials to use to build fractal volume, either from mixing models or normal materials
            mixingmodel = next((x for x in G.mixingmodels if x.ID == tmp[12]), None)
            material = next((x for x in G.materials if x.ID == tmp[12]), None)
            nbins = round_value(tmp[11])

            if mixingmodel:
                if nbins == 1:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' must be used with more than one material from the mixing model.')
                # Create materials from mixing model as number of bins now known from fractal_box command
                mixingmodel.calculate_debye_properties(nbins, G)
            elif not material:
                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' mixing model or material with ID {} does not exist'.format(tmp[12]))

            volume = FractalVolume(xs, xf, ys, yf, zs, zf, float(tmp[7]))
            volume.ID = tmp[13]
            volume.operatingonID = tmp[12]
            volume.nbins = nbins
            volume.seed = seed
            volume.weighting = np.array([float(tmp[8]), float(tmp[9]), float(tmp[10])])
            volume.averaging = averagefractalbox

            if G.messages:
                if volume.averaging:
                    dielectricsmoothing = 'on'
                else:
                    dielectricsmoothing = 'off'
                tqdm.write('Fractal box {} from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with {}, fractal dimension {:g}, fractal weightings {:g}, {:g}, {:g}, fractal seeding {}, with {} material(s) created, dielectric smoothing is {}.'.format(volume.ID, xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, volume.operatingonID, volume.dimension, volume.weighting[0], volume.weighting[1], volume.weighting[2], volume.seed, volume.nbins, dielectricsmoothing))

            G.fractalvolumes.append(volume)

            # Search and process any modifiers for the fractal box
            for object in geometry:
                tmp = object.split()

                if tmp[0] == '#add_surface_roughness:':
                    if len(tmp) < 13:
                        raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least twelve parameters')
                    elif len(tmp) == 13:
                        seed = None
                    elif len(tmp) == 14:
                        seed = int(tmp[13])
                    else:
                        raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

                    # Only process rough surfaces for this fractal volume
                    if tmp[12] == volume.ID:
                        xs = round_value(float(tmp[1]) / G.dx)
                        xf = round_value(float(tmp[4]) / G.dx)
                        ys = round_value(float(tmp[2]) / G.dy)
                        yf = round_value(float(tmp[5]) / G.dy)
                        zs = round_value(float(tmp[3]) / G.dz)
                        zf = round_value(float(tmp[6]) / G.dz)

                        if xs < 0 or xs > G.nx:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
                        if xf < 0 or xf > G.nx:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper x-coordinate {:g}m is not within the model domain'.format(xf * G.dx))
                        if ys < 0 or ys > G.ny:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower y-coordinate {:g}m is not within the model domain'.format(ys * G.dy))
                        if yf < 0 or yf > G.ny:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper y-coordinate {:g}m is not within the model domain'.format(yf * G.dy))
                        if zs < 0 or zs > G.nz:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower z-coordinate {:g}m is not within the model domain'.format(zs * G.dz))
                        if zf < 0 or zf > G.nz:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper z-coordinate {:g}m is not within the model domain'.format(zf * G.dz))
                        if xs > xf or ys > yf or zs > zf:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')
                        if float(tmp[7]) < 0:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal dimension')
                        if float(tmp[8]) < 0:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal weighting in the first direction of the surface')
                        if float(tmp[9]) < 0:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal weighting in the second direction of the surface')

                        # Check for valid orientations
                        if xs == xf:
                            if ys == yf or zs == zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if xs != volume.xs and xs != volume.xf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
                            fractalrange = (round_value(float(tmp[10]) / G.dx), round_value(float(tmp[11]) / G.dx))
                            # xminus surface
                            if xs == volume.xs:
                                if fractalrange[0] < 0 or fractalrange[1] > volume.xf:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply fractal surface to fractal box as it would exceed either the upper coordinates of the fractal box or the domain in the x direction')
                                requestedsurface = 'xminus'
                            # xplus surface
                            elif xf == volume.xf:
                                if fractalrange[0] < volume.xs or fractalrange[1] > G.nx:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply fractal surface to fractal box as it would exceed either the lower coordinates of the fractal box or the domain in the x direction')
                                requestedsurface = 'xplus'

                        elif ys == yf:
                            if xs == xf or zs == zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if ys != volume.ys and ys != volume.yf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
                            fractalrange = (round_value(float(tmp[10]) / G.dy), round_value(float(tmp[11]) / G.dy))
                            # yminus surface
                            if ys == volume.ys:
                                if fractalrange[0] < 0 or fractalrange[1] > volume.yf:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply fractal surface to fractal box as it would exceed either the upper coordinates of the fractal box or the domain in the y direction')
                                requestedsurface = 'yminus'
                            # yplus surface
                            elif yf == volume.yf:
                                if fractalrange[0] < volume.ys or fractalrange[1] > G.ny:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply fractal surface to fractal box as it would exceed either the lower coordinates of the fractal box or the domain in the y direction')
                                requestedsurface = 'yplus'

                        elif zs == zf:
                            if xs == xf or ys == yf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if zs != volume.zs and zs != volume.zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
                            fractalrange = (round_value(float(tmp[10]) / G.dz), round_value(float(tmp[11]) / G.dz))
                            # zminus surface
                            if zs == volume.zs:
                                if fractalrange[0] < 0 or fractalrange[1] > volume.zf:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply fractal surface to fractal box as it would exceed either the upper coordinates of the fractal box or the domain in the x direction')
                                requestedsurface = 'zminus'
                            # zplus surface
                            elif zf == volume.zf:
                                if fractalrange[0] < volume.zs or fractalrange[1] > G.nz:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply fractal surface to fractal box as it would exceed either the lower coordinates of the fractal box or the domain in the z direction')
                                requestedsurface = 'zplus'

                        else:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')

                        surface = FractalSurface(xs, xf, ys, yf, zs, zf, float(tmp[7]))
                        surface.surfaceID = requestedsurface
                        surface.fractalrange = fractalrange
                        surface.operatingonID = volume.ID
                        surface.seed = seed
                        surface.weighting = np.array([float(tmp[8]), float(tmp[9])])

                        # List of existing surfaces IDs
                        existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
                        if surface.surfaceID in existingsurfaceIDs:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' has already been used on the {} surface'.format(surface.surfaceID))

                        surface.generate_fractal_surface(G)
                        volume.fractalsurfaces.append(surface)

                        if G.messages:
                            tqdm.write('Fractal surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal weightings {:g}, {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, surface.dimension, surface.weighting[0], surface.weighting[1], surface.seed, float(tmp[10]), float(tmp[11]), surface.operatingonID))

                if tmp[0] == '#add_surface_water:':
                    if len(tmp) != 9:
                        raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly eight parameters')

                    # Only process surfaces for this fractal volume
                    if tmp[8] == volume.ID:
                        xs = round_value(float(tmp[1]) / G.dx)
                        xf = round_value(float(tmp[4]) / G.dx)
                        ys = round_value(float(tmp[2]) / G.dy)
                        yf = round_value(float(tmp[5]) / G.dy)
                        zs = round_value(float(tmp[3]) / G.dz)
                        zf = round_value(float(tmp[6]) / G.dz)
                        depth = float(tmp[7])

                        if xs < 0 or xs > G.nx:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
                        if xf < 0 or xf > G.nx:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper x-coordinate {:g}m is not within the model domain'.format(xf * G.dx))
                        if ys < 0 or ys > G.ny:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower y-coordinate {:g}m is not within the model domain'.format(ys * G.dy))
                        if yf < 0 or yf > G.ny:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper y-coordinate {:g}m is not within the model domain'.format(yf * G.dy))
                        if zs < 0 or zs > G.nz:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower z-coordinate {:g}m is not within the model domain'.format(zs * G.dz))
                        if zf < 0 or zf > G.nz:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper z-coordinate {:g}m is not within the model domain'.format(zf * G.dz))
                        if xs > xf or ys > yf or zs > zf:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')
                        if depth <= 0:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the depth of water')

                        # Check for valid orientations
                        if xs == xf:
                            if ys == yf or zs == zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if xs != volume.xs and xs != volume.xf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
                            # xminus surface
                            if xs == volume.xs:
                                requestedsurface = 'xminus'
                            # xplus surface
                            elif xf == volume.xf:
                                requestedsurface = 'xplus'
                            filldepthcells = round_value(depth / G.dx)
                            filldepth = filldepthcells * G.dx

                        elif ys == yf:
                            if xs == xf or zs == zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if ys != volume.ys and ys != volume.yf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
                            # yminus surface
                            if ys == volume.ys:
                                requestedsurface = 'yminus'
                            # yplus surface
                            elif yf == volume.yf:
                                requestedsurface = 'yplus'
                            filldepthcells = round_value(depth / G.dy)
                            filldepth = filldepthcells * G.dy

                        elif zs == zf:
                            if xs == xf or ys == yf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if zs != volume.zs and zs != volume.zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
                            # zminus surface
                            if zs == volume.zs:
                                requestedsurface = 'zminus'
                            # zplus surface
                            elif zf == volume.zf:
                                requestedsurface = 'zplus'
                            filldepthcells = round_value(depth / G.dz)
                            filldepth = filldepthcells * G.dz

                        else:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')

                        surface = next((x for x in volume.fractalsurfaces if x.surfaceID == requestedsurface), None)
                        if not surface:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' specified surface {} does not have a rough surface applied'.format(requestedsurface))

                        surface.filldepth = filldepthcells

                        # Check that requested fill depth falls within range of surface roughness
                        if surface.filldepth < surface.fractalrange[0] or surface.filldepth > surface.fractalrange[1]:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a value for the depth of water that lies with the range of the requested surface roughness')

                        # Check to see if water has been already defined as a material
                        if not any(x.ID == 'water' for x in G.materials):
                            m = Material(len(G.materials), 'water')
                            m.averagable = False
                            m.type = 'builtin, debye'
                            m.er = Material.watereri
                            m.deltaer.append(Material.waterdeltaer)
                            m.tau.append(Material.watertau)
                            G.materials.append(m)
                            if Material.maxpoles == 0:
                                Material.maxpoles = 1

                        # Check if time step for model is suitable for using water
                        water = next((x for x in G.materials if x.ID == 'water'))
                        testwater = next((x for x in water.tau if x < G.dt), None)
                        if testwater:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the time step for the model to be less than the relaxation time required to model water.')

                        if G.messages:
                            tqdm.write('Water on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with depth {:g}m, added to {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, filldepth, surface.operatingonID))

                if tmp[0] == '#add_grass:':
                    if len(tmp) < 12:
                        raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eleven parameters')
                    elif len(tmp) == 12:
                        seed = None
                    elif len(tmp) == 13:
                        seed = int(tmp[12])
                    else:
                        raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')

                    # Only process grass for this fractal volume
                    if tmp[11] == volume.ID:
                        xs = round_value(float(tmp[1]) / G.dx)
                        xf = round_value(float(tmp[4]) / G.dx)
                        ys = round_value(float(tmp[2]) / G.dy)
                        yf = round_value(float(tmp[5]) / G.dy)
                        zs = round_value(float(tmp[3]) / G.dz)
                        zf = round_value(float(tmp[6]) / G.dz)
                        numblades = int(tmp[10])

                        if xs < 0 or xs > G.nx:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
                        if xf < 0 or xf > G.nx:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper x-coordinate {:g}m is not within the model domain'.format(xf * G.dx))
                        if ys < 0 or ys > G.ny:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower y-coordinate {:g}m is not within the model domain'.format(ys * G.dy))
                        if yf < 0 or yf > G.ny:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper y-coordinate {:g}m is not within the model domain'.format(yf * G.dy))
                        if zs < 0 or zs > G.nz:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower z-coordinate {:g}m is not within the model domain'.format(zs * G.dz))
                        if zf < 0 or zf > G.nz:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the upper z-coordinate {:g}m is not within the model domain'.format(zf * G.dz))
                        if xs > xf or ys > yf or zs > zf:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')
                        if float(tmp[7]) < 0:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal dimension')
                        if float(tmp[8]) < 0 or float(tmp[9]) < 0:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the minimum and maximum heights for grass blades')

                        # Check for valid orientations
                        if xs == xf:
                            if ys == yf or zs == zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if xs != volume.xs and xs != volume.xf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' must specify external surfaces on a fractal box')
                            fractalrange = (round_value(float(tmp[8]) / G.dx), round_value(float(tmp[9]) / G.dx))
                            # xminus surface
                            if xs == volume.xs:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' grass can only be specified on surfaces in the positive axis direction')
                            # xplus surface
                            elif xf == volume.xf:
                                if fractalrange[1] > G.nx:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply grass to fractal box as it would exceed the domain size in the x direction')
                                requestedsurface = 'xplus'

                        elif ys == yf:
                            if xs == xf or zs == zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if ys != volume.ys and ys != volume.yf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' must specify external surfaces on a fractal box')
                            fractalrange = (round_value(float(tmp[8]) / G.dy), round_value(float(tmp[9]) / G.dy))
                            # yminus surface
                            if ys == volume.ys:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' grass can only be specified on surfaces in the positive axis direction')
                            # yplus surface
                            elif yf == volume.yf:
                                if fractalrange[1] > G.ny:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply grass to fractal box as it would exceed the domain size in the y direction')
                                requestedsurface = 'yplus'

                        elif zs == zf:
                            if xs == xf or ys == yf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
                            if zs != volume.zs and zs != volume.zf:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' must specify external surfaces on a fractal box')
                            fractalrange = (round_value(float(tmp[8]) / G.dz), round_value(float(tmp[9]) / G.dz))
                            # zminus surface
                            if zs == volume.zs:
                                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' grass can only be specified on surfaces in the positive axis direction')
                            # zplus surface
                            elif zf == volume.zf:
                                if fractalrange[1] > G.nz:
                                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply grass to fractal box as it would exceed the domain size in the z direction')
                                requestedsurface = 'zplus'

                        else:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')

                        surface = FractalSurface(xs, xf, ys, yf, zs, zf, float(tmp[7]))
                        surface.ID = 'grass'
                        surface.surfaceID = requestedsurface
                        surface.seed = seed

                        # Set the fractal range to scale the fractal distribution between zero and one
                        surface.fractalrange = (0, 1)
                        surface.operatingonID = volume.ID
                        surface.generate_fractal_surface(G)
                        if numblades > surface.fractalsurface.shape[0] * surface.fractalsurface.shape[1]:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the specified surface is not large enough for the number of grass blades/roots specified')

                        # Scale the distribution so that the summation is equal to one, i.e. a probability distribution
                        surface.fractalsurface = surface.fractalsurface / np.sum(surface.fractalsurface)

                        # Set location of grass blades using probability distribution
                        # Create 1D vector of probability values from the 2D surface
                        probability1D = np.cumsum(np.ravel(surface.fractalsurface))

                        # Create random numbers between zero and one for the number of blades of grass
                        R = np.random.RandomState(surface.seed)
                        A = R.random_sample(numblades)

                        # Locate the random numbers in the bins created by the 1D vector of probability values, and convert the 1D index back into a x, y index for the original surface.
                        bladesindex = np.unravel_index(np.digitize(A, probability1D), (surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))

                        # Set the fractal range to minimum and maximum heights of the grass blades
                        surface.fractalrange = fractalrange

                        # Set the fractal surface using the pre-calculated spatial distribution and a random height
                        surface.fractalsurface = np.zeros((surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))
                        for i in range(len(bladesindex[0])):
                            surface.fractalsurface[bladesindex[0][i], bladesindex[1][i]] = R.randint(surface.fractalrange[0], surface.fractalrange[1], size=1)

                        # Create grass geometry parameters
                        g = Grass(numblades)
                        g.seed = surface.seed
                        surface.grass.append(g)

                        # Check to see if grass has been already defined as a material
                        if not any(x.ID == 'grass' for x in G.materials):
                            m = Material(len(G.materials), 'grass')
                            m.averagable = False
                            m.type = 'builtin, debye'
                            m.er = Material.grasseri
                            m.deltaer.append(Material.grassdeltaer)
                            m.tau.append(Material.grasstau)
                            G.materials.append(m)
                            if Material.maxpoles == 0:
                                Material.maxpoles = 1

                        # Check if time step for model is suitable for using grass
                        grass = next((x for x in G.materials if x.ID == 'grass'))
                        testgrass = next((x for x in grass.tau if x < G.dt), None)
                        if testgrass:
                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the time step for the model to be less than the relaxation time required to model grass.')

                        volume.fractalsurfaces.append(surface)

                        if G.messages:
                            tqdm.write('{} blades of grass on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(numblades, xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, surface.dimension, surface.seed, float(tmp[8]), float(tmp[9]), surface.operatingonID))

            # Process any modifications to the original fractal box then generate it
            if volume.fractalsurfaces:
                volume.originalxs = volume.xs
                volume.originalxf = volume.xf
                volume.originalys = volume.ys
                volume.originalyf = volume.yf
                volume.originalzs = volume.zs
                volume.originalzf = volume.zf

                # Extend the volume to accomodate any rough surfaces, grass, or roots
                for surface in volume.fractalsurfaces:
                    if surface.surfaceID == 'xminus':
                        if surface.fractalrange[0] < volume.xs:
                            volume.nx += volume.xs - surface.fractalrange[0]
                            volume.xs = surface.fractalrange[0]
                    elif surface.surfaceID == 'xplus':
                        if surface.fractalrange[1] > volume.xf:
                            volume.nx += surface.fractalrange[1] - volume.xf
                            volume.xf = surface.fractalrange[1]
                    elif surface.surfaceID == 'yminus':
                        if surface.fractalrange[0] < volume.ys:
                            volume.ny += volume.ys - surface.fractalrange[0]
                            volume.ys = surface.fractalrange[0]
                    elif surface.surfaceID == 'yplus':
                        if surface.fractalrange[1] > volume.yf:
                            volume.ny += surface.fractalrange[1] - volume.yf
                            volume.yf = surface.fractalrange[1]
                    elif surface.surfaceID == 'zminus':
                        if surface.fractalrange[0] < volume.zs:
                            volume.nz += volume.zs - surface.fractalrange[0]
                            volume.zs = surface.fractalrange[0]
                    elif surface.surfaceID == 'zplus':
                        if surface.fractalrange[1] > volume.zf:
                            volume.nz += surface.fractalrange[1] - volume.zf
                            volume.zf = surface.fractalrange[1]

                # If there is only 1 bin then a normal material is being used, otherwise a mixing model
                if volume.nbins == 1:
                    volume.fractalvolume = np.ones((volume.nx, volume.ny, volume.nz), dtype=floattype)
                    materialnumID = next(x.numID for x in G.materials if x.ID == volume.operatingonID)
                    volume.fractalvolume *= materialnumID
                else:
                    volume.generate_fractal_volume(G)
                    volume.fractalvolume += mixingmodel.startmaterialnum

                volume.generate_volume_mask()

                # Apply any rough surfaces and add any surface water to the 3D mask array
                for surface in volume.fractalsurfaces:
                    if surface.surfaceID == 'xminus':
                        for i in range(surface.fractalrange[0], surface.fractalrange[1]):
                            for j in range(surface.ys, surface.yf):
                                for k in range(surface.zs, surface.zf):
                                    if i > surface.fractalsurface[j - surface.ys, k - surface.zs]:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
                                    elif surface.filldepth > 0 and i > surface.filldepth:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
                                    else:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0

                    elif surface.surfaceID == 'xplus':
                        if not surface.ID:
                            for i in range(surface.fractalrange[0], surface.fractalrange[1]):
                                for j in range(surface.ys, surface.yf):
                                    for k in range(surface.zs, surface.zf):
                                        if i < surface.fractalsurface[j - surface.ys, k - surface.zs]:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
                                        elif surface.filldepth > 0 and i < surface.filldepth:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
                                        else:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
                        elif surface.ID == 'grass':
                            g = surface.grass[0]
                            # Build the blades of the grass
                            blade = 0
                            for j in range(surface.ys, surface.yf):
                                for k in range(surface.zs, surface.zf):
                                    if surface.fractalsurface[j - surface.ys, k - surface.zs] > 0:
                                        height = 0
                                        for i in range(volume.xs, surface.fractalrange[1]):
                                            if i < surface.fractalsurface[j - surface.ys, k - surface.zs] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
                                                y, z = g.calculate_blade_geometry(blade, height)
                                                # Add y, z coordinates to existing location
                                                yy = int(j - volume.ys + y)
                                                zz = int(k - volume.zs + z)
                                                # If these coordinates are outwith fractal volume stop building the blade, otherwise set the mask for grass
                                                if yy < 0 or yy >= volume.mask.shape[1] or zz < 0 or zz >= volume.mask.shape[2]:
                                                    break
                                                else:
                                                    volume.mask[i - volume.xs, yy, zz] = 3
                                                    height += 1
                                        blade += 1

                            # Build the roots of the grass
                            root = 0
                            for j in range(surface.ys, surface.yf):
                                for k in range(surface.zs, surface.zf):
                                    if surface.fractalsurface[j - surface.ys, k - surface.zs] > 0:
                                        depth = 0
                                        i = volume.xf - 1
                                        while i > volume.xs:
                                            if i > volume.originalxf - (surface.fractalsurface[j - surface.ys, k - surface.zs] - volume.originalxf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
                                                y, z = g.calculate_root_geometry(root, depth)
                                                # Add y, z coordinates to existing location
                                                yy = int(j - volume.ys + y)
                                                zz = int(k - volume.zs + z)
                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
                                                if yy < 0 or yy >= volume.mask.shape[1] or zz < 0 or zz >= volume.mask.shape[2]:
                                                    break
                                                else:
                                                    volume.mask[i - volume.xs, yy, zz] = 3
                                                    depth += 1
                                            i -= 1
                                        root += 1

                    elif surface.surfaceID == 'yminus':
                        for i in range(surface.xs, surface.xf):
                            for j in range(surface.fractalrange[0], surface.fractalrange[1]):
                                for k in range(surface.zs, surface.zf):
                                    if j > surface.fractalsurface[i - surface.xs, k - surface.zs]:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
                                    elif surface.filldepth > 0 and j > surface.filldepth:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
                                    else:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0

                    elif surface.surfaceID == 'yplus':
                        if not surface.ID:
                            for i in range(surface.xs, surface.xf):
                                for j in range(surface.fractalrange[0], surface.fractalrange[1]):
                                    for k in range(surface.zs, surface.zf):
                                        if j < surface.fractalsurface[i - surface.xs, k - surface.zs]:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
                                        elif surface.filldepth > 0 and j < surface.filldepth:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
                                        else:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
                        elif surface.ID == 'grass':
                            g = surface.grass[0]
                            # Build the blades of the grass
                            blade = 0
                            for i in range(surface.xs, surface.xf):
                                for k in range(surface.zs, surface.zf):
                                    if surface.fractalsurface[i - surface.xs, k - surface.zs] > 0:
                                        height = 0
                                        for j in range(volume.ys, surface.fractalrange[1]):
                                            if j < surface.fractalsurface[i - surface.xs, k - surface.zs] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
                                                x, z = g.calculate_blade_geometry(blade, height)
                                                # Add x, z coordinates to existing location
                                                xx = int(i - volume.xs + x)
                                                zz = int(k - volume.zs + z)
                                                # If these coordinates are outwith fractal volume stop building the blade, otherwise set the mask for grass
                                                if xx < 0 or xx >= volume.mask.shape[0] or zz < 0 or zz >= volume.mask.shape[2]:
                                                    break
                                                else:
                                                    volume.mask[xx, j - volume.ys, zz] = 3
                                                    height += 1
                                        blade += 1

                            # Build the roots of the grass
                            root = 0
                            for i in range(surface.xs, surface.xf):
                                for k in range(surface.zs, surface.zf):
                                    if surface.fractalsurface[i - surface.xs, k - surface.zs] > 0:
                                        depth = 0
                                        j = volume.yf - 1
                                        while j > volume.ys:
                                            if j > volume.originalyf - (surface.fractalsurface[i - surface.xs, k - surface.zs] - volume.originalyf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
                                                x, z = g.calculate_root_geometry(root, depth)
                                                # Add x, z coordinates to existing location
                                                xx = int(i - volume.xs + x)
                                                zz = int(k - volume.zs + z)
                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
                                                if xx < 0 or xx >= volume.mask.shape[0] or zz < 0 or zz >= volume.mask.shape[2]:
                                                    break
                                                else:
                                                    volume.mask[xx, j - volume.ys, zz] = 3
                                                    depth += 1
                                            j -= 1
                                        root += 1

                    elif surface.surfaceID == 'zminus':
                        for i in range(surface.xs, surface.xf):
                            for j in range(surface.ys, surface.yf):
                                for k in range(surface.fractalrange[0], surface.fractalrange[1]):
                                    if k > surface.fractalsurface[i - surface.xs, j - surface.ys]:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
                                    elif surface.filldepth > 0 and k > surface.filldepth:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
                                    else:
                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0

                    elif surface.surfaceID == 'zplus':
                        if not surface.ID:
                            for i in range(surface.xs, surface.xf):
                                for j in range(surface.ys, surface.yf):
                                    for k in range(surface.fractalrange[0], surface.fractalrange[1]):
                                        if k < surface.fractalsurface[i - surface.xs, j - surface.ys]:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
                                        elif surface.filldepth > 0 and k < surface.filldepth:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
                                        else:
                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
                        elif surface.ID == 'grass':
                            g = surface.grass[0]
                            # Build the blades of the grass
                            blade = 0
                            for i in range(surface.xs, surface.xf):
                                for j in range(surface.ys, surface.yf):
                                    if surface.fractalsurface[i - surface.xs, j - surface.ys] > 0:
                                        height = 0
                                        for k in range(volume.zs, surface.fractalrange[1]):
                                            if k < surface.fractalsurface[i - surface.xs, j - surface.ys] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
                                                x, y = g.calculate_blade_geometry(blade, height)
                                                # Add x, y coordinates to existing location
                                                xx = int(i - volume.xs + x)
                                                yy = int(j - volume.ys + y)
                                                # If these coordinates are outwith the fractal volume stop building the blade, otherwise set the mask for grass
                                                if xx < 0 or xx >= volume.mask.shape[0] or yy < 0 or yy >= volume.mask.shape[1]:
                                                    break
                                                else:
                                                    volume.mask[xx, yy, k - volume.zs] = 3
                                                    height += 1
                                        blade += 1

                            # Build the roots of the grass
                            root = 0
                            for i in range(surface.xs, surface.xf):
                                for j in range(surface.ys, surface.yf):
                                    if surface.fractalsurface[i - surface.xs, j - surface.ys] > 0:
                                        depth = 0
                                        k = volume.zf - 1
                                        while k > volume.zs:
                                            if k > volume.originalzf - (surface.fractalsurface[i - surface.xs, j - surface.ys] - volume.originalzf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
                                                x, y = g.calculate_root_geometry(root, depth)
                                                # Add x, y coordinates to existing location
                                                xx = int(i - volume.xs + x)
                                                yy = int(j - volume.ys + y)
                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
                                                if xx < 0 or xx >= volume.mask.shape[0] or yy < 0 or yy >= volume.mask.shape[1]:
                                                    break
                                                else:
                                                    volume.mask[xx, yy, k - volume.zs] = 3
                                                    depth += 1
                                            k -= 1
                                        root += 1

                # Build voxels from any true values of the 3D mask array
                waternumID = next((x.numID for x in G.materials if x.ID == 'water'), 0)
                grassnumID = next((x.numID for x in G.materials if x.ID == 'grass'), 0)
                data = volume.fractalvolume.astype('int16', order='C')
                mask = volume.mask.copy(order='C')
                build_voxels_from_array_mask(volume.xs, volume.ys, volume.zs, waternumID, grassnumID, volume.averaging, mask, data, G.solid, G.rigidE, G.rigidH, G.ID)

            else:
                if volume.nbins == 1:
                    raise CmdInputError("'" + ' '.join(tmp) + "'" + ' is being used with a single material and no modifications, therefore please use a #box command instead.')
                else:
                    volume.generate_fractal_volume(G)
                    volume.fractalvolume += mixingmodel.startmaterialnum

                data = volume.fractalvolume.astype('int16', order='C')
                build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0, volume.averaging, data, G.solid, G.rigidE, G.rigidH, G.ID) def process_python_include_code(inputfile, usernamespace):
    """Looks for and processes any Python code found in the input file.
    It will ignore any lines that are comments, i.e. begin with a
    double hash (##), and any blank lines. It will also ignore any
    lines that do not begin with a hash (#) after it has processed
    Python commands. It will also process any include file commands
    and insert the contents of the included file at that location.

    Args:
        inputfile (object): File object for input file.
        usernamespace (dict): Namespace that can be accessed by user
                in any Python code blocks in input file.

    Returns:
        processedlines (list): Input commands after Python processing.
    """

    # Strip out any newline characters and comments that must begin with double hashes
    inputlines = [line.rstrip() for line in inputfile if(not line.startswith('##') and line.rstrip('\n'))]

    # Rewind input file in preparation for any subsequent reading function
    inputfile.seek(0)

    # List to hold final processed commands
    processedlines = []

    x = 0
    while(x < len(inputlines)):

        # Process any Python code
        if(inputlines[x].startswith('#python:')):

            # String to hold Python code to be executed
            pythoncode = ''
            x += 1
            while not inputlines[x].startswith('#end_python:'):
                # Add all code in current code block to string
                pythoncode += inputlines[x] + '\n'
                x += 1
                if x == len(inputlines):
                    raise CmdInputError('Cannot find the end of the Python code block, i.e. missing #end_python: command.')
            # Compile code for faster execution
            pythoncompiledcode = compile(pythoncode, '<string>', 'exec')
            # Redirect stdout to a text stream
            sys.stdout = result = StringIO()
            # Execute code block & make available only usernamespace
            exec(pythoncompiledcode, usernamespace)
            # String containing buffer of executed code
            codeout = result.getvalue().split('\n')
            result.close()

            # Reset stdio
            sys.stdout = sys.__stdout__

            # Separate commands from any other generated output
            hashcmds = []
            pythonout = []
            for line in codeout:
                if line.startswith('#'):
                    hashcmds.append(line + '\n')
                elif line:
                    pythonout.append(line)

            # Add commands to a list
            processedlines.extend(hashcmds)

            # Print any generated output that is not commands
            if pythonout:
                print('Python messages (from stdout/stderr): {}\n'.format(pythonout))

        # Add any other commands to list
        elif(inputlines[x].startswith('#')):
            # Add gprMax command to list
            inputlines[x] += ('\n')
            processedlines.append(inputlines[x])

        x += 1

    # Process any include file commands
    processedlines = process_include_files(processedlines, inputfile)

    return processedlines def write_processed_file(processedlines, appendmodelnumber, G):
    """
    Writes an input file after any Python code and include commands
    in the original input file have been processed.

    Args:
        processedlines (list): Input commands after after processing any
            Python code and include commands.
        appendmodelnumber (str): Text to append to filename.
        G (class): Grid class instance - holds essential parameters describing the model.
    """

    processedfile = os.path.join(G.inputdirectory, os.path.splitext(G.inputfilename)[0] + appendmodelnumber + '_processed.in')

    with open(processedfile, 'w') as f:
        for item in processedlines:
            f.write('{}'.format(item))

    print('Written input commands, after processing any Python code and include commands, to file: {}\n'.format(processedfile)) def check_cmd_names(processedlines, checkessential=True):
    """
    Checks the validity of commands, i.e. are they gprMax commands,
        and that all essential commands are present.

    Args:
        processedlines (list): Input commands after Python processing.
        checkessential (boolean): Perform check to see that all essential commands are present.

    Returns:
        singlecmds (dict): Commands that can only occur once in the model.
        multiplecmds (dict): Commands that can have multiple instances in the model.
        geometry (list): Geometry commands in the model.
    """

    # Dictionaries of available commands
    # Essential commands neccessary to run a gprMax model
    essentialcmds = ['#domain', '#dx_dy_dz', '#time_window']

    # Commands that there should only be one instance of in a model
    singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_formulation', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#taguchi', '#end_taguchi', '#output_dir'], None)

    # Commands that there can be multiple instances of in a model - these will be lists within the dictionary
    multiplecmds = {key: [] for key in ['#geometry_view', '#geometry_objects_write', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#transmission_line', '#rx', '#rx_array', '#snapshot', '#pml_cfs', '#include_file']}

    # Geometry object building commands that there can be multiple instances
    # of in a model - these will be lists within the dictionary
    geometrycmds = ['#geometry_objects_read', '#edge', '#plate', '#triangle', '#box', '#sphere', '#cylinder', '#cylindrical_sector', '#fractal_box', '#add_surface_roughness', '#add_surface_water', '#add_grass']
    # List to store all geometry object commands in order from input file
    geometry = []

    # Check if command names are valid, if essential commands are present, and
    # add command parameters to appropriate dictionary values or lists
    countessentialcmds = 0
    lindex = 0
    while(lindex < len(processedlines)):
        cmd = processedlines[lindex].split(':')
        cmdname = cmd[0]
        cmdparams = cmd[1]

        # Check if there is space between command name and parameters, i.e.
        # check first character of parameter string. Ignore case when there
        # are no parameters for a command, e.g. for #taguchi:
        if ' ' not in cmdparams[0] and len(cmdparams.strip('\n')) != 0:
            raise CmdInputError('There must be a space between the command name and parameters in ' + processedlines[lindex])

        # Check if command name is valid
        if cmdname not in essentialcmds and cmdname not in singlecmds and cmdname not in multiplecmds and cmdname not in geometrycmds:
            raise CmdInputError('Your input file contains an invalid command: ' + cmdname)

        # Count essential commands
        if cmdname in essentialcmds:
            countessentialcmds += 1

        # Assign command parameters as values to dictionary keys
        if cmdname in singlecmds:
            if singlecmds[cmdname] is None:
                singlecmds[cmdname] = cmd[1].strip(' \t\n')
            else:
                raise CmdInputError('You can only have a single instance of ' + cmdname + ' in your model')

        elif cmdname in multiplecmds:
            multiplecmds[cmdname].append(cmd[1].strip(' \t\n'))

        elif cmdname in geometrycmds:
            geometry.append(processedlines[lindex].strip(' \t\n'))

        lindex += 1

    if checkessential:
        if (countessentialcmds < len(essentialcmds)):
            raise CmdInputError('Your input file is missing essential commands required to run a model. Essential commands are: ' + ', '.join(essentialcmds))

    return singlecmds, multiplecmds, geometry def process_multicmds(multicmds, G):
    """
    Checks the validity of command parameters and creates instances of
        classes of parameters.

    Args:
        multicmds (dict): Commands that can have multiple instances in the model.
        G (class): Grid class instance - holds essential parameters describing the model.
    """

    # Check if coordinates are within the bounds of the grid
    def check_coordinates(x, y, z, name=''):
        try:
            G.within_bounds(x=x, y=y, z=z)
        except ValueError as err:
            s = "'{}: {} ' {} {}-coordinate is not within the model domain".format(cmdname, ' '.join(tmp), name, err.args[0])
            raise CmdInputError(s)

    # Waveform definitions
    cmdname = '#waveform'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 4:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly four parameters')
            if tmp[0].lower() not in Waveform.types:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' must have one of the following types {}'.format(','.join(Waveform.types)))
            if float(tmp[2]) <= 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires an excitation frequency value of greater than zero')
            if any(x.ID == tmp[3] for x in G.waveforms):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' with ID {} already exists'.format(tmp[3]))

            w = Waveform()
            w.ID = tmp[3]
            w.type = tmp[0].lower()
            w.amp = float(tmp[1])
            w.freq = float(tmp[2])

            if G.messages:
                print('Waveform {} of type {} with maximum amplitude scaling {:g}, frequency {:g}Hz created.'.format(w.ID, w.type, w.amp, w.freq))

            G.waveforms.append(w)

    # Voltage source
    cmdname = '#voltage_source'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) < 6:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least six parameters')

            # Check polarity & position parameters
            polarisation = tmp[0].lower()
            if polarisation not in ('x', 'y', 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
            if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
            elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
            elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')

            xcoord = G.calculate_coord('x', tmp[1])
            ycoord = G.calculate_coord('y', tmp[2])
            zcoord = G.calculate_coord('z', tmp[3])
            resistance = float(tmp[4])

            check_coordinates(xcoord, ycoord, zcoord)
            if xcoord < G.pmlthickness['x0'] or xcoord > G.nx - G.pmlthickness['xmax'] or ycoord < G.pmlthickness['y0'] or ycoord > G.ny - G.pmlthickness['ymax'] or zcoord < G.pmlthickness['z0'] or zcoord > G.nz - G.pmlthickness['zmax']:
                print(Fore.RED + "WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.' + Style.RESET_ALL)
            if resistance < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a source resistance of zero or greater')

            # Check if there is a waveformID in the waveforms list
            if not any(x.ID == tmp[5] for x in G.waveforms):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' there is no waveform with the identifier {}'.format(tmp[5]))

            v = VoltageSource()
            v.polarisation = polarisation
            v.xcoord = xcoord
            v.ycoord = ycoord
            v.zcoord = zcoord
            v.ID = v.__class__.__name__ + '(' + str(v.xcoord) + ',' + str(v.ycoord) + ',' + str(v.zcoord) + ')'
            v.resistance = resistance
            v.waveformID = tmp[5]

            if len(tmp) > 6:
                # Check source start & source remove time parameters
                start = float(tmp[6])
                stop = float(tmp[7])
                if start < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' delay of the initiation of the source should not be less than zero')
                if stop < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time to remove the source should not be less than zero')
                if stop - start <= 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' duration of the source should not be zero or less')
                v.start = start
                if stop > G.timewindow:
                    v.stop = G.timewindow
                else:
                    v.stop = stop
                startstop = ' start time {:g} secs, finish time {:g} secs '.format(v.start, v.stop)
            else:
                v.start = 0
                v.stop = G.timewindow
                startstop = ' '

            v.calculate_waveform_values(G)

            if G.messages:
                print('Voltage source with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(v.polarisation, v.xcoord * G.dx, v.ycoord * G.dy, v.zcoord * G.dz, v.resistance) + startstop + 'using waveform {} created.'.format(v.waveformID))

            G.voltagesources.append(v)

    # Hertzian dipole
    cmdname = '#hertzian_dipole'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) < 5:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least five parameters')

            # Check polarity & position parameters
            polarisation = tmp[0].lower()
            if polarisation not in ('x', 'y', 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
            if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
            elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
            elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')

            xcoord = G.calculate_coord('x', tmp[1])
            ycoord = G.calculate_coord('y', tmp[2])
            zcoord = G.calculate_coord('z', tmp[3])
            check_coordinates(xcoord, ycoord, zcoord)
            if xcoord < G.pmlthickness['x0'] or xcoord > G.nx - G.pmlthickness['xmax'] or ycoord < G.pmlthickness['y0'] or ycoord > G.ny - G.pmlthickness['ymax'] or zcoord < G.pmlthickness['z0'] or zcoord > G.nz - G.pmlthickness['zmax']:
                print(Fore.RED + "WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.' + Style.RESET_ALL)

            # Check if there is a waveformID in the waveforms list
            if not any(x.ID == tmp[4] for x in G.waveforms):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' there is no waveform with the identifier {}'.format(tmp[4]))

            h = HertzianDipole()
            h.polarisation = polarisation

            # Set length of dipole to grid size in polarisation direction
            if h.polarisation == 'x':
                h.dl = G.dx
            elif h.polarisation == 'y':
                h.dl = G.dy
            elif h.polarisation == 'z':
                h.dl = G.dz

            h.xcoord = xcoord
            h.ycoord = ycoord
            h.zcoord = zcoord
            h.xcoordorigin = xcoord
            h.ycoordorigin = ycoord
            h.zcoordorigin = zcoord
            h.ID = h.__class__.__name__ + '(' + str(h.xcoord) + ',' + str(h.ycoord) + ',' + str(h.zcoord) + ')'
            h.waveformID = tmp[4]

            if len(tmp) > 5:
                # Check source start & source remove time parameters
                start = float(tmp[5])
                stop = float(tmp[6])
                if start < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' delay of the initiation of the source should not be less than zero')
                if stop < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time to remove the source should not be less than zero')
                if stop - start <= 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' duration of the source should not be zero or less')
                h.start = start
                if stop > G.timewindow:
                    h.stop = G.timewindow
                else:
                    h.stop = stop
                startstop = ' start time {:g} secs, finish time {:g} secs '.format(h.start, h.stop)
            else:
                h.start = 0
                h.stop = G.timewindow
                startstop = ' '

            h.calculate_waveform_values(G)

            if G.messages:
                if G.mode == '2D':
                    print('Hertzian dipole is a line source in 2D with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * G.dx, h.ycoord * G.dy, h.zcoord * G.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))
                else:
                    print('Hertzian dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * G.dx, h.ycoord * G.dy, h.zcoord * G.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))

            G.hertziandipoles.append(h)

    # Magnetic dipole
    cmdname = '#magnetic_dipole'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) < 5:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least five parameters')

            # Check polarity & position parameters
            polarisation = tmp[0].lower()
            if polarisation not in ('x', 'y', 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
            if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
            elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
            elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')

            xcoord = G.calculate_coord('x', tmp[1])
            ycoord = G.calculate_coord('y', tmp[2])
            zcoord = G.calculate_coord('z', tmp[3])
            check_coordinates(xcoord, ycoord, zcoord)
            if xcoord < G.pmlthickness['x0'] or xcoord > G.nx - G.pmlthickness['xmax'] or ycoord < G.pmlthickness['y0'] or ycoord > G.ny - G.pmlthickness['ymax'] or zcoord < G.pmlthickness['z0'] or zcoord > G.nz - G.pmlthickness['zmax']:
                print(Fore.RED + "WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.' + Style.RESET_ALL)

            # Check if there is a waveformID in the waveforms list
            if not any(x.ID == tmp[4] for x in G.waveforms):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' there is no waveform with the identifier {}'.format(tmp[4]))

            m = MagneticDipole()
            m.polarisation = polarisation
            m.xcoord = xcoord
            m.ycoord = ycoord
            m.zcoord = zcoord
            m.xcoordorigin = xcoord
            m.ycoordorigin = ycoord
            m.zcoordorigin = zcoord
            m.ID = m.__class__.__name__ + '(' + str(m.xcoord) + ',' + str(m.ycoord) + ',' + str(m.zcoord) + ')'
            m.waveformID = tmp[4]

            if len(tmp) > 5:
                # Check source start & source remove time parameters
                start = float(tmp[5])
                stop = float(tmp[6])
                if start < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' delay of the initiation of the source should not be less than zero')
                if stop < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time to remove the source should not be less than zero')
                if stop - start <= 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' duration of the source should not be zero or less')
                m.start = start
                if stop > G.timewindow:
                    m.stop = G.timewindow
                else:
                    m.stop = stop
                startstop = ' start time {:g} secs, finish time {:g} secs '.format(m.start, m.stop)
            else:
                m.start = 0
                m.stop = G.timewindow
                startstop = ' '

            m.calculate_waveform_values(G)

            if G.messages:
                print('Magnetic dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(m.polarisation, m.xcoord * G.dx, m.ycoord * G.dy, m.zcoord * G.dz) + startstop + 'using waveform {} created.'.format(m.waveformID))

            G.magneticdipoles.append(m)

    # Transmission line
    cmdname = '#transmission_line'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) < 6:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least six parameters')

            # Warn about using a transmission line on GPU
            if G.gpu is not None:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' A #transmission_line cannot currently be used with GPU solving. Consider using a #voltage_source instead.')

            # Check polarity & position parameters
            polarisation = tmp[0].lower()
            if polarisation not in ('x', 'y', 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
            if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
            elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
            elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')

            xcoord = G.calculate_coord('x', tmp[1])
            ycoord = G.calculate_coord('y', tmp[2])
            zcoord = G.calculate_coord('z', tmp[3])
            resistance = float(tmp[4])

            check_coordinates(xcoord, ycoord, zcoord)
            if xcoord < G.pmlthickness['x0'] or xcoord > G.nx - G.pmlthickness['xmax'] or ycoord < G.pmlthickness['y0'] or ycoord > G.ny - G.pmlthickness['ymax'] or zcoord < G.pmlthickness['z0'] or zcoord > G.nz - G.pmlthickness['zmax']:
                print(Fore.RED + "WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.' + Style.RESET_ALL)
            if resistance <= 0 or resistance >= z0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' 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 == tmp[5] for x in G.waveforms):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' there is no waveform with the identifier {}'.format(tmp[5]))

            t = TransmissionLine(G)
            t.polarisation = polarisation
            t.xcoord = xcoord
            t.ycoord = ycoord
            t.zcoord = zcoord
            t.ID = t.__class__.__name__ + '(' + str(t.xcoord) + ',' + str(t.ycoord) + ',' + str(t.zcoord) + ')'
            t.resistance = resistance
            t.waveformID = tmp[5]

            if len(tmp) > 6:
                # Check source start & source remove time parameters
                start = float(tmp[6])
                stop = float(tmp[7])
                if start < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' delay of the initiation of the source should not be less than zero')
                if stop < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time to remove the source should not be less than zero')
                if stop - start <= 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' duration of the source should not be zero or less')
                t.start = start
                if stop > G.timewindow:
                    t.stop = G.timewindow
                else:
                    t.stop = stop
                startstop = ' start time {:g} secs, finish time {:g} secs '.format(t.start, t.stop)
            else:
                t.start = 0
                t.stop = G.timewindow
                startstop = ' '

            t.calculate_waveform_values(G)
            t.calculate_incident_V_I(G)

            if G.messages:
                print('Transmission line with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(t.polarisation, t.xcoord * G.dx, t.ycoord * G.dy, t.zcoord * G.dz, t.resistance) + startstop + 'using waveform {} created.'.format(t.waveformID))

            G.transmissionlines.append(t)

    # Receiver
    cmdname = '#rx'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 3 and len(tmp) < 5:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' has an incorrect number of parameters')

            # Check position parameters
            xcoord = round_value(float(tmp[0]) / G.dx)
            ycoord = round_value(float(tmp[1]) / G.dy)
            zcoord = round_value(float(tmp[2]) / G.dz)
            check_coordinates(xcoord, ycoord, zcoord)
            if xcoord < G.pmlthickness['x0'] or xcoord > G.nx - G.pmlthickness['xmax'] or ycoord < G.pmlthickness['y0'] or ycoord > G.ny - G.pmlthickness['ymax'] or zcoord < G.pmlthickness['z0'] or zcoord > G.nz - G.pmlthickness['zmax']:
                print(Fore.RED + "WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.' + Style.RESET_ALL)

            r = Rx()
            r.xcoord = xcoord
            r.ycoord = ycoord
            r.zcoord = zcoord
            r.xcoordorigin = xcoord
            r.ycoordorigin = ycoord
            r.zcoordorigin = zcoord

            # If no ID or outputs are specified, use default
            if len(tmp) == 3:
                r.ID = r.__class__.__name__ + '(' + str(r.xcoord) + ',' + str(r.ycoord) + ',' + str(r.zcoord) + ')'
                for key in Rx.defaultoutputs:
                    r.outputs[key] = np.zeros(G.iterations, dtype=floattype)
            else:
                r.ID = tmp[3]
                # Get allowable outputs
                if G.gpu is not None:
                    allowableoutputs = Rx.gpu_allowableoutputs
                else:
                    allowableoutputs = Rx.allowableoutputs
                # Check and add field output names
                for field in tmp[4::]:
                    if field in allowableoutputs:
                        r.outputs[field] = np.zeros(G.iterations, dtype=floattype)
                    else:
                        raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' contains an output type that is not allowable. Allowable outputs in current context are {}'.format(allowableoutputs))

            if G.messages:
                print('Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs)))

            G.rxs.append(r)

    # Receiver array
    cmdname = '#rx_array'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 9:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly nine parameters')

            xs = G.calculate_coord('x', tmp[0])
            ys = G.calculate_coord('y', tmp[1])
            zs = G.calculate_coord('z', tmp[2])

            xf = G.calculate_coord('x', tmp[3])
            yf = G.calculate_coord('y', tmp[4])
            zf = G.calculate_coord('z', tmp[5])

            dx = G.calculate_coord('x', tmp[6])
            dy = G.calculate_coord('y', tmp[7])
            dz = G.calculate_coord('z', tmp[8])

            check_coordinates(xs, ys, zs, name='lower')
            check_coordinates(xf, yf, zf, name='upper')

            if xcoord < G.pmlthickness['x0'] or xcoord > G.nx - G.pmlthickness['xmax'] or ycoord < G.pmlthickness['y0'] or ycoord > G.ny - G.pmlthickness['ymax'] or zcoord < G.pmlthickness['z0'] or zcoord > G.nz - G.pmlthickness['zmax']:
                print(Fore.RED + "WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.' + Style.RESET_ALL)
            if xs > xf or ys > yf or zs > zf:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')
            if dx < 0 or dy < 0 or dz < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than zero')
            if dx < 1:
                if dx == 0:
                    dx = 1
                else:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than the spatial discretisation')
            if dy < 1:
                if dy == 0:
                    dy = 1
                else:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than the spatial discretisation')
            if dz < 1:
                if dz == 0:
                    dz = 1
                else:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than the spatial discretisation')

            if G.messages:
                print('Receiver array {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g}m'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, dx * G.dx, dy * G.dy, dz * G.dz))

            for x in range(xs, xf + 1, dx):
                for y in range(ys, yf + 1, dy):
                    for z in range(zs, zf + 1, dz):
                        r = Rx()
                        r.xcoord = x
                        r.ycoord = y
                        r.zcoord = z
                        r.xcoordorigin = x
                        r.ycoordorigin = y
                        r.zcoordorigin = z
                        r.ID = r.__class__.__name__ + '(' + str(x) + ',' + str(y) + ',' + str(z) + ')'
                        for key in Rx.defaultoutputs:
                            r.outputs[key] = np.zeros(G.iterations, dtype=floattype)
                        if G.messages:
                            print('  Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs)))
                        G.rxs.append(r)

    # Snapshot
    cmdname = '#snapshot'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 11:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly eleven parameters')

            xs = G.calculate_coord('x', tmp[0])
            ys = G.calculate_coord('y', tmp[1])
            zs = G.calculate_coord('z', tmp[2])

            xf = G.calculate_coord('x', tmp[3])
            yf = G.calculate_coord('y', tmp[4])
            zf = G.calculate_coord('z', tmp[5])

            dx = G.calculate_coord('x', tmp[6])
            dy = G.calculate_coord('y', tmp[7])
            dz = G.calculate_coord('z', tmp[8])

            # If number of iterations given
            try:
                time = int(tmp[9])
            # If real floating point value given
            except ValueError:
                time = float(tmp[9])
                if time > 0:
                    time = round_value((time / G.dt)) + 1
                else:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time value must be greater than zero')

            check_coordinates(xs, ys, zs, name='lower')
            check_coordinates(xf, yf, zf, name='upper')

            if xs >= xf or ys >= yf or zs >= zf:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')
            if dx < 0 or dy < 0 or dz < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than zero')
            if dx < 1 or dy < 1 or dz < 1:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than the spatial discretisation')
            if time <= 0 or time > G.iterations:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time value is not valid')

            s = Snapshot(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, tmp[10])

            if G.messages:
                print('Snapshot from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, at {:g} secs with filename {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, dx * G.dx, dy * G.dy, dz * G.dz, s.time * G.dt, s.basefilename))

            G.snapshots.append(s)

    # Materials
    cmdname = '#material'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 5:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly five parameters')
            if float(tmp[0]) < 1:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value of one or greater for static (DC) permittivity')
            if tmp[1] != 'inf':
                se = float(tmp[1])
                if se < 0:
                    raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for conductivity')
            else:
                se = float('inf')
            if float(tmp[2]) < 1:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value of one or greater for permeability')
            if float(tmp[3]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for magnetic conductivity')
            if any(x.ID == tmp[4] for x in G.materials):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' with ID {} already exists'.format(tmp[4]))

            # Create a new instance of the Material class material (start index after pec & free_space)
            m = Material(len(G.materials), tmp[4])
            m.er = float(tmp[0])
            m.se = se
            m.mr = float(tmp[2])
            m.sm = float(tmp[3])

            # Set material averaging to False if infinite conductivity, i.e. pec
            if m.se == float('inf'):
                m.averagable = False

            if G.messages:
                tqdm.write('Material {} with eps_r={:g}, sigma={:g} S/m; mu_r={:g}, sigma*={:g} Ohm/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm))

            # Append the new material object to the materials list
            G.materials.append(m)

    cmdname = '#add_dispersion_debye'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()

            if len(tmp) < 4:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least four parameters')
            if int(tmp[0]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for number of poles')
            poles = int(tmp[0])
            materialsrequested = tmp[(2 * poles) + 1:len(tmp)]

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            for material in materials:
                material.type = 'debye'
                material.poles = poles
                material.averagable = False
                for pole in range(1, 2 * poles, 2):
                    # N.B Not checking if relaxation times are greater than time-step
                    if float(tmp[pole]) > 0:
                        material.deltaer.append(float(tmp[pole]))
                        material.tau.append(float(tmp[pole + 1]))
                    else:
                        raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the permittivity difference.')
                if material.poles > Material.maxpoles:
                    Material.maxpoles = material.poles

                if G.messages:
                    tqdm.write('Debye disperion added to {} with delta_eps_r={}, and tau={} secs created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau)))

    cmdname = '#add_dispersion_lorentz'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()

            if len(tmp) < 5:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least five parameters')
            if int(tmp[0]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for number of poles')
            poles = int(tmp[0])
            materialsrequested = tmp[(3 * poles) + 1:len(tmp)]

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            for material in materials:
                material.type = 'lorentz'
                material.poles = poles
                material.averagable = False
                for pole in range(1, 3 * poles, 3):
                    if float(tmp[pole]) > 0 and float(tmp[pole + 1]) > G.dt and float(tmp[pole + 2]) > G.dt:
                        material.deltaer.append(float(tmp[pole]))
                        material.tau.append(float(tmp[pole + 1]))
                        material.alpha.append(float(tmp[pole + 2]))
                    else:
                        raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the permittivity difference and frequencies, and associated times that are greater than the time step for the model.')
                if material.poles > Material.maxpoles:
                    Material.maxpoles = material.poles

                if G.messages:
                    tqdm.write('Lorentz disperion added to {} with delta_eps_r={}, omega={} secs, and gamma={} created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))

    cmdname = '#add_dispersion_drude'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()

            if len(tmp) < 5:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least five parameters')
            if int(tmp[0]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for number of poles')
            poles = int(tmp[0])
            materialsrequested = tmp[(3 * poles) + 1:len(tmp)]

            # Look up requested materials in existing list of material instances
            materials = [y for x in materialsrequested for y in G.materials if y.ID == x]

            if len(materials) != len(materialsrequested):
                notfound = [x for x in materialsrequested if x not in materials]
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound))

            for material in materials:
                material.type = 'drude'
                material.poles = poles
                material.averagable = False
                for pole in range(1, 2 * poles, 2):
                    if float(tmp[pole]) > 0 and float(tmp[pole + 1]) > G.dt:
                        material.tau.append(float(tmp[pole]))
                        material.alpha.append(float(tmp[pole + 1]))
                    else:
                        raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the frequencies, and associated times that are greater than the time step for the model.')
                if material.poles > Material.maxpoles:
                    Material.maxpoles = material.poles

                if G.messages:
                    tqdm.write('Drude disperion added to {} with omega={} secs, and gamma={} secs created.'.format(material.ID, ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))

    cmdname = '#soil_peplinski'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 7:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at exactly seven parameters')
            if float(tmp[0]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for the sand fraction')
            if float(tmp[1]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for the clay fraction')
            if float(tmp[2]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for the bulk density')
            if float(tmp[3]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for the sand particle density')
            if float(tmp[4]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for the lower limit of the water volumetric fraction')
            if float(tmp[5]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for the upper limit of the water volumetric fraction')
            if any(x.ID == tmp[6] for x in G.mixingmodels):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' with ID {} already exists'.format(tmp[6]))

            # Create a new instance of the Material class material (start index after pec & free_space)
            s = PeplinskiSoil(tmp[6], float(tmp[0]), float(tmp[1]), float(tmp[2]), float(tmp[3]), (float(tmp[4]), float(tmp[5])))

            if G.messages:
                print('Mixing model (Peplinski) used to create {} with sand fraction {:g}, clay fraction {:g}, bulk density {:g}g/cm3, sand particle density {:g}g/cm3, and water volumetric fraction {:g} to {:g} created.'.format(s.ID, s.S, s.C, s.rb, s.rs, s.mu[0], s.mu[1]))

            # Append the new material object to the materials list
            G.mixingmodels.append(s)

    # Geometry views (creates VTK-based geometry files)
    cmdname = '#geometry_view'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 11:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly eleven parameters')

            xs = G.calculate_coord('x', tmp[0])
            ys = G.calculate_coord('y', tmp[1])
            zs = G.calculate_coord('z', tmp[2])

            xf = G.calculate_coord('x', tmp[3])
            yf = G.calculate_coord('y', tmp[4])
            zf = G.calculate_coord('z', tmp[5])

            dx = G.calculate_coord('x', tmp[6])
            dy = G.calculate_coord('y', tmp[7])
            dz = G.calculate_coord('z', tmp[8])

            check_coordinates(xs, ys, zs, name='lower')
            check_coordinates(xf, yf, zf, name='upper')

            if xs >= xf or ys >= yf or zs >= zf:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')
            if dx < 0 or dy < 0 or dz < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than zero')
            if dx > G.nx or dy > G.ny or dz > G.nz:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should be less than the domain size')
            if dx < 1 or dy < 1 or dz < 1:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than the spatial discretisation')
            if tmp[10].lower() != 'n' and tmp[10].lower() != 'f':
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires type to be either n (normal) or f (fine)')
            if tmp[10].lower() == 'f' and (dx * G.dx != G.dx or dy * G.dy != G.dy or dz * G.dz != G.dz):
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires the spatial discretisation for the geometry view to be the same as the model for geometry view of type f (fine)')

            # Set type of geometry file
            if tmp[10].lower() == 'n':
                fileext = '.vti'
            else:
                fileext = '.vtp'

            g = GeometryView(xs, ys, zs, xf, yf, zf, dx, dy, dz, tmp[9], fileext)

            if G.messages:
                print('Geometry view from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, with filename base {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, dx * G.dx, dy * G.dy, dz * G.dz, g.basefilename))

            # Append the new GeometryView object to the geometry views list
            G.geometryviews.append(g)

    # Geometry object(s) output
    cmdname = '#geometry_objects_write'
    if multicmds[cmdname] is not None:
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 7:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly seven parameters')

            xs = G.calculate_coord('x', tmp[0])
            ys = G.calculate_coord('y', tmp[1])
            zs = G.calculate_coord('z', tmp[2])

            xf = G.calculate_coord('x', tmp[3])
            yf = G.calculate_coord('y', tmp[4])
            zf = G.calculate_coord('z', tmp[5])

            check_coordinates(xs, ys, zs, name='lower')
            check_coordinates(xf, yf, zf, name='upper')

            if xs >= xf or ys >= yf or zs >= zf:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the lower coordinates should be less than the upper coordinates')

            g = GeometryObjects(xs, ys, zs, xf, yf, zf, tmp[6])

            if G.messages:
                print('Geometry objects in the volume from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, will be written to {}, with materials written to {}'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, g.filename, g.materialsfilename))

            # Append the new GeometryView object to the geometry objects to write list
            G.geometryobjectswrite.append(g)

    # Complex frequency shifted (CFS) PML parameter
    cmdname = '#pml_cfs'
    if multicmds[cmdname] is not None:
        if len(multicmds[cmdname]) > 2:
            raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' can only be used up to two times, for up to a 2nd order PML')
        for cmdinstance in multicmds[cmdname]:
            tmp = cmdinstance.split()
            if len(tmp) != 12:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly twelve parameters')
            if tmp[0] not in CFSParameter.scalingprofiles.keys() or tmp[4] not in CFSParameter.scalingprofiles.keys() or tmp[8] not in CFSParameter.scalingprofiles.keys():
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' must have scaling type {}'.format(','.join(CFSParameter.scalingprofiles.keys())))
            if tmp[1] not in CFSParameter.scalingdirections or tmp[5] not in CFSParameter.scalingdirections or tmp[9] not in CFSParameter.scalingdirections:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' must have scaling type {}'.format(','.join(CFSParameter.scalingdirections)))
            if float(tmp[2]) < 0 or float(tmp[3]) < 0 or float(tmp[6]) < 0 or float(tmp[7]) < 0 or float(tmp[10]) < 0:
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' minimum and maximum scaling values must be positive')
            if float(tmp[6]) < 1 and G.pmlformulation == 'HORIPML':
                raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' minimum scaling value for kappa must be greater than or equal to one')

            cfsalpha = CFSParameter()
            cfsalpha.ID = 'alpha'
            cfsalpha.scalingprofile = tmp[0]
            cfsalpha.scalingdirection = tmp[1]
            cfsalpha.min = float(tmp[2])
            cfsalpha.max = float(tmp[3])
            cfskappa = CFSParameter()
            cfskappa.ID = 'kappa'
            cfskappa.scalingprofile = tmp[4]
            cfskappa.scalingdirection = tmp[5]
            cfskappa.min = float(tmp[6])
            cfskappa.max = float(tmp[7])
            cfssigma = CFSParameter()
            cfssigma.ID = 'sigma'
            cfssigma.scalingprofile = tmp[8]
            cfssigma.scalingdirection = tmp[9]
            cfssigma.min = float(tmp[10])
            if tmp[11] == 'None':
                cfssigma.max = None
            else:
                cfssigma.max = float(tmp[11])
            cfs = CFS()
            cfs.alpha = cfsalpha
            cfs.kappa = cfskappa
            cfs.sigma = cfssigma

            if G.messages:
                print('PML CFS parameters: alpha (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), kappa (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), sigma (scaling: {}, scaling direction: {}, min: {:g}, max: {}) created.'.format(cfsalpha.scalingprofile, cfsalpha.scalingdirection, cfsalpha.min, cfsalpha.max, cfskappa.scalingprofile, cfskappa.scalingdirection, cfskappa.min, cfskappa.max, cfssigma.scalingprofile, cfssigma.scalingdirection, cfssigma.min, cfssigma.max))

            G.cfs.append(cfs) def process_singlecmds(singlecmds, G):
    """Checks the validity of command parameters and creates instances of classes of parameters.

    Args:
        singlecmds (dict): Commands that can only occur once in the model.
        G (class): Grid class instance - holds essential parameters describing the model.
    """

    # Check validity of command parameters in order needed
    # messages
    cmd = '#messages'
    if singlecmds[cmd] is not None:
        tmp = singlecmds[cmd].split()
        if len(tmp) != 1:
            raise CmdInputError(cmd + ' requires exactly one parameter')
        if singlecmds[cmd].lower() == 'y':
            G.messages = True
        elif singlecmds[cmd].lower() == 'n':
            G.messages = False
        else:
            raise CmdInputError(cmd + ' requires input values of either y or n')

    # Title
    cmd = '#title'
    if singlecmds[cmd] is not None:
        G.title = singlecmds[cmd]
        if G.messages:
            print('Model title: {}'.format(G.title))

    # Get information about host machine
    hostinfo = get_host_info()

    # Number of threads (OpenMP) to use
    cmd = '#num_threads'
    if sys.platform == 'darwin':
        os.environ['OMP_WAIT_POLICY'] = 'ACTIVE'  # Should waiting threads consume CPU power (can drastically effect performance)
    os.environ['OMP_DYNAMIC'] = 'FALSE'  # Number of threads may be adjusted by the run time environment to best utilize system resources
    os.environ['OMP_PLACES'] = 'cores'  # Each place corresponds to a single core (having one or more hardware threads)
    os.environ['OMP_PROC_BIND'] = 'TRUE'  # Bind threads to physical cores
    # os.environ['OMP_DISPLAY_ENV'] = 'TRUE' # Prints OMP version and environment variables (useful for debug)

    # Catch bug with Windows Subsystem for Linux (https://github.com/Microsoft/BashOnWindows/issues/785)
    if 'Microsoft' in G.hostinfo['osversion']:
        os.environ['KMP_AFFINITY'] = 'disabled'
        del os.environ['OMP_PLACES']
        del os.environ['OMP_PROC_BIND']

    if singlecmds[cmd] is not None:
        tmp = tuple(int(x) for x in singlecmds[cmd].split())
        if len(tmp) != 1:
            raise CmdInputError(cmd + ' requires exactly one parameter to specify the number of threads to use')
        if tmp[0] < 1:
            raise CmdInputError(cmd + ' requires the value to be an integer not less than one')
        G.nthreads = tmp[0]
        os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
    elif os.environ.get('OMP_NUM_THREADS'):
        G.nthreads = int(os.environ.get('OMP_NUM_THREADS'))
    else:
        # Set number of threads to number of physical CPU cores
        G.nthreads = hostinfo['physicalcores']
        os.environ['OMP_NUM_THREADS'] = str(G.nthreads)

    if G.messages:
        print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
    if G.nthreads > G.hostinfo['physicalcores']:
        print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)

    # Print information about any GPU in use
    if G.messages:
        if G.gpu is not None:
            print('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name))

    # Spatial discretisation
    cmd = '#dx_dy_dz'
    tmp = [float(x) for x in singlecmds[cmd].split()]
    if len(tmp) != 3:
        raise CmdInputError(cmd + ' requires exactly three parameters')
    if tmp[0] <= 0:
        raise CmdInputError(cmd + ' requires the x-direction spatial step to be greater than zero')
    if tmp[1] <= 0:
        raise CmdInputError(cmd + ' requires the y-direction spatial step to be greater than zero')
    if tmp[2] <= 0:
        raise CmdInputError(cmd + ' requires the z-direction spatial step to be greater than zero')
    G.dx = tmp[0]
    G.dy = tmp[1]
    G.dz = tmp[2]
    if G.messages:
        print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(G.dx, G.dy, G.dz))

    # Domain
    cmd = '#domain'
    tmp = [float(x) for x in singlecmds[cmd].split()]
    if len(tmp) != 3:
        raise CmdInputError(cmd + ' requires exactly three parameters')
    G.nx = round_value(tmp[0] / G.dx)
    G.ny = round_value(tmp[1] / G.dy)
    G.nz = round_value(tmp[2] / G.dz)
    if G.nx == 0 or G.ny == 0 or G.nz == 0:
        raise CmdInputError(cmd + ' requires at least one cell in every dimension')
    if G.messages:
        print('Domain size: {:g} x {:g} x {:g}m ({:d} x {:d} x {:d} = {:g} cells)'.format(tmp[0], tmp[1], tmp[2], G.nx, G.ny, G.nz, (G.nx * G.ny * G.nz)))

    # Time step CFL limit (either 2D or 3D); switch off appropriate PMLs for 2D
    if G.nx == 1:
        G.dt = 1 / (c * np.sqrt((1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
        G.mode = '2D TMx'
        G.pmlthickness['x0'] = 0
        G.pmlthickness['xmax'] = 0
    elif G.ny == 1:
        G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dz) * (1 / G.dz)))
        G.mode = '2D TMy'
        G.pmlthickness['y0'] = 0
        G.pmlthickness['ymax'] = 0
    elif G.nz == 1:
        G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy)))
        G.mode = '2D TMz'
        G.pmlthickness['z0'] = 0
        G.pmlthickness['zmax'] = 0
    else:
        G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
        G.mode = '3D'

    # Round down time step to nearest float with precision one less than hardware maximum.
    # Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
    G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1)

    if G.messages:
        print('Mode: {}'.format(G.mode))
        print('Time step (at CFL limit): {:g} secs'.format(G.dt))

    # Time step stability factor
    cmd = '#time_step_stability_factor'
    if singlecmds[cmd] is not None:
        tmp = tuple(float(x) for x in singlecmds[cmd].split())
        if len(tmp) != 1:
            raise CmdInputError(cmd + ' requires exactly one parameter')
        if tmp[0] <= 0 or tmp[0] > 1:
            raise CmdInputError(cmd + ' requires the value of the time step stability factor to be between zero and one')
        G.dt = G.dt * tmp[0]
        if G.messages:
            print('Time step (modified): {:g} secs'.format(G.dt))

    # Time window
    cmd = '#time_window'
    tmp = singlecmds[cmd].split()
    if len(tmp) != 1:
        raise CmdInputError(cmd + ' requires exactly one parameter to specify the time window. Either in seconds or number of iterations.')
    tmp = tmp[0].lower()

    # If number of iterations given
    # The +/- 1 used in calculating the number of iterations is to account for
    # the fact that the solver (iterations) loop runs from 0 to < G.iterations
    try:
        tmp = int(tmp)
        G.timewindow = (tmp - 1) * G.dt
        G.iterations = tmp
    # If real floating point value given
    except ValueError:
        tmp = float(tmp)
        if tmp > 0:
            G.timewindow = tmp
            G.iterations = int(np.ceil(tmp / G.dt)) + 1
        else:
            raise CmdInputError(cmd + ' must have a value greater than zero')
    if G.messages:
        print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))

    # PML cells
    cmd = '#pml_cells'
    if singlecmds[cmd] is not None:
        tmp = singlecmds[cmd].split()
        if len(tmp) != 1 and len(tmp) != 6:
            raise CmdInputError(cmd + ' requires either one or six parameter(s)')
        if len(tmp) == 1:
            for key in G.pmlthickness.keys():
                G.pmlthickness[key] = int(tmp[0])
        else:
            G.pmlthickness['x0'] = int(tmp[0])
            G.pmlthickness['y0'] = int(tmp[1])
            G.pmlthickness['z0'] = int(tmp[2])
            G.pmlthickness['xmax'] = int(tmp[3])
            G.pmlthickness['ymax'] = int(tmp[4])
            G.pmlthickness['zmax'] = int(tmp[5])
    if 2 * G.pmlthickness['x0'] >= G.nx or 2 * G.pmlthickness['y0'] >= G.ny or 2 * G.pmlthickness['z0'] >= G.nz or 2 * G.pmlthickness['xmax'] >= G.nx or 2 * G.pmlthickness['ymax'] >= G.ny or 2 * G.pmlthickness['zmax'] >= G.nz:
        raise CmdInputError(cmd + ' has too many cells for the domain size')

    # PML formulation
    cmd = '#pml_formulation'
    if singlecmds[cmd] is not None:
        tmp = singlecmds[cmd].split()
        if len(tmp) != 1:
            raise CmdInputError(cmd + ' requires exactly one parameter')
        if singlecmds[cmd].upper() in PML.formulations:
            G.pmlformulation = singlecmds[cmd].upper()
        else:
            raise CmdInputError(cmd + ' PML formulation is not found')

    # src_steps
    cmd = '#src_steps'
    if singlecmds[cmd] is not None:
        tmp = singlecmds[cmd].split()
        if len(tmp) != 3:
            raise CmdInputError(cmd + ' requires exactly three parameters')
        G.srcsteps[0] = round_value(float(tmp[0]) / G.dx)
        G.srcsteps[1] = round_value(float(tmp[1]) / G.dy)
        G.srcsteps[2] = round_value(float(tmp[2]) / G.dz)
        if G.messages:
            print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcsteps[0] * G.dx, G.srcsteps[1] * G.dy, G.srcsteps[2] * G.dz))

    # rx_steps
    cmd = '#rx_steps'
    if singlecmds[cmd] is not None:
        tmp = singlecmds[cmd].split()
        if len(tmp) != 3:
            raise CmdInputError(cmd + ' requires exactly three parameters')
        G.rxsteps[0] = round_value(float(tmp[0]) / G.dx)
        G.rxsteps[1] = round_value(float(tmp[1]) / G.dy)
        G.rxsteps[2] = round_value(float(tmp[2]) / G.dz)
        if G.messages:
            print('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxsteps[0] * G.dx, G.rxsteps[1] * G.dy, G.rxsteps[2] * G.dz))

    # Excitation file for user-defined source waveforms
    cmd = '#excitation_file'
    if singlecmds[cmd] is not None:
        tmp = singlecmds[cmd].split()
        if len(tmp) != 1 and len(tmp) != 3:
            raise CmdInputError(cmd + ' requires either one or three parameter(s)')
        excitationfile = tmp[0]

        # Optional parameters passed directly to scipy.interpolate.interp1d
        kwargs = dict()
        if len(tmp) > 1:
            kwargs['kind'] = tmp[1]
            kwargs['fill_value'] = tmp[2]
        else:
            args, varargs, keywords, defaults = inspect.getargspec(interpolate.interp1d)
            kwargs = dict(zip(reversed(args), reversed(defaults)))

        # 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))

        if G.messages:
            print('\nExcitation file: {}'.format(excitationfile))

        # 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)

        # 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('Waveform with ID {} already exists'.format(waveformIDs[waveform]))
            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]

            # 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):
                tmp = np.zeros(len(waveformtime))
                tmp[:len(singlewaveformvalues)] = singlewaveformvalues
                singlewaveformvalues = tmp

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

            if G.messages:
                print('User waveform {} created using {} and, if required, interpolation parameters (kind: {}, fill value: {}).'.format(w.ID, timestr, kwargs['kind'], kwargs['fill_value']))

            G.waveforms.append(w)

    # Set the output directory
    cmd = '#output_dir'
    if singlecmds[cmd] is not None:
        outputdir = singlecmds[cmd]
        G.outputdirectory = outputdir class Material(object):
    """Materials, their properties and update coefficients."""

    # Maximum number of dispersive material poles in a model
    maxpoles = 0

    # Properties of water from: http://dx.doi.org/10.1109/TGRS.2006.873208
    waterer = 80.1
    watereri = 4.9
    waterdeltaer = waterer - watereri
    watertau = 9.231e-12

    # Properties of grass from: http://dx.doi.org/10.1007/BF00902994
    grasser = 18.5087
    grasseri = 12.7174
    grassdeltaer = grasser - grasseri
    grasstau = 1.0793e-11

    def __init__(self, numID, ID):
        """
        Args:
            numID (int): Numeric identifier of the material.
            ID (str): Name of the material.
        """

        self.numID = numID
        self.ID = ID
        self.type = ''
        # Default material averaging
        self.averagable = True

        # Default material constitutive parameters (free_space)
        self.er = 1.0
        self.se = 0.0
        self.mr = 1.0
        self.sm = 0.0

        # Parameters for dispersive materials
        self.poles = 0
        self.deltaer = []
        self.tau = []
        self.alpha = []

    def calculate_update_coeffsH(self, G):
        """Calculates the magnetic update coefficients of the material.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        HA = (m0 * self.mr / G.dt) + 0.5 * self.sm
        HB = (m0 * self.mr / G.dt) - 0.5 * self.sm
        self.DA = HB / HA
        self.DBx = (1 / G.dx) * 1 / HA
        self.DBy = (1 / G.dy) * 1 / HA
        self.DBz = (1 / G.dz) * 1 / HA
        self.srcm = 1 / HA

    def calculate_update_coeffsE(self, G):
        """Calculates the electric update coefficients of the material.

        Args:
            G (class): Grid class instance - holds essential parameters
                    describing the model.
        """

        # The implementation of the dispersive material modelling comes from the
        # derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549
        if self.maxpoles > 0:
            self.w = np.zeros(self.maxpoles, dtype=complextype)
            self.q = np.zeros(self.maxpoles, dtype=complextype)
            self.zt = np.zeros(self.maxpoles, dtype=complextype)
            self.zt2 = np.zeros(self.maxpoles, dtype=complextype)
            self.eqt = np.zeros(self.maxpoles, dtype=complextype)
            self.eqt2 = np.zeros(self.maxpoles, dtype=complextype)

            for x in range(self.poles):
                if 'debye' in self.type:
                    self.w[x] = self.deltaer[x] / self.tau[x]
                    self.q[x] = -1 / self.tau[x]
                elif 'lorentz' in self.type:
                    # tau for Lorentz materials are pole frequencies
                    # alpha for Lorentz materials are the damping coefficients
                    wp2 = (2 * np.pi * self.tau[x])**2
                    self.w[x] = -1j * ((wp2 * self.deltaer[x]) / np.sqrt(wp2 - self.alpha[x]**2))
                    self.q[x] = -self.alpha[x] + (1j * np.sqrt(wp2 - self.alpha[x]**2))
                elif 'drude' in self.type:
                    # tau for Drude materials are pole frequencies
                    # alpha for Drude materials are the inverse of relaxation times
                    wp2 = (2 * np.pi * self.tau[x])**2
                    self.se += wp2 / self.alpha[x]
                    self.w[x] = - (wp2 / self.alpha[x])
                    self.q[x] = - self.alpha[x]

                self.eqt[x] = np.exp(self.q[x] * G.dt)
                self.eqt2[x] = np.exp(self.q[x] * (G.dt / 2))
                self.zt[x] = (self.w[x] / self.q[x]) * (1 - self.eqt[x]) / G.dt
                self.zt2[x] = (self.w[x] / self.q[x]) * (1 - self.eqt2[x])

            EA = (e0 * self.er / G.dt) + 0.5 * self.se - (e0 / G.dt) * np.sum(self.zt2.real)
            EB = (e0 * self.er / G.dt) - 0.5 * self.se - (e0 / G.dt) * np.sum(self.zt2.real)

        else:
            EA = (e0 * self.er / G.dt) + 0.5 * self.se
            EB = (e0 * self.er / G.dt) - 0.5 * self.se

        if self.ID == 'pec' or self.se == float('inf'):
            self.CA = 0
            self.CBx = 0
            self.CBy = 0
            self.CBz = 0
            self.srce = 0
        else:
            self.CA = EB / EA
            self.CBx = (1 / G.dx) * 1 / EA
            self.CBy = (1 / G.dy) * 1 / EA
            self.CBz = (1 / G.dz) * 1 / EA
            self.srce = 1 / EA

    def calculate_er(self, freq):
        """
        Calculates the complex relative permittivity of the material at a specific frequency.

        Args:
            freq (float): Frequency used to calculate complex relative permittivity.

        Returns:
            er (float): Complex relative permittivity.
        """

        # Permittivity at infinite frequency if the material is dispersive
        er = self.er

        if self.poles > 0:
            w = 2 * np.pi * freq
            er += self.se / (1j * w * e0)
            if 'debye' in self.type:
                for pole in range(self.poles):
                    er += self.deltaer[pole] / (1 + 1j * w * self.tau[pole])
            elif 'lorentz' in self.type:
                for pole in range(self.poles):
                    er += (self.deltaer[pole] * self.tau[pole]**2) / (self.tau[pole]**2 + 2j * w * self.alpha[pole] - w**2)
            elif 'drude' in self.type:
                ersum = 0
                for pole in range(self.poles):
                    ersum += self.tau[pole]**2 / (w**2 - 1j * w * self.alpha[pole])
                    er -= ersum

        return er def process_materials(G):
    """
    Process complete list of materials - calculate update coefficients,
        store in arrays, and build text list of materials/properties

    Args:
        G (class): Grid class instance - holds essential parameters describing the model.

    Returns:
        materialsdata (list): List of material IDs, names, and properties to print a table.
    """

    if Material.maxpoles == 0:
        materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]', '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
    else:
        materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]', 'Delta\neps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]', 'gamma\n[Hz]', '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]

    for material in G.materials:
        # Calculate update coefficients for material
        material.calculate_update_coeffsE(G)
        material.calculate_update_coeffsH(G)

        # Store all update coefficients together
        G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
        G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm

        # Store coefficients for any dispersive materials
        if Material.maxpoles > 0:
            z = 0
            for pole in range(Material.maxpoles):
                G.updatecoeffsdispersive[material.numID, z:z + 3] = e0 * material.eqt2[pole], material.eqt[pole], material.zt[pole]
                z += 3

        # Construct information on material properties for printing table
        materialtext = []
        materialtext.append(str(material.numID))
        materialtext.append(material.ID[:50] if len(material.ID) > 50 else material.ID)
        materialtext.append(material.type)
        materialtext.append('{:g}'.format(material.er))
        materialtext.append('{:g}'.format(material.se))
        if Material.maxpoles > 0:
            if 'debye' in material.type:
                materialtext.append('\n'.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
                materialtext.append('\n'.join('{:g}'.format(tau) for tau in material.tau))
                materialtext.extend(['', '', ''])
            elif 'lorentz' in material.type:
                materialtext.append(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
                materialtext.append('')
                materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
                materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha))
                materialtext.append('')
            elif 'drude' in material.type:
                materialtext.extend(['', ''])
                materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
                materialtext.append('')
                materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha))
            else:
                materialtext.extend(['', '', '', '', ''])

        materialtext.append('{:g}'.format(material.mr))
        materialtext.append('{:g}'.format(material.sm))
        materialtext.append(material.averagable)
        materialsdata.append(materialtext)

    return materialsdata class CFS(object):
    """CFS term for PML."""

    def __init__(self):
        """
        Args:
            alpha (CFSParameter): alpha parameter for CFS.
            kappa (CFSParameter): kappa parameter for CFS.
            sigma (CFSParameter): sigma parameter for CFS.
        """

        self.alpha = CFSParameter(ID='alpha', scalingprofile='constant')
        self.kappa = CFSParameter(ID='kappa', scalingprofile='constant', min=1, max=1)
        self.sigma = CFSParameter(ID='sigma', scalingprofile='quartic', min=0, max=None)

    def calculate_sigmamax(self, d, er, mr, G):
        """Calculates an optimum value for sigma max based on underlying material properties.

        Args:
            d (float): dx, dy, or dz in direction of PML.
            er (float): Average permittivity of underlying material.
            mr (float): Average permeability of underlying material.
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        # Calculation of the maximum value of sigma from http://dx.doi.org/10.1109/8.546249
        m = CFSParameter.scalingprofiles[self.sigma.scalingprofile]
        self.sigma.max = (0.8 * (m + 1)) / (z0 * d * np.sqrt(er * mr))

    def scaling_polynomial(self, order, Evalues, Hvalues):
        """Applies the polynomial to be used for the scaling profile for electric and magnetic PML updates.

        Args:
            order (int): Order of polynomial for scaling profile.
            Evalues (float): numpy array holding scaling profile values for electric PML update.
            Hvalues (float): numpy array holding scaling profile values for magnetic PML update.

        Returns:
            Evalues (float): numpy array holding scaling profile values for electric PML update.
            Hvalues (float): numpy array holding scaling profile values for magnetic PML update.
        """

        tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2 * len(Evalues)) / (len(Evalues) - 1)) ** order
        Evalues = tmp[0:-1:2]
        Hvalues = tmp[1::2]

        return Evalues, Hvalues

    def calculate_values(self, thickness, parameter):
        """Calculates values for electric and magnetic PML updates based on profile type and minimum and maximum values.

        Args:
            thickness (int): Thickness of PML in cells.
            parameter (CFSParameter): Instance of CFSParameter

        Returns:
            Evalues (float): numpy array holding profile value for electric PML update.
            Hvalues (float): numpy array holding profile value for magnetic PML update.
        """

        # Extra cell of thickness added to allow correct scaling of electric and magnetic values
        Evalues = np.zeros(thickness + 1, dtype=floattype)
        Hvalues = np.zeros(thickness + 1, dtype=floattype)

        if parameter.scalingprofile == 'constant':
            Evalues += parameter.max
            Hvalues += parameter.max

        elif parameter.scaling == 'polynomial':
            Evalues, Hvalues = self.scaling_polynomial(CFSParameter.scalingprofiles[parameter.scalingprofile], Evalues, Hvalues)
            if parameter.ID == 'alpha':
                Evalues = Evalues * (self.alpha.max - self.alpha.min) + self.alpha.min
                Hvalues = Hvalues * (self.alpha.max - self.alpha.min) + self.alpha.min
            elif parameter.ID == 'kappa':
                Evalues = Evalues * (self.kappa.max - self.kappa.min) + self.kappa.min
                Hvalues = Hvalues * (self.kappa.max - self.kappa.min) + self.kappa.min
            elif parameter.ID == 'sigma':
                Evalues = Evalues * (self.sigma.max - self.sigma.min) + self.sigma.min
                Hvalues = Hvalues * (self.sigma.max - self.sigma.min) + self.sigma.min

        if parameter.scalingdirection == 'reverse':
            Evalues = Evalues[::-1]
            Hvalues = Hvalues[::-1]
            # Magnetic values must be shifted one element to the left after reversal
            Hvalues = np.roll(Hvalues, -1)

        # Extra cell of thickness not required and therefore removed after scaling
        Evalues = Evalues[:-1]
        Hvalues = Hvalues[:-1]

        return Evalues, Hvalues class PML(object):
    """Perfectly Matched Layer (PML) Absorbing Boundary Conditions (ABC)"""

    # Available PML formulations:
    # Higher Order RIPML (HORIPML) see: https://doi.org/10.1109/TAP.2011.2180344
    # Multipole RIPML (MRIPML) see: https://doi.org/10.1109/TAP.2018.2823864
    formulations = ['HORIPML', 'MRIPML']

    # PML slabs IDs at boundaries of domain.
    boundaryIDs = ['x0', 'y0', 'z0', 'xmax', 'ymax', 'zmax']

    # Indicates direction of increasing absorption
    # xminus, yminus, zminus - absorption increases in negative direction of x-axis, y-axis, or z-axis
    # xplus, yplus, zplus - absorption increases in positive direction of x-axis, y-axis, or z-axis
    directions = ['xminus', 'yminus', 'zminus', 'xplus', 'yplus', 'zplus']

    def __init__(self, G, ID=None, direction=None, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0):
        """
        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
            ID (str): Identifier for PML slab.
            direction (str): Direction of increasing absorption.
            xs, xf, ys, yf, zs, zf (float): Extent of the PML slab.
        """

        self.ID = ID
        self.direction = direction
        self.xs = xs
        self.xf = xf
        self.ys = ys
        self.yf = yf
        self.zs = zs
        self.zf = zf
        self.nx = xf - xs
        self.ny = yf - ys
        self.nz = zf - zs

        # Spatial discretisation and thickness
        if self.direction[0] == 'x':
            self.d = G.dx
            self.thickness = self.nx
        elif self.direction[0] == 'y':
            self.d = G.dy
            self.thickness = self.ny
        elif self.direction[0] == 'z':
            self.d = G.dz
            self.thickness = self.nz

        self.CFS = G.cfs

        if G.gpu is None:
            self.initialise_field_arrays()

    def initialise_field_arrays(self):
        """Initialise arrays to store fields in PML."""

        if self.direction[0] == 'x':
            self.EPhi1 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=floattype)
            self.EPhi2 = np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=floattype)
            self.HPhi1 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype)
            self.HPhi2 = np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=floattype)
        elif self.direction[0] == 'y':
            self.EPhi1 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=floattype)
            self.EPhi2 = np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=floattype)
            self.HPhi1 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype)
            self.HPhi2 = np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=floattype)
        elif self.direction[0] == 'z':
            self.EPhi1 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=floattype)
            self.EPhi2 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=floattype)
            self.HPhi1 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype)
            self.HPhi2 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype)

    def calculate_update_coeffs(self, er, mr, G):
        """Calculates electric and magnetic update coefficients for the PML.

        Args:
            er (float): Average permittivity of underlying material
            mr (float): Average permeability of underlying material
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        self.ERA = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
        self.ERB = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
        self.ERE = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
        self.ERF = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
        self.HRA = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
        self.HRB = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
        self.HRE = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
        self.HRF = np.zeros((len(self.CFS), self.thickness), dtype=floattype)

        for x, cfs in enumerate(self.CFS):
            if not cfs.sigma.max:
                cfs.calculate_sigmamax(self.d, er, mr, G)
            Ealpha, Halpha = cfs.calculate_values(self.thickness, cfs.alpha)
            Ekappa, Hkappa = cfs.calculate_values(self.thickness, cfs.kappa)
            Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma)

            # Define different parameters depending on PML formulation
            if G.pmlformulation == 'HORIPML':
                # HORIPML electric update coefficients
                tmp = (2 * e0 * Ekappa) + G.dt * (Ealpha * Ekappa + Esigma)
                self.ERA[x, :] = (2 * e0 + G.dt * Ealpha) / tmp
                self.ERB[x, :] = (2 * e0 * Ekappa) / tmp
                self.ERE[x, :] = ((2 * e0 * Ekappa) - G.dt * (Ealpha * Ekappa + Esigma)) / tmp
                self.ERF[x, :] = (2 * Esigma * G.dt) / (Ekappa * tmp)

                # HORIPML magnetic update coefficients
                tmp = (2 * e0 * Hkappa) + G.dt * (Halpha * Hkappa + Hsigma)
                self.HRA[x, :] = (2 * e0 + G.dt * Halpha) / tmp
                self.HRB[x, :] = (2 * e0 * Hkappa) / tmp
                self.HRE[x, :] = ((2 * e0 * Hkappa) - G.dt * (Halpha * Hkappa + Hsigma)) / tmp
                self.HRF[x, :] = (2 * Hsigma * G.dt) / (Hkappa * tmp)

            elif G.pmlformulation == 'MRIPML':
                tmp = 2 * e0 + G.dt * Ealpha
                self.ERA[x, :] = Ekappa + (G.dt * Esigma) / tmp
                self.ERB[x, :] = (2 * e0) / tmp
                self.ERE[x, :] = ((2 * e0) - G.dt * Ealpha) / tmp
                self.ERF[x, :] = (2 * Esigma * G.dt) / tmp

                # MRIPML magnetic update coefficients
                tmp = 2 * e0 + G.dt * Halpha
                self.HRA[x, :] = Hkappa + (G.dt * Hsigma) / tmp
                self.HRB[x, :] = (2 * e0) / tmp
                self.HRE[x, :] = ((2 * e0) - G.dt * Halpha) / tmp
                self.HRF[x, :] = (2 * Hsigma * G.dt) / tmp

    def update_electric(self, G):
        """This functions updates electric field components with the PML correction.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        pmlmodule = 'gprMax.pml_updates.pml_updates_electric_' + G.pmlformulation + '_ext'
        func = getattr(import_module(pmlmodule), 'order' + str(len(self.CFS)) + '_' + self.direction)
        func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.EPhi1, self.EPhi2, self.ERA, self.ERB, self.ERE, self.ERF, self.d)

    def update_magnetic(self, G):
        """This functions updates magnetic field components with the PML correction.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        pmlmodule = 'gprMax.pml_updates.pml_updates_magnetic_' + G.pmlformulation + '_ext'
        func = getattr(import_module(pmlmodule), 'order' + str(len(self.CFS)) + '_' + self.direction)
        func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.HPhi1, self.HPhi2, self.HRA, self.HRB, self.HRE, self.HRF, self.d)

    def gpu_set_blocks_per_grid(self, G):
        """Set the blocks per grid size used for updating the PML field arrays on a GPU.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        self.bpg = (int(np.ceil(((self.EPhi1_gpu.shape[1] + 1) * (self.EPhi1_gpu.shape[2] + 1) * (self.EPhi1_gpu.shape[3] + 1)) / G.tpb[0])), 1, 1)

    def gpu_initialise_arrays(self):
        """Initialise PML field and coefficient arrays on GPU."""

        import pycuda.gpuarray as gpuarray

        self.ERA_gpu = gpuarray.to_gpu(self.ERA)
        self.ERB_gpu = gpuarray.to_gpu(self.ERB)
        self.ERE_gpu = gpuarray.to_gpu(self.ERE)
        self.ERF_gpu = gpuarray.to_gpu(self.ERF)
        self.HRA_gpu = gpuarray.to_gpu(self.HRA)
        self.HRB_gpu = gpuarray.to_gpu(self.HRB)
        self.HRE_gpu = gpuarray.to_gpu(self.HRE)
        self.HRF_gpu = gpuarray.to_gpu(self.HRF)

        if self.direction[0] == 'x':
            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=floattype))
            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=floattype))
            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype))
            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=floattype))
        elif self.direction[0] == 'y':
            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=floattype))
            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=floattype))
            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype))
            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=floattype))
        elif self.direction[0] == 'z':
            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=floattype))
            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=floattype))
            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype))
            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype))

    def gpu_get_update_funcs(self, kernelselectric, kernelsmagnetic):
        """Get 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_gpu = kernelselectric.get_function('order' + str(len(self.CFS)) + '_' + self.direction)
        self.update_magnetic_gpu = kernelsmagnetic.get_function('order' + str(len(self.CFS)) + '_' + self.direction)

    def gpu_update_electric(self, G):
        """This functions updates electric field components with the PML correction on the GPU.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        self.update_electric_gpu(np.int32(self.xs), np.int32(self.xf), np.int32(self.ys), np.int32(self.yf), np.int32(self.zs), np.int32(self.zf), np.int32(self.EPhi1_gpu.shape[1]), np.int32(self.EPhi1_gpu.shape[2]), np.int32(self.EPhi1_gpu.shape[3]), np.int32(self.EPhi2_gpu.shape[1]), np.int32(self.EPhi2_gpu.shape[2]), np.int32(self.EPhi2_gpu.shape[3]), np.int32(self.thickness), G.ID_gpu.gpudata, G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata, G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata, self.EPhi1_gpu.gpudata, self.EPhi2_gpu.gpudata, self.ERA_gpu.gpudata, self.ERB_gpu.gpudata, self.ERE_gpu.gpudata, self.ERF_gpu.gpudata, floattype(self.d), block=G.tpb, grid=self.bpg)

    def gpu_update_magnetic(self, G):
        """This functions updates magnetic field components with the PML correction on the GPU.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        self.update_magnetic_gpu(np.int32(self.xs), np.int32(self.xf), np.int32(self.ys), np.int32(self.yf), np.int32(self.zs), np.int32(self.zf), np.int32(self.HPhi1_gpu.shape[1]), np.int32(self.HPhi1_gpu.shape[2]), np.int32(self.HPhi1_gpu.shape[3]), np.int32(self.HPhi2_gpu.shape[1]), np.int32(self.HPhi2_gpu.shape[2]), np.int32(self.HPhi2_gpu.shape[3]), np.int32(self.thickness), G.ID_gpu.gpudata, G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata, G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata, self.HPhi1_gpu.gpudata, self.HPhi2_gpu.gpudata, self.HRA_gpu.gpudata, self.HRB_gpu.gpudata, self.HRE_gpu.gpudata, self.HRF_gpu.gpudata, floattype(self.d), block=G.tpb, grid=self.bpg) def build_pmls(G, pbar):
    """
    This function builds instances of the PML and calculates the initial
        parameters and coefficients including setting profile
        (based on underlying material er and mr from solid array).

    Args:
        G (class): Grid class instance - holds essential parameters describing the model.
        pbar (class): Progress bar class instance.
    """

    for key, value in G.pmlthickness.items():
        if value > 0:
            sumer = 0  # Sum of relative permittivities in PML slab
            summr = 0  # Sum of relative permeabilities in PML slab

            if key[0] == 'x':
                if key == 'x0':
                    pml = PML(G, ID=key, direction='xminus', xf=value, yf=G.ny, zf=G.nz)
                elif key == 'xmax':
                    pml = PML(G, ID=key, direction='xplus', xs=G.nx - value, xf=G.nx, yf=G.ny, zf=G.nz)
                G.pmls.append(pml)
                for j in range(G.ny):
                    for k in range(G.nz):
                        numID = G.solid[pml.xs, j, k]
                        material = next(x for x in G.materials if x.numID == numID)
                        sumer += material.er
                        summr += material.mr
                averageer = sumer / (G.ny * G.nz)
                averagemr = summr / (G.ny * G.nz)

            elif key[0] == 'y':
                if key == 'y0':
                    pml = PML(G, ID=key, direction='yminus', yf=value, xf=G.nx, zf=G.nz)
                elif key == 'ymax':
                    pml = PML(G, ID=key, direction='yplus', ys=G.ny - value, xf=G.nx, yf=G.ny, zf=G.nz)
                G.pmls.append(pml)
                for i in range(G.nx):
                    for k in range(G.nz):
                        numID = G.solid[i, pml.ys, k]
                        material = next(x for x in G.materials if x.numID == numID)
                        sumer += material.er
                        summr += material.mr
                averageer = sumer / (G.nx * G.nz)
                averagemr = summr / (G.nx * G.nz)

            elif key[0] == 'z':
                if key == 'z0':
                    pml = PML(G, ID=key, direction='zminus', zf=value, xf=G.nx, yf=G.ny)
                elif key == 'zmax':
                    pml = PML(G, ID=key, direction='zplus', zs=G.nz - value, xf=G.nx, yf=G.ny, zf=G.nz)
                G.pmls.append(pml)
                for i in range(G.nx):
                    for j in range(G.ny):
                        numID = G.solid[i, j, pml.zs]
                        material = next(x for x in G.materials if x.numID == numID)
                        sumer += material.er
                        summr += material.mr
                averageer = sumer / (G.nx * G.ny)
                averagemr = summr / (G.nx * G.ny)

            pml.calculate_update_coeffs(averageer, averagemr, G)
            pbar.update() def gpu_initialise_rx_arrays(G):
    """Initialise arrays on GPU for receiver coordinates and to store field components for receivers.

    Args:
        G (class): Grid class instance - holds essential parameters describing the model.
    """

    import pycuda.gpuarray as gpuarray

    # Array to store receiver coordinates on GPU
    rxcoords = np.zeros((len(G.rxs), 3), dtype=np.int32)
    for i, rx in enumerate(G.rxs):
        rxcoords[i, 0] = rx.xcoord
        rxcoords[i, 1] = rx.ycoord
        rxcoords[i, 2] = rx.zcoord

    # Array to store field components for receivers on GPU - rows are field components; columns are iterations; pages are receivers
    rxs = np.zeros((len(Rx.gpu_allowableoutputs), G.iterations, len(G.rxs)), dtype=floattype)

    # Copy arrays to GPU
    rxcoords_gpu = gpuarray.to_gpu(rxcoords)
    rxs_gpu = gpuarray.to_gpu(rxs)

    return rxcoords_gpu, rxs_gpu def gpu_get_rx_array(rxs_gpu, rxcoords_gpu, G):
    """Copy output from receivers array used on GPU back to receiver objects.

    Args:
        rxs_gpu (float): numpy array of receiver data from GPU - rows are field components; columns are iterations; pages are receivers.
        rxcoords_gpu (float): numpy array of receiver coordinates from GPU.
        G (class): Grid class instance - holds essential parameters describing the model.
    """

    for rx in G.rxs:
        for rxgpu in range(len(G.rxs)):
            if rx.xcoord == rxcoords_gpu[rxgpu, 0] and rx.ycoord == rxcoords_gpu[rxgpu, 1] and rx.zcoord == rxcoords_gpu[rxgpu, 2]:
                rx.outputs['Ex'] = rxs_gpu[0, :, rxgpu]
                rx.outputs['Ey'] = rxs_gpu[1, :, rxgpu]
                rx.outputs['Ez'] = rxs_gpu[2, :, rxgpu]
                rx.outputs['Hx'] = rxs_gpu[3, :, rxgpu]
                rx.outputs['Hy'] = rxs_gpu[4, :, rxgpu]
                rx.outputs['Hz'] = rxs_gpu[5, :, rxgpu] class Snapshot(object):
    """Snapshots of the electric and magnetic field values."""

    # Dimensions of largest requested snapshot
    nx_max = 0
    ny_max = 0
    nz_max = 0

    # GPU - threads per block
    tpb = (1, 1, 1)
    # GPU - blocks per grid - set according to largest requested snapshot
    bpg = None

    # Set string for byte order
    if sys.byteorder == 'little':
        byteorder = 'LittleEndian'
    else:
        byteorder = 'BigEndian'

    # Set format text and string depending on float type
    if np.dtype(floattype).name == 'float32':
        floatname = 'Float32'
        floatstring = 'f'
    elif np.dtype(floattype).name == 'float64':
        floatname = 'Float64'
        floatstring = 'd'

    def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None, dx=None, dy=None, dz=None, time=None, filename=None):
        """
        Args:
            xs, xf, ys, yf, zs, zf (int): Extent of the volume in cells.
            dx, dy, dz (int): Spatial discretisation in cells.
            time (int): Iteration number to take the snapshot on.
            filename (str): Filename to save to.
        """

        self.xs = xs
        self.ys = ys
        self.zs = zs
        self.xf = xf
        self.yf = yf
        self.zf = zf
        self.dx = dx
        self.dy = dy
        self.dz = dz
        self.nx = round_value((self.xf - self.xs) / self.dx)
        self.ny = round_value((self.yf - self.ys) / self.dy)
        self.nz = round_value((self.zf - self.zs) / self.dz)
        self.sx = slice(self.xs, self.xf + self.dx, self.dx)
        self.sy = slice(self.ys, self.yf + self.dy, self.dy)
        self.sz = slice(self.zs, self.zf + self.dz, self.dz)
        self.ncells = self.nx * self.ny * self.nz
        self.datasizefield = 3 * np.dtype(floattype).itemsize * self.ncells
        self.vtkdatawritesize = 2 * self.datasizefield + 2 * np.dtype(np.uint32).itemsize
        self.time = time
        self.basefilename = filename

    def store(self, G):
        """Store (in memory) electric and magnetic field values for snapshot.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        # Memory views of field arrays to dimensions required for the snapshot
        Exslice = np.ascontiguousarray(G.Ex[self.sx, self.sy, self.sz])
        Eyslice = np.ascontiguousarray(G.Ey[self.sx, self.sy, self.sz])
        Ezslice = np.ascontiguousarray(G.Ez[self.sx, self.sy, self.sz])
        Hxslice = np.ascontiguousarray(G.Hx[self.sx, self.sy, self.sz])
        Hyslice = np.ascontiguousarray(G.Hy[self.sx, self.sy, self.sz])
        Hzslice = np.ascontiguousarray(G.Hz[self.sx, self.sy, self.sz])

        # Create arrays to hold the field data for snapshot
        Exsnap = np.zeros((self.nx, self.ny, self.nz), dtype=floattype)
        Eysnap = np.zeros((self.nx, self.ny, self.nz), dtype=floattype)
        Ezsnap = np.zeros((self.nx, self.ny, self.nz), dtype=floattype)
        Hxsnap = np.zeros((self.nx, self.ny, self.nz), dtype=floattype)
        Hysnap = np.zeros((self.nx, self.ny, self.nz), dtype=floattype)
        Hzsnap = np.zeros((self.nx, self.ny, self.nz), dtype=floattype)

        # Calculate field values at points (comes from averaging field components in cells)
        calculate_snapshot_fields(
            self.nx,
            self.ny,
            self.nz,
            Exslice,
            Eyslice,
            Ezslice,
            Hxslice,
            Hyslice,
            Hzslice,
            Exsnap,
            Eysnap,
            Ezsnap,
            Hxsnap,
            Hysnap,
            Hzsnap)

        # Convert to format for Paraview
        self.electric = np.stack((Exsnap, Eysnap, Ezsnap)).reshape(-1, order='F')
        self.magnetic = np.stack((Hxsnap, Hysnap, Hzsnap)).reshape(-1, order='F')

    def write_vtk_imagedata(self, pbar, G):
        """Write snapshot data to a VTK ImageData (.vti) file.

            N.B. No Python 3 support for VTK at time of writing (03/2015)

        Args:
            pbar (class): Progress bar class instance.
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        hfield_offset = 3 * np.dtype(floattype).itemsize * self.ncells + np.dtype(np.uint32).itemsize

        self.filehandle = open(self.filename, 'wb')
        self.filehandle.write('<?xml version="1.0"?>\n'.encode('utf-8'))
        self.filehandle.write('<VTKFile type="ImageData" version="1.0" byte_order="{}">\n'.format(Snapshot.byteorder).encode('utf-8'))
        self.filehandle.write('<ImageData WholeExtent="{} {} {} {} {} {}" Origin="0 0 0" Spacing="{:.3} {:.3} {:.3}">\n'.format(self.xs, round_value(self.xf / self.dx), self.ys, round_value(self.yf / self.dy), self.zs, round_value(self.zf / self.dz), self.dx * G.dx, self.dy * G.dy, self.dz * G.dz).encode('utf-8'))
        self.filehandle.write('<Piece Extent="{} {} {} {} {} {}">\n'.format(self.xs, round_value(self.xf / self.dx), self.ys, round_value(self.yf / self.dy), self.zs, round_value(self.zf / self.dz)).encode('utf-8'))
        self.filehandle.write('<CellData Vectors="E-field H-field">\n'.encode('utf-8'))
        self.filehandle.write('<DataArray type="{}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'.format(Snapshot.floatname).encode('utf-8'))
        self.filehandle.write('<DataArray type="{}" Name="H-field" NumberOfComponents="3" format="appended" offset="{}" />\n'.format(Snapshot.floatname, hfield_offset).encode('utf-8'))
        self.filehandle.write('</CellData>\n</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))

        # Write number of bytes of appended data as UInt32
        self.filehandle.write(pack('I', self.datasizefield))
        pbar.update(n=4)
        self.electric.tofile(self.filehandle)
        pbar.update(n=self.datasizefield)

        # Write number of bytes of appended data as UInt32
        self.filehandle.write(pack('I', self.datasizefield))
        pbar.update(n=4)
        self.magnetic.tofile(self.filehandle)
        pbar.update(n=self.datasizefield)

        self.filehandle.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))
        self.filehandle.close() def gpu_initialise_snapshot_array(G):
    """Initialise array on GPU for to store field data for snapshots.

    Args:
        G (class): Grid class instance - holds essential parameters describing the model.

    Returns:
        snapE_gpu, snapH_gpu (float): numpy arrays of snapshot data on GPU.
    """

    import pycuda.gpuarray as gpuarray

    # Get dimensions of largest requested snapshot
    for snap in G.snapshots:
        if snap.nx > Snapshot.nx_max:
            Snapshot.nx_max = snap.nx
        if snap.ny > Snapshot.ny_max:
            Snapshot.ny_max = snap.ny
        if snap.nz > Snapshot.nz_max:
            Snapshot.nz_max = snap.nz

    # GPU - blocks per grid - according to largest requested snapshot
    Snapshot.bpg = (int(np.ceil(((Snapshot.nx_max) * (Snapshot.ny_max) * (Snapshot.nz_max)) / Snapshot.tpb[0])), 1, 1)

    # 4D arrays to store snapshots on GPU, e.g. snapEx(time, x, y, z)
    numsnaps = 1 if G.snapsgpu2cpu else len(G.snapshots)
    snapEx = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)
    snapEy = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)
    snapEz = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)
    snapHx = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)
    snapHy = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)
    snapHz = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)

    # Copy arrays to GPU
    snapEx_gpu = gpuarray.to_gpu(snapEx)
    snapEy_gpu = gpuarray.to_gpu(snapEy)
    snapEz_gpu = gpuarray.to_gpu(snapEz)
    snapHx_gpu = gpuarray.to_gpu(snapHx)
    snapHy_gpu = gpuarray.to_gpu(snapHy)
    snapHz_gpu = gpuarray.to_gpu(snapHz)

    return snapEx_gpu, snapEy_gpu, snapEz_gpu, snapHx_gpu, snapHy_gpu, snapHz_gpu def gpu_get_snapshot_array(snapEx_gpu, snapEy_gpu, snapEz_gpu, snapHx_gpu, snapHy_gpu, snapHz_gpu, i, snap):
    """Copy snapshot array used on GPU back to snapshot objects and store in format for Paraview.

    Args:
        snapE_gpu, snapH_gpu (float): numpy arrays of snapshot data from GPU.
        i (int): index for snapshot data on GPU array.
        snap (class): Snapshot class instance
    """

    snap.electric = np.stack((snapEx_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
                              snapEy_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
                              snapEz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F')
    snap.magnetic = np.stack((snapHx_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
                              snapHy_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
                              snapHz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F') def gpu_initialise_src_arrays(sources, G):
    """Initialise arrays on GPU for source coordinates/polarisation, 
        other source information, and source waveform values.

    Args:
        sources (list): List of sources of one class, e.g. HertzianDipoles.
        G (class): Grid class instance - holds essential parameters describing the model.

    Returns:
        srcinfo1_gpu (int): numpy array of source cell coordinates and polarisation information.
        srcinfo2_gpu (float): numpy array of other source information, e.g. length, resistance etc...
        srcwaves_gpu (float): numpy array of source waveform values.
    """

    import pycuda.gpuarray as gpuarray

    srcinfo1 = np.zeros((len(sources), 4), dtype=np.int32)
    srcinfo2 = np.zeros((len(sources)), dtype=floattype)
    srcwaves = np.zeros((len(sources), G.iterations), dtype=floattype)
    for i, src in enumerate(sources):
        srcinfo1[i, 0] = src.xcoord
        srcinfo1[i, 1] = src.ycoord
        srcinfo1[i, 2] = src.zcoord

        if src.polarisation == 'x':
            srcinfo1[i, 3] = 0
        elif src.polarisation == 'y':
            srcinfo1[i, 3] = 1
        elif src.polarisation == 'z':
            srcinfo1[i, 3] = 2

        if src.__class__.__name__ == 'HertzianDipole':
            srcinfo2[i] = src.dl
            srcwaves[i, :] = src.waveformvalues_wholestep
        elif src.__class__.__name__ == 'VoltageSource':
            if src.resistance:
                srcinfo2[i] = src.resistance
                srcwaves[i, :] = src.waveformvalues_wholestep
            else:
                srcinfo2[i] = 0
                srcwaves[i, :] = src.waveformvalues_halfstep
        elif src.__class__.__name__ == 'MagneticDipole':
            srcwaves[i, :] = src.waveformvalues_halfstep

    srcinfo1_gpu = gpuarray.to_gpu(srcinfo1)
    srcinfo2_gpu = gpuarray.to_gpu(srcinfo2)
    srcwaves_gpu = gpuarray.to_gpu(srcwaves)

    return srcinfo1_gpu, srcinfo2_gpu, srcwaves_gpu def get_host_info():
    """Get information about the machine, CPU, RAM, and OS.

    Returns:
        hostinfo (dict): Manufacturer and model of machine; description of CPU
                type, speed, cores; RAM; name and version of operating system.
    """

    # Default to 'unknown' if any of the detection fails
    manufacturer = model = cpuID = sockets = threadspercore = 'unknown'

    # Windows
    if sys.platform == 'win32':
        # Manufacturer/model
        try:
            manufacturer = subprocess.check_output("wmic csproduct get vendor", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            manufacturer = manufacturer.split('\n')
            if len(manufacturer) > 1:
                manufacturer = manufacturer[1]
            else:
                manufacturer = manufacturer[0]
            model = subprocess.check_output("wmic computersystem get model", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            model = model.split('\n')
            if len(model) > 1:
                model = model[1]
            else:
                model = model[0]
        except subprocess.CalledProcessError:
            pass
        machineID = manufacturer + ' ' + model

        # CPU information
        try:
            allcpuinfo = subprocess.check_output("wmic cpu get Name", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            allcpuinfo = allcpuinfo.split('\n')
            sockets = 0
            for line in allcpuinfo:
                if 'CPU' in line:
                    cpuID = line.strip()
                    sockets += 1
        except subprocess.CalledProcessError:
            pass

        physicalcores = psutil.cpu_count(logical=False)
        logicalcores = psutil.cpu_count(logical=True)

        # OS version
        if platform.machine().endswith('64'):
            osbit = ' (64-bit)'
        else:
            osbit = ' (32-bit)'
        osversion = 'Windows ' + platform.release() + osbit

    # Mac OS X/macOS
    elif sys.platform == 'darwin':
        # Manufacturer/model
        manufacturer = 'Apple'
        try:
            model = subprocess.check_output("sysctl -n hw.model", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
        except subprocess.CalledProcessError:
            pass
        machineID = manufacturer + ' ' + model

        # CPU information
        try:
            sockets = subprocess.check_output("sysctl -n hw.packages", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            sockets = int(sockets)
            cpuID = subprocess.check_output("sysctl -n machdep.cpu.brand_string", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            cpuID = ' '.join(cpuID.split())
            physicalcores = subprocess.check_output("sysctl -n hw.physicalcpu", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            physicalcores = int(physicalcores)
            logicalcores = subprocess.check_output("sysctl -n hw.logicalcpu", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            logicalcores = int(logicalcores)
        except subprocess.CalledProcessError:
            pass

        # OS version
        if int(platform.mac_ver()[0].split('.')[1]) < 12:
            osversion = 'Mac OS X (' + platform.mac_ver()[0] + ')'
        else:
            osversion = 'macOS (' + platform.mac_ver()[0] + ')'

    # Linux
    elif sys.platform == 'linux':
        # Manufacturer/model
        try:
            manufacturer = subprocess.check_output("cat /sys/class/dmi/id/sys_vendor", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
            model = subprocess.check_output("cat /sys/class/dmi/id/product_name", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
        except subprocess.CalledProcessError:
            pass
        machineID = manufacturer + ' ' + model

        # CPU information
        try:
            my_env = os.environ.copy()
            my_env["LC_ALL"] = "C"
            cpuIDinfo = subprocess.check_output("cat /proc/cpuinfo", shell=True, env=my_env, stderr=subprocess.STDOUT).decode('utf-8').strip()
            for line in cpuIDinfo.split('\n'):
                if re.search('model name', line):
                    cpuID = re.sub('.*model name.*:', '', line, 1).strip()
            allcpuinfo = subprocess.check_output("lscpu", shell=True, env=my_env, stderr=subprocess.STDOUT).decode('utf-8').strip()
            for line in allcpuinfo.split('\n'):
                if 'Socket(s)' in line:
                    sockets = int(re.sub("\D", "", line.strip()))
                if 'Thread(s) per core' in line:
                    threadspercore = int(re.sub("\D", "", line.strip()))
                if 'Core(s) per socket' in line:
                    corespersocket = int(re.sub("\D", "", line.strip()))
        except subprocess.CalledProcessError:
            pass

        physicalcores = sockets * corespersocket
        logicalcores = sockets * corespersocket * threadspercore

        # OS version
        osversion = platform.platform()

    # Dictionary of host information
    hostinfo = {}
    hostinfo['hostname'] = platform.node()
    hostinfo['machineID'] = machineID.strip()
    hostinfo['sockets'] = sockets
    hostinfo['cpuID'] = cpuID
    hostinfo['osversion'] = osversion

    # Hyperthreading
    if logicalcores != physicalcores:
        hostinfo['hyperthreading'] = True
    else:
        hostinfo['hyperthreading'] = False

    hostinfo['logicalcores'] = logicalcores
    # Number of physical CPU cores, i.e. avoid hyperthreading with OpenMP
    hostinfo['physicalcores'] = physicalcores

    # Handle case where cpu_count returns None on some machines
    if not hostinfo['physicalcores']:
        hostinfo['physicalcores'] = hostinfo['logicalcores']
    hostinfo['ram'] = psutil.virtual_memory().total

    return hostinfo def get_terminal_width():
    """Get/set width of terminal being used.

    Returns:
        terminalwidth (int): Terminal width
    """

    terminalwidth = get_terminal_size()[0]
    if terminalwidth == 0:
        terminalwidth = 100

    return terminalwidth def human_size(size, a_kilobyte_is_1024_bytes=False):
    """Convert a file size to human-readable form.

    Args:
        size (int): file size in bytes.
        a_kilobyte_is_1024_bytes (boolean) - true for multiples of 1024, false for multiples of 1000.

    Returns:
        Human-readable (string).
    """

    suffixes = {1000: ['KB', 'MB', 'GB', 'TB', 'PB', 'EB', 'ZB', 'YB'], 1024: ['KiB', 'MiB', 'GiB', 'TiB', 'PiB', 'EiB', 'ZiB', 'YiB']}

    if size < 0:
        raise ValueError('Number must be non-negative.')

    multiple = 1024 if a_kilobyte_is_1024_bytes else 1000
    for suffix in suffixes[multiple]:
        size /= multiple
        if size < multiple:
            return '{:.3g}{}'.format(size, suffix)

    raise ValueError('Number is too large.') @contextmanager
def open_path_file(path_or_file):
    """
    Accepts either a path as a string or a file object and returns a file
    object (http://stackoverflow.com/a/6783680).

    Args:
        path_or_file: path as a string or a file object.

    Returns:
        f (object): File object.
    """

    if isinstance(path_or_file, str):
        f = file_to_close = codecs.open(path_or_file, 'r', encoding='utf-8')
    else:
        f = path_or_file
        file_to_close = None

    try:
        yield f
    finally:
        if file_to_close:
            file_to_close.close() def round32(value):
    """Rounds up to nearest multiple of 32."""
    return int(32 * np.ceil(float(value) / 32)) def timer():
    """Function to return the current process wide time in fractional seconds."""
    return perf_counter()
import datetime
import itertools
import os
import psutil
import sys
import numpy as np
    import pycuda.driver as drv
from importlib import import_module
from colorama import init
from colorama import Fore
from colorama import Style
from terminaltables import AsciiTable
from tqdm import tqdm
from gprMax.constants import floattype
from gprMax.constants import complextype
from gprMax.constants import cudafloattype
from gprMax.constants import cudacomplextype
from gprMax.exceptions import GeneralError
from gprMax.fields_outputs import store_outputs
from gprMax.fields_outputs import kernel_template_store_outputs
from gprMax.fields_outputs import write_hdf5_outputfile
from gprMax.fields_updates_ext import update_electric
from gprMax.fields_updates_ext import update_magnetic
from gprMax.fields_updates_ext import update_electric_dispersive_multipole_A
from gprMax.fields_updates_ext import update_electric_dispersive_multipole_B
from gprMax.fields_updates_ext import update_electric_dispersive_1pole_A
from gprMax.fields_updates_ext import update_electric_dispersive_1pole_B
from gprMax.fields_updates_gpu import kernels_template_fields
from gprMax.grid import FDTDGrid
from gprMax.grid import dispersion_analysis
from gprMax.input_cmds_geometry import process_geometrycmds
from gprMax.input_cmds_file import process_python_include_code
from gprMax.input_cmds_file import write_processed_file
from gprMax.input_cmds_file import check_cmd_names
from gprMax.input_cmds_multiuse import process_multicmds
from gprMax.input_cmds_singleuse import process_singlecmds
from gprMax.materials import Material
from gprMax.materials import process_materials
from gprMax.pml import CFS
from gprMax.pml import PML
from gprMax.pml import build_pmls
from gprMax.receivers import gpu_initialise_rx_arrays
from gprMax.receivers import gpu_get_rx_array
from gprMax.snapshots import Snapshot
from gprMax.snapshots import gpu_initialise_snapshot_array
from gprMax.snapshots import gpu_get_snapshot_array
from gprMax.snapshots_gpu import kernel_template_store_snapshot
from gprMax.sources import gpu_initialise_src_arrays
from gprMax.source_updates_gpu import kernels_template_sources
from gprMax.utilities import get_host_info
from gprMax.utilities import get_terminal_width
from gprMax.utilities import human_size
from gprMax.utilities import open_path_file
from gprMax.utilities import round32
from gprMax.utilities import timer
from gprMax.yee_cell_build_ext import build_electric_components
from gprMax.yee_cell_build_ext import build_magnetic_components
    from pycuda.compiler import SourceModule
    G.gpu_initialise_arrays()

    # PML updates
    if G.pmls:
        # Prepare kernels
        pmlmodulelectric = 'gprMax.pml_updates.pml_updates_electric_' + G.pmlformulation + '_gpu'
        kernelelectricfunc = getattr(import_module(pmlmodulelectric), 'kernels_template_pml_electric_' + G.pmlformulation)
        pmlmodulemagnetic = 'gprMax.pml_updates.pml_updates_magnetic_' + G.pmlformulation + '_gpu'
        kernelmagneticfunc = getattr(import_module(pmlmodulemagnetic), 'kernels_template_pml_magnetic_' + G.pmlformulation)
        kernels_pml_electric = SourceModule(kernelelectricfunc.substitute(REAL=cudafloattype, N_updatecoeffsE=G.updatecoeffsE.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NX_FIELDS=G.nx + 1, NY_FIELDS=G.ny + 1, NZ_FIELDS=G.nz + 1, NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3]), options=compiler_opts)
        kernels_pml_magnetic = SourceModule(kernelmagneticfunc.substitute(REAL=cudafloattype, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsH.shape[1], NX_FIELDS=G.nx + 1, NY_FIELDS=G.ny + 1, NZ_FIELDS=G.nz + 1, NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3]), options=compiler_opts)
        # Copy material coefficient arrays to constant memory of GPU (must be <64KB) for PML kernels
        updatecoeffsE = kernels_pml_electric.get_global('updatecoeffsE')[0]
        updatecoeffsH = kernels_pml_magnetic.get_global('updatecoeffsH')[0]
        drv.memcpy_htod(updatecoeffsE, G.updatecoeffsE)
        drv.memcpy_htod(updatecoeffsH, G.updatecoeffsH)
        # Set block per grid, initialise arrays on GPU, and get kernel functions
        for pml in G.pmls:
            pml.gpu_initialise_arrays()
            pml.gpu_get_update_funcs(kernels_pml_electric, kernels_pml_magnetic)
            pml.gpu_set_blocks_per_grid(G)

    # Receivers
    if G.rxs:
        # Initialise arrays on GPU
        rxcoords_gpu, rxs_gpu = gpu_initialise_rx_arrays(G)
        # Prepare kernel and get kernel function
        kernel_store_outputs = SourceModule(kernel_template_store_outputs.substitute(REAL=cudafloattype, NY_RXCOORDS=3, NX_RXS=6, NY_RXS=G.iterations, NZ_RXS=len(G.rxs), NX_FIELDS=G.nx + 1, NY_FIELDS=G.ny + 1, NZ_FIELDS=G.nz + 1), options=compiler_opts)
        store_outputs_gpu = kernel_store_outputs.get_function("store_outputs")

    # Sources - initialise arrays on GPU, prepare kernel and get kernel functions
    if G.voltagesources + G.hertziandipoles + G.magneticdipoles:
        kernels_sources = SourceModule(kernels_template_sources.substitute(REAL=cudafloattype, N_updatecoeffsE=G.updatecoeffsE.size, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NY_SRCINFO=4, NY_SRCWAVES=G.iterations, NX_FIELDS=G.nx + 1, NY_FIELDS=G.ny + 1, NZ_FIELDS=G.nz + 1, NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3]), options=compiler_opts)
        # Copy material coefficient arrays to constant memory of GPU (must be <64KB) for source kernels
        updatecoeffsE = kernels_sources.get_global('updatecoeffsE')[0]
        updatecoeffsH = kernels_sources.get_global('updatecoeffsH')[0]
        drv.memcpy_htod(updatecoeffsE, G.updatecoeffsE)
        drv.memcpy_htod(updatecoeffsH, G.updatecoeffsH)
        if G.hertziandipoles:
            srcinfo1_hertzian_gpu, srcinfo2_hertzian_gpu, srcwaves_hertzian_gpu = gpu_initialise_src_arrays(G.hertziandipoles, G)
            update_hertzian_dipole_gpu = kernels_sources.get_function("update_hertzian_dipole")
        if G.magneticdipoles:
            srcinfo1_magnetic_gpu, srcinfo2_magnetic_gpu, srcwaves_magnetic_gpu = gpu_initialise_src_arrays(G.magneticdipoles, G)
            update_magnetic_dipole_gpu = kernels_sources.get_function("update_magnetic_dipole")
        if G.voltagesources:
            srcinfo1_voltage_gpu, srcinfo2_voltage_gpu, srcwaves_voltage_gpu = gpu_initialise_src_arrays(G.voltagesources, G)
            update_voltage_source_gpu = kernels_sources.get_function("update_voltage_source")

    # Snapshots - initialise arrays on GPU, prepare kernel and get kernel functions
    if G.snapshots:
        # Initialise arrays on GPU
        snapEx_gpu, snapEy_gpu, snapEz_gpu, snapHx_gpu, snapHy_gpu, snapHz_gpu = gpu_initialise_snapshot_array(G)
        # Prepare kernel and get kernel function
        kernel_store_snapshot = SourceModule(kernel_template_store_snapshot.substitute(REAL=cudafloattype, NX_SNAPS=Snapshot.nx_max, NY_SNAPS=Snapshot.ny_max, NZ_SNAPS=Snapshot.nz_max, NX_FIELDS=G.nx + 1, NY_FIELDS=G.ny + 1, NZ_FIELDS=G.nz + 1), options=compiler_opts)
        store_snapshot_gpu = kernel_store_snapshot.get_function("store_snapshot")

    # Iteration loop timer
    iterstart = drv.Event()
    iterend = drv.Event()
    iterstart.record()

    for iteration in tqdm(range(G.iterations), desc='Running simulation, model ' + str(currentmodelrun) + '/' + str(modelend), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not G.progressbars):

        # Get GPU memory usage on final iteration
        if iteration == G.iterations - 1:
            memsolve = drv.mem_get_info()[1] - drv.mem_get_info()[0]

        # Store field component values for every receiver
        if G.rxs:
            store_outputs_gpu(np.int32(len(G.rxs)), np.int32(iteration),
                              rxcoords_gpu.gpudata, rxs_gpu.gpudata,
                              G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata,
                              G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,
                              block=(1, 1, 1), grid=(round32(len(G.rxs)), 1, 1))

        # Store any snapshots
        for i, snap in enumerate(G.snapshots):
            if snap.time == iteration + 1:
                if not G.snapsgpu2cpu:
                    store_snapshot_gpu(np.int32(i), np.int32(snap.xs),
                                       np.int32(snap.xf), np.int32(snap.ys),
                                       np.int32(snap.yf), np.int32(snap.zs),
                                       np.int32(snap.zf), np.int32(snap.dx),
                                       np.int32(snap.dy), np.int32(snap.dz),
                                       G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata,
                                       G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,
                                       snapEx_gpu.gpudata, snapEy_gpu.gpudata, snapEz_gpu.gpudata,
                                       snapHx_gpu.gpudata, snapHy_gpu.gpudata, snapHz_gpu.gpudata,
                                       block=Snapshot.tpb, grid=Snapshot.bpg)
                else:
                    store_snapshot_gpu(np.int32(0), np.int32(snap.xs),
                                       np.int32(snap.xf), np.int32(snap.ys),
                                       np.int32(snap.yf), np.int32(snap.zs),
                                       np.int32(snap.zf), np.int32(snap.dx),
                                       np.int32(snap.dy), np.int32(snap.dz),
                                       G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata,
                                       G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,
                                       snapEx_gpu.gpudata, snapEy_gpu.gpudata, snapEz_gpu.gpudata,
                                       snapHx_gpu.gpudata, snapHy_gpu.gpudata, snapHz_gpu.gpudata,
                                       block=Snapshot.tpb, grid=Snapshot.bpg)
                    gpu_get_snapshot_array(snapEx_gpu.get(), snapEy_gpu.get(), snapEz_gpu.get(),
                                           snapHx_gpu.get(), snapHy_gpu.get(), snapHz_gpu.get(), 0, snap)

        # Update magnetic field components
        update_h_gpu(np.int32(G.nx), np.int32(G.ny), np.int32(G.nz),
                     G.ID_gpu.gpudata, G.Hx_gpu.gpudata, G.Hy_gpu.gpudata,
                     G.Hz_gpu.gpudata, G.Ex_gpu.gpudata, G.Ey_gpu.gpudata,
                     G.Ez_gpu.gpudata, block=G.tpb, grid=G.bpg)

        # Update magnetic field components with the PML correction
        for pml in G.pmls:
            pml.gpu_update_magnetic(G)

        # Update magnetic field components for magetic dipole sources
        if G.magneticdipoles:
            update_magnetic_dipole_gpu(np.int32(len(G.magneticdipoles)), np.int32(iteration),
                                       floattype(G.dx), floattype(G.dy), floattype(G.dz),
                                       srcinfo1_magnetic_gpu.gpudata, srcinfo2_magnetic_gpu.gpudata,
                                       srcwaves_magnetic_gpu.gpudata, G.ID_gpu.gpudata,
                                       G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,

                                       block=(1, 1, 1), grid=(round32(len(G.magneticdipoles)), 1, 1))