diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index e1aa521c..17987f1d 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -68,7 +68,7 @@ def main(): def api(inputfile, n=1, mpi=False, benchmark=False, geometry_only=False, geometry_fixed=False, write_processed=False, opt_taguchi=False): """If installed as a module this is the entry point.""" - + class ImportArguments: pass @@ -88,11 +88,11 @@ def api(inputfile, n=1, mpi=False, benchmark=False, geometry_only=False, geometr def run_main(args): """Top-level function that controls what mode of simulation (standard/optimsation/benchmark etc...) is run. - + Args: args (dict): Namespace with input arguments from command line or api. """ - + numbermodelruns = args.n inputdirectory = os.path.dirname(os.path.abspath(args.inputfile)) inputfile = os.path.abspath(os.path.join(inputdirectory, os.path.basename(args.inputfile))) diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index ebedad37..e751536a 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/input_cmds_file.py @@ -25,32 +25,32 @@ from gprMax.utilities import ListStream 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 commands and insert the contents of the included file at that location. - + Args: inputfile (str): Name of the input file to open. 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. """ - + userpython = False - + with open(inputfile, 'r') as f: # Strip out any newline characters and comments that must begin with double hashes inputlines = [line.rstrip() for line in f if(not line.startswith('##') and line.rstrip('\n'))] - + # List to hold final processed commands processedlines = [] - + x = 0 while(x < len(inputlines)): - + # Process any Python code if(inputlines[x].startswith('#python:')): # Save stdout location to restore later stdout = sys.stdout - + # String to hold Python code to be executed pythoncode = '' x += 1 @@ -66,36 +66,36 @@ def process_python_include_code(inputfile, usernamespace): sys.stdout = codeout = ListStream() # Execute code block & make available only usernamespace exec(pythoncompiledcode, usernamespace) - + # Now strip out any lines that don't begin with a hash command codeproc = [line + ('\n') for line in codeout.data if(line.startswith('#'))] - + # Add processed Python code to list processedlines.extend(codeproc) - + # Reset stdio sys.stdout = stdout - + # Process any include commands elif(inputlines[x].startswith('#include_file:')): includefile = inputlines[x].split() - + if len(includefile) != 2: raise CmdInputError('#include_file requires exactly one parameter') - + includefile = includefile[1] - + # See if file exists at specified path and if not try input file directory if not os.path.isfile(includefile): includefile = os.path.join(usernamespace['inputdirectory'], includefile) - + with open(includefile, 'r') as f: # Strip out any newline characters and comments that must begin with double hashes includelines = [includeline.rstrip() + '\n' for includeline in f if(not includeline.startswith('##') and includeline.rstrip('\n'))] - + # Add lines from include file to list processedlines.extend(includelines) - + # Add any other commands to list elif(inputlines[x].startswith('#')): # Add gprMax command to list @@ -103,20 +103,20 @@ def process_python_include_code(inputfile, usernamespace): processedlines.append(inputlines[x]) x += 1 - + return processedlines def write_processed_file(inputfile, modelrun, numbermodelruns, processedlines): """Writes an input file after any Python code and include commands in the original input file have been processed. - + Args: inputfile (str): Name of the input file to open. modelrun (int): Current model run number. numbermodelruns (int): Total number of model runs. processedlines (list): Input commands after after processing any Python code and include commands. """ - + if numbermodelruns == 1: processedfile = os.path.splitext(inputfile)[0] + '_processed.in' else: @@ -131,11 +131,11 @@ def write_processed_file(inputfile, modelrun, numbermodelruns, processedlines): 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. @@ -150,13 +150,13 @@ def check_cmd_names(processedlines, checkessential=True): singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#taguchi', '#end_taguchi'], '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', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#transmission_line', '#tem_transmission_line', '#rx', '#rx_box', '#snapshot', '#pml_cfs']} - + multiplecmds = {key: [] for key in ['#geometry_view', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#transmission_line', '#rx', '#rx_box', '#snapshot', '#pml_cfs']} + # Geometry object building commands that there can be multiple instances of in a model - these will be lists within the dictionary geometrycmds = ['#geometry_objects_file', '#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 @@ -171,22 +171,22 @@ def check_cmd_names(processedlines, checkessential=True): # Count essential commands if cmdname in essentialcmds: countessentialcmds += 1 - + # Assign command parameters as values to dictionary keys if cmdname in singlecmds: if singlecmds[cmdname] == 'None': singlecmds[cmdname] = cmd[1].strip(' \t\n') else: raise CmdInputError('You can only have one ' + cmdname + ' commmand 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 gprMax commands required to run a model. Essential commands are: ' + ', '.join(essentialcmds)) diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 33a6a4ee..b69c3b5d 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -1,4 +1,4 @@ - # Copyright (C) 2015-2016: The University of Edinburgh +# Copyright (C) 2015-2016: The University of Edinburgh # Authors: Craig Warren and Antonis Giannopoulos # # This file is part of gprMax. @@ -16,765 +16,708 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from .constants import z0 -from .exceptions import CmdInputError -from .geometry_views import GeometryView -from .materials import Material, PeplinskiSoil -from .pml import CFSParameter, CFS -from .receivers import Rx -from .snapshots import Snapshot -from .sources import ( - VoltageSource, HertzianDipole, MagneticDipole, TransmissionLine, TEMTransmissionLine) -from .utilities import round_value -from .waveforms import Waveform - -from collections import OrderedDict - - -def check_coordinates(x, y, z, G, cmdname, tmp, name=''): - # Check if coordinates are within the bounds of the grid - 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) - - -def create_cmd_msg(cmdname, cmdinstance): - """Function to create a formatted str relating to a command - and its parameters - """ - return "'{}:{}'".format(cmdname, cmdinstance) - - -def create_transmission_line(cmdname, cmdinstance, G): - - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - tmp = cmdinstance.split() - - # Check correct number of args - if len(tmp) < 6: - raise CmdInputError( - "{} requires at least six parameters".format(cmd_msg)) - - # Check polarity & position parameters - if tmp[0].lower() not in ('x', 'y', 'z', 'zx', 'xz'): - raise CmdInputError( - "{} polarisation must be x, y, z or zx".format(cmd_msg)) - - 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, G, cmdname, tmp) - - if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: - print("WARNING: {} sources and receivers should not normally be \ - positioned within the PML.".format(cmd_msg)) - if resistance <= 0 or resistance > z0: - raise CmdInputError( - "{} requires a resistance greater than zero and less than the impedance \ - of free space, i.e. 376.73 Ohms".format(cmd_msg)) - - # Check if there is a waveformID in the waveforms list - if not any(x.ID == tmp[5] for x in G.waveforms): - raise CmdInputError("{} there is no waveform with the identifier \ - {}".format(cmd_msg, tmp[4])) - - # Can create several types of transmission line - if cmdname == '#transmission_line': - t = TransmissionLine(G) - elif cmdname == '#tem_transmission_line': - t = TEMTransmissionLine(G) - else: - raise CmdInputError("{} is not a valid transmission line type") - - t.polarisation = tmp[0] - t.xcoord = xcoord - t.ycoord = ycoord - t.zcoord = zcoord - t.setID() - t.resistance = resistance - t.waveformID = tmp[5] - t.calculate_incident_V_I(G) - - if len(tmp) > 6: - # Check source start & source remove time parameters - start = float(tmp[6]) - stop = float(tmp[7]) - if start < 0: - raise CmdInputError("{} delay of the initiation of the source should \ - not be less than zero".format(cmd_msg)) - if stop < 0: - raise CmdInputError("{} time to remove the source should not be less \ - than zero".format(cmd_msg)) - if stop - start <= 0: - raise CmdInputError("{} duration of the source should not be zero or \ - less".format(cmd_msg)) - 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 = ' ' - - 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) - - -def create_voltage_source(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) < 6: - raise CmdInputError("{} requires at least six parameters".format( - cmd_msg)) - - # Check polarity & position parameters - if tmp[0].lower() not in ('x', 'y', 'z'): - raise CmdInputError("{} polarisation must be x, y, or z".format( - cmd_msg)) - 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, G, cmdname, tmp) - if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: - print("WARNING: {} sources and receivers should not normally be \ - positioned within the PML.".format(cmd_msg)) - if resistance < 0: - raise CmdInputError("{} requires a source resistance of zero or \ - greater".format(cmd_msg)) - - # Check if there is a waveformID in the waveforms list - if not any(x.ID == tmp[5] for x in G.waveforms): - raise CmdInputError("{} there is no waveform with the \ - identifier {}".format(cmd_msg, tmp[5])) - - v = VoltageSource() - v.polarisation = tmp[0] - v.xcoord = xcoord - v.ycoord = ycoord - v.zcoord = zcoord - v.ID = 'VoltageSource(' + 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("{} delay of the initiation of the source should \ - not be less than zero".format(cmd_msg)) - if stop < 0: - raise CmdInputError("{} time to remove the source should not be less \ - than zero".format(cmd_msg)) - if stop - start <= 0: - raise CmdInputError("{} duration of the source should not be zero or \ - less".format(cmd_msg)) - 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 = ' ' - - 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) - - -def create_waveforms(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) != 4: - raise CmdInputError("{} requires exactly four parameters".format( - cmd_msg)) - if tmp[0].lower() not in Waveform.types: - raise CmdInputError("{} must have one of the following types \ - {}".format(cmd_msg, ','.join(Waveform.types))) - if float(tmp[2]) <= 0: - raise CmdInputError("{} requires an excitation frequency value of greater \ - than zero".format(cmd_msg)) - if any(x.ID == tmp[3] for x in G.waveforms): - raise CmdInputError("{} with ID {} already exists".format( - cmd_msg, 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 amplitude {:g}, frequency {:g}Hz created.'.format(w.ID, w.type, w.amp, w.freq)) - - G.waveforms.append(w) - - -def create_hertzian_dipole(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) < 5: - raise CmdInputError("{} requires at least five parameters".format(cmd_msg)) - - # Check polarity & position parameters - if tmp[0].lower() not in ('x', 'y', 'z'): - raise CmdInputError("{} polarisation must be x, y, or z".format(cmd_msg)) - 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, G, cmdname, tmp) - if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: - print("WARNING: {} sources and receivers should not normally be positioned within the PML.".format(cmd_msg)) - - # Check if there is a waveformID in the waveforms list - if not any(x.ID == tmp[4] for x in G.waveforms): - raise CmdInputError("{} there is no waveform with the identifier {}".format(cmd_msg, tmp[4])) - - h = HertzianDipole() - h.polarisation = tmp[0] - h.xcoord = xcoord - h.ycoord = ycoord - h.zcoord = zcoord - h.xcoordbase = xcoord - h.ycoordbase = ycoord - h.zcoordbase = zcoord - h.ID = 'HertzianDipole(' + 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("{} delay of the initiation of the source should not be less than zero".format(cmd_msg)) - if stop < 0: - raise CmdInputError("{} time to remove the source should not be less than zero".format(cmd_msg)) - if stop - start <= 0: - raise CmdInputError("{} duration of the source should not be zero or less".format(cmd_msg)) - 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 = ' ' - - if G.messages: - 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) - - -def create_magnetic_dipole(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) < 5: - raise CmdInputError("{} requires at least five parameters".format(cmd_msg)) - - # Check polarity & position parameters - if tmp[0].lower() not in ('x', 'y', 'z'): - raise CmdInputError("{} polarisation must be x, y, or z".format(cmd_msg)) - 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, G, cmdname, tmp) - if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: - print("WARNING: {} sources and receivers should not normally be positioned within the PML.".format(cmd_msg)) - - # Check if there is a waveformID in the waveforms list - if not any(x.ID == tmp[4] for x in G.waveforms): - raise CmdInputError("{} there is no waveform with the identifier {}".format(cmd_msg, tmp[4])) - - m = MagneticDipole() - m.polarisation = tmp[0] - m.xcoord = xcoord - m.ycoord = ycoord - m.zcoord = zcoord - m.xcoordbase = xcoord - m.ycoordbase = ycoord - m.zcoordbase = zcoord - m.ID = 'MagneticDipole(' + 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("{} delay of the initiation of the source should not be less than zero".format(cmd_msg)) - if stop < 0: - raise CmdInputError("{} time to remove the source should not be less than zero".format(cmd_msg)) - if stop - start <= 0: - raise CmdInputError("{} duration of the source should not be zero or less".format(cmd_msg)) - 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 = ' ' - - 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) - - -def create_rx(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) != 3 and len(tmp) < 5: - raise CmdInputError("{} has an incorrect number of parameters".format(cmd_msg)) - - # 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[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: - print("WARNING: {} sources and receivers should not normally be positioned within the PML.\n".format(cmd_msg)) - - r = Rx() - r.xcoord = xcoord - r.ycoord = ycoord - r.zcoord = zcoord - r.xcoordbase = xcoord - r.ycoordbase = ycoord - r.zcoordbase = zcoord - - # If no ID or outputs are specified, use default i.e Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy, Iz - if len(tmp) == 3: - r.ID = 'Rx(' + str(r.xcoord) + ',' + str(r.ycoord) + ',' + str(r.zcoord) + ')' - r.outputs = Rx.availableoutputs[0:9] - else: - r.ID = tmp[3] - # Check and add field output names - for field in tmp[4::]: - if field in Rx.availableoutputs: - r.outputs.append(field) - else: - raise CmdInputError("{} contains an output type that is not available".format(cmd_msg)) - - if G.messages: - print('Receiver at {:g}m, {:g}m, {:g}m with output(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs))) - - G.rxs.append(r) - - -def create_rx_box(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) != 9: - raise CmdInputError("{} requires exactly nine parameters".format(cmd_msg)) - - 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, G, cmdname, tmp, name='lower') - check_coordinates(xf, yf, zf, G, cmdname, tmp, name='upper') - - if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: - print("WARNING: {} sources and receivers should not normally be positioned within the PML.".format(cmd_msg)) - if xs >= xf or ys >= yf or zs >= zf: - raise CmdInputError("{} the lower coordinates should be less than the upper coordinates".format(cmd_msg)) - if dx < 0 or dy < 0 or dz < 0: - raise CmdInputError("{} the step size should not be less than zero".format(cmd_msg)) - if dx < G.dx or dy < G.dy or dz < G.dz: - raise CmdInputError("{} the step size should not be less than the spatial discretisation".format(cmd_msg)) - - for x in range(xs, xf, dx): - for y in range(ys, yf, dy): - for z in range(zs, zf, dz): - r = Rx() - r.xcoord = x - r.ycoord = y - r.zcoord = z - r.xcoordbase = x - r.ycoordbase = y - r.zcoordbase = z - r.ID = 'Rx(' + str(x) + ',' + str(y) + ',' + str(z) + ')' - G.rxs.append(r) - - if G.messages: - print('Receiver box {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g} 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)) - - -def create_snapshot(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) != 11: - raise CmdInputError("{} requires exactly eleven parameters".format(cmd_msg)) - - 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: - time = float(tmp[9]) - if time > 0: - time = round_value((time / G.dt)) + 1 - else: - raise CmdInputError("{} time value must be greater than zero".format(cmd_msg)) - - check_coordinates(xs, ys, zs, G, cmdname, tmp, name='lower') - check_coordinates(xf, yf, zf, G, cmdname, tmp, name='upper') - - if xs >= xf or ys >= yf or zs >= zf: - raise CmdInputError("{} the lower coordinates should be less than the upper coordinates".format(cmd_msg)) - if dx < 0 or dy < 0 or dz < 0: - raise CmdInputError("{} the step size should not be less than zero".format(cmd_msg)) - if dx < G.dx or dy < G.dy or dz < G.dz: - raise CmdInputError("{} the step size should not be less than the spatial discretisation".format(cmd_msg)) - if time <= 0 or time > G.iterations: - raise CmdInputError("{} time value is not valid".format(cmd_msg)) - - 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, dx * G.dy, dx * G.dz, s.time * G.dt, s.basefilename)) - - G.snapshots.append(s) - - -def create_material(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) != 5: - raise CmdInputError("{} requires exactly five parameters".format(cmd_msg)) - if float(tmp[0]) < 0: - raise CmdInputError("{} requires a positive value for static (DC) permittivity".format(cmd_msg)) - if float(tmp[1]) < 0: - raise CmdInputError("{} requires a positive value for conductivity".format(cmd_msg)) - if float(tmp[2]) < 0: - raise CmdInputError("{} requires a positive value for permeability".format(cmd_msg)) - if float(tmp[3]) < 0: - raise CmdInputError("{} requires a positive value for magnetic conductivity".format(cmd_msg)) - if any(x.ID == tmp[4] for x in G.materials): - raise CmdInputError("{} with ID {} already exists'.format(tmp[4]".format(cmd_msg)) - - # Create a new instance of the Material class material (start index after pec & free_space) - m = Material(len(G.materials), tmp[4], G) - m.er = float(tmp[0]) - m.se = float(tmp[1]) - m.mr = float(tmp[2]) - m.sm = float(tmp[3]) - - if G.messages: - print('Material {} with epsr={:g}, sig={:g} S/m; mur={:g}, sig*={:g} S/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) - - -def create_add_dispersion_debye(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - - if len(tmp) < 4: - raise CmdInputError("{} requires at least four parameters".format(cmd_msg)) - if int(tmp[0]) < 0: - raise CmdInputError("{} requires a positive value for number of poles".format(cmd_msg)) - 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("{} material(s) {} do not exist".format(cmd_msg, notfound)) - - for material in materials: - material.type = 'debye' - material.poles = poles - material.average = False - for pole in range(1, 2 * poles, 2): - if float(tmp[pole]) > 0 and float(tmp[pole + 1]) > G.dt: - material.deltaer.append(float(tmp[pole])) - material.tau.append(float(tmp[pole + 1])) - else: - raise CmdInputError("{} requires positive values for the permittivity \ - difference and relaxation times, and relaxation times that are greater \ - than the time step for the model.".format(cmd_msg)) - if material.poles > Material.maxpoles: - Material.maxpoles = material.poles - - if G.messages: - print('Debye disperion added to {} with delta_epsr={}, 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))) - - -def create_add_dispersion_lorentz(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - - if len(tmp) < 5: - raise CmdInputError("{} requires at least five parameters".format(cmd_msg)) - if int(tmp[0]) < 0: - raise CmdInputError("{} requires a positive value for number of poles".format(cmd_msg)) - 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("{} material(s) {} do not exist".format(notfound, cmd_msg)) - - for material in materials: - material.type = 'lorentz' - material.poles = poles - material.average = 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("{} requires positive values for the permittivity difference and frequencies, and associated times that are greater than the time step for the model.".format(cmd_msg)) - if material.poles > Material.maxpoles: - Material.maxpoles = material.poles - - if G.messages: - print('Lorentz disperion added to {} with delta_epsr={}, 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))) - - -def create_add_dispersion_drude(cmdname, cmdinstance, G): - - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - - if len(tmp) < 5: - raise CmdInputError("{} requires at least five parameters".format(cmd_msg)) - if int(tmp[0]) < 0: - raise CmdInputError("{} requires a positive value for number of poles".format(cmd_msg)) - 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("{} material(s) {} do not exist".format(cmd_msg, notfound)) - - for material in materials: - material.type = 'drude' - material.poles = poles - material.average = 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("{} requires positive values for the frequencies, and associated times that are greater than the time step for the model.".format(cmd_msg)) - if material.poles > Material.maxpoles: - Material.maxpoles = material.poles - - if G.messages: - print('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))) - - -def create_soil_peplinski(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) != 7: - raise CmdInputError("{} requires at exactly seven parameters".format(cmd_msg)) - if float(tmp[0]) < 0: - raise CmdInputError("{} requires a positive value for the sand fraction".format(cmd_msg)) - if float(tmp[1]) < 0: - raise CmdInputError("{} requires a positive value for the clay fraction".format(cmd_msg)) - if float(tmp[2]) < 0: - raise CmdInputError("{} requires a positive value for the bulk density".format(cmd_msg)) - if float(tmp[3]) < 0: - raise CmdInputError("{} requires a positive value for the sand particle density".format(cmd_msg)) - if float(tmp[4]) < 0: - raise CmdInputError("{} requires a positive value for the lower limit of the water volumetric fraction".format(cmd_msg)) - if float(tmp[5]) < 0: - raise CmdInputError("{} requires a positive value for the upper limit of the water volumetric fraction".format(cmd_msg)) - if any(x.ID == tmp[6] for x in G.mixingmodels): - raise CmdInputError("{} with ID {} already exists".format(cmd_msg, 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) - - -def create_geometry_view(cmdname, cmdinstance, G): - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - if len(tmp) != 11: - raise CmdInputError("{} requires exactly eleven parameters".format(cmd_msg)) - - 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, G, cmdname, tmp, name='lower') - check_coordinates(xf, yf, zf, G, cmdname, tmp, name='upper') - - if xs >= xf or ys >= yf or zs >= zf: - raise CmdInputError("{} the lower coordinates should be less than the upper coordinates".format(cmd_msg)) - if dx < 0 or dy < 0 or dz < 0: - raise CmdInputError("{} the step size should not be less than zero".format(cmd_msg)) - if dx > G.nx or dy > G.ny or dz > G.nz: - raise CmdInputError("{} the step size should be less than the domain size".format(cmd_msg)) - if dx < G.dx or dy < G.dy or dz < G.dz: - raise CmdInputError("{} the step size should not be less than the spatial discretisation".format(cmd_msg)) - if tmp[10].lower() != 'n' and tmp[10].lower() != 'f': - raise CmdInputError("{} requires type to be either n (normal) or f (fine)".format(cmd_msg)) - 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("{} requires the spatial discretisation for the geometry view to be the same as the model for geometry view of type f (fine)".format(cmd_msg)) - - g = GeometryView(xs, ys, zs, xf, yf, zf, dx, dy, dz, tmp[9], tmp[10].lower()) - - 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, 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, g.basefilename)) - - # Append the new GeometryView object to the geometry views list - G.geometryviews.append(g) - - -def create_pml_cfs(cmdname, cmdinstance, G): - - tmp = cmdinstance.split() - cmd_msg = create_cmd_msg(cmdname, cmdinstance) - - if len(G.cfs) == 2: - raise CmdInputError("{} can only be used up to two times, for up to a 2nd order PML".format(cmd_msg)) - - if len(tmp) != 12: - raise CmdInputError("{} requires exactly twelve parameters".format(cmd_msg)) - 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("{} must have scaling type {}'.format(','.join(CFSParameter.scalingprofiles.keys())".format(cmd_msg)) - if tmp[1] not in CFSParameter.scalingdirections or tmp[5] not in CFSParameter.scalingdirections or tmp[9] not in CFSParameter.scalingdirections: - raise CmdInputError("{} must have scaling type {}'.format(','.join(CFSParameter.scalingprofiles.keys())".format(cmd_msg)) - 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("{} minimum and maximum scaling values must be greater than zero".format(cmd_msg)) - if float(tmp[6]) < 1: - raise CmdInputError("{} minimum scaling value for kappa must be greater than zero".format(cmd_msg)) - - 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: {:g}) 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) +from gprMax.constants import z0 +from gprMax.exceptions import CmdInputError +from gprMax.geometry_views import GeometryView +from gprMax.materials import Material, PeplinskiSoil +from gprMax.pml import CFSParameter, CFS +from gprMax.receivers import Rx +from gprMax.snapshots import Snapshot +from gprMax.sources import VoltageSource, HertzianDipole, MagneticDipole, TransmissionLine +from gprMax.utilities import round_value +from gprMax.waveforms import Waveform def process_multicmds(multicmds, G): - """Checks the validity of command parameters and creates instances of - classes of parameters. + """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. + multicmds (dict): Commands that can have multiple instances in the model. + G (class): Grid class instance - holds essential parameters describing the model. """ - # Order of commands is important as some cmds have dependencies. - cmds = OrderedDict() + # 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) - cmds['#waveform'] = create_waveforms - cmds['#voltage_source'] = create_voltage_source - cmds['#hertzian_dipole'] = create_hertzian_dipole - cmds['#magnetic_dipole'] = create_magnetic_dipole - cmds['#transmission_line'] = create_transmission_line - cmds['#tem_transmission_line'] = create_transmission_line - cmds['#rx'] = create_rx - cmds['#rx_box'] = create_rx_box - cmds['#snapshot'] = create_snapshot - cmds['#material'] = create_material - cmds['#add_dispersion_debye'] = create_add_dispersion_debye - cmds['#add_dispersion_lorentz'] = create_add_dispersion_lorentz - cmds['#add_dispersion_drude'] = create_add_dispersion_drude - cmds['#soil_peplinski'] = create_soil_peplinski - cmds['#geometry_view'] = create_geometry_view - cmds['#pml_cfs'] = create_pml_cfs + # Waveform definitions + cmdname = '#waveform' + if multicmds[cmdname] != '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 amplitude {: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] != '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 + if tmp[0].lower() not in ('x', 'y', 'z'): + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z') + 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[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: + print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n') + 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= tmp[0] + v.xcoord = xcoord + v.ycoord = ycoord + v.zcoord = zcoord + v.ID = 'VoltageSource(' + 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 = ' ' + + 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] != '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 + if tmp[0].lower() not in ('x', 'y', 'z'): + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z') + 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[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: + print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n') + + # 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 = tmp[0] + h.xcoord = xcoord + h.ycoord = ycoord + h.zcoord = zcoord + h.xcoordbase = xcoord + h.ycoordbase = ycoord + h.zcoordbase = zcoord + h.ID = 'HertzianDipole(' + 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 = ' ' + + if G.messages: + 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] != '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 + if tmp[0].lower() not in ('x', 'y', 'z'): + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z') + 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[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: + print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n') + + # 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 = tmp[0] + m.xcoord = xcoord + m.ycoord = ycoord + m.zcoord = zcoord + m.xcoordbase = xcoord + m.ycoordbase = ycoord + m.zcoordbase = zcoord + m.ID = 'MagneticDipole(' + 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 = ' ' + + 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] != '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 + if tmp[0].lower() not in ('x', 'y', 'z'): + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z') + + 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[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: + print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n') + 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[4])) + + t = TransmissionLine(G) + t.polarisation = tmp[0] + t.xcoord = xcoord + t.ycoord = ycoord + t.zcoord = zcoord + t.ID = 'TransmissionLine(' + str(t.xcoord) + ',' + str(t.ycoord) + ',' + str(t.zcoord) + ')' + t.resistance = resistance + t.waveformID = tmp[5] + t.calculate_incident_V_I(G) + + 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 = ' ' + + 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] != '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[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: + print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n') + + r = Rx() + r.xcoord = xcoord + r.ycoord = ycoord + r.zcoord = zcoord + r.xcoordbase = xcoord + r.ycoordbase = ycoord + r.zcoordbase = zcoord + + # If no ID or outputs are specified, use default i.e Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy, Iz + if len(tmp) == 3: + r.ID = 'Rx(' + str(r.xcoord) + ',' + str(r.ycoord) + ',' + str(r.zcoord) + ')' + r.outputs = Rx.availableoutputs[0:9] + else: + r.ID = tmp[3] + # Check and add field output names + for field in tmp[4::]: + if field in Rx.availableoutputs: + r.outputs.append(field) + else: + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' contains an output type that is not available') + + if G.messages: + print('Receiver at {:g}m, {:g}m, {:g}m with output(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs))) + + G.rxs.append(r) + + + # Receiver box + cmdname = '#rx_box' + if multicmds[cmdname] != '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[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: + print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n') + 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.dx or dy < G.dy or dz < G.dz: + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than the spatial discretisation') + + for x in range(xs, xf, dx): + for y in range(ys, yf, dy): + for z in range(zs, zf, dz): + r = Rx() + r.xcoord = x + r.ycoord = y + r.zcoord = z + r.xcoordbase = x + r.ycoordbase = y + r.zcoordbase = z + r.ID = 'Rx(' + str(x) + ',' + str(y) + ',' + str(z) + ')' + G.rxs.append(r) + + if G.messages: + print('Receiver box {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g} 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)) + + + # Snapshot + cmdname = '#snapshot' + if multicmds[cmdname] != '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: + 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 < G.dx or dy < G.dy or dz < G.dz: + 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, dx * G.dy, dx * G.dz, s.time * G.dt, s.basefilename)) + + G.snapshots.append(s) + + # Materials + cmdname = '#material' + if multicmds[cmdname] != '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]) < 0: + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for static (DC) permittivity') + if float(tmp[1]) < 0: + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value for conductivity') + if float(tmp[2]) < 0: + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires a positive value 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], G) + m.er = float(tmp[0]) + m.se = float(tmp[1]) + m.mr = float(tmp[2]) + m.sm = float(tmp[3]) + + if G.messages: + print('Material {} with epsr={:g}, sig={:g} S/m; mur={:g}, sig*={:g} S/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] != '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.average = False + for pole in range(1, 2 * poles, 2): + if float(tmp[pole]) > 0 and float(tmp[pole + 1]) > G.dt: + 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 and relaxation times, and relaxation times that are greater than the time step for the model.') + if material.poles > Material.maxpoles: + Material.maxpoles = material.poles + + if G.messages: + print('Debye disperion added to {} with delta_epsr={}, 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] != '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.average = 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: + print('Lorentz disperion added to {} with delta_epsr={}, 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] != '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.average = 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: + print('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] != '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] != '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 < G.dx or dy < G.dy or dz < G.dz: + 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)') + + g = GeometryView(xs, ys, zs, xf, yf, zf, dx, dy, dz, tmp[9], tmp[10].lower()) + + 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, 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, g.basefilename)) + + # Append the new GeometryView object to the geometry views list + G.geometryviews.append(g) + + + # Complex frequency shifted (CFS) PML parameter + cmdname = '#pml_cfs' + if multicmds[cmdname] != '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.scalingprofiles.keys()))) + 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 greater than zero') + if float(tmp[6]) < 1: + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' minimum scaling value for kappa must be greater than zero') + + 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: {:g}) 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) - for cmdname, func in cmds.items(): - if multicmds[cmdname] != 'None': - for cmdinstance in multicmds[cmdname]: - func(cmdname, cmdinstance, G) diff --git a/gprMax/materials.py b/gprMax/materials.py index 01b65a59..e4590943 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -20,34 +20,28 @@ import numpy as np from gprMax.constants import e0, m0, floattype, complextype -# Look up for material ids -material_ids = { - 'pec': 0, - 'freespace': 1 -} - class Material(object): """Materials, their properties and update coefficients.""" - + # Maximum number of dispersive material poles in a model maxpoles = 0 - + # Types of material types = ['standard', 'debye', 'lorentz', 'drude'] - + # 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, G): """ Args: @@ -55,7 +49,7 @@ class Material(object): ID (str): Name of the material. G (class): Grid class instance - holds essential parameters describing the model. """ - + self.numID = numID self.ID = ID self.type = 'standard' @@ -67,20 +61,20 @@ class Material(object): 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 @@ -92,11 +86,11 @@ class Material(object): # Calculate electric update coefficients 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) @@ -105,7 +99,7 @@ class Material(object): 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 self.type == 'debye': self.w[x] = self.deltaer[x] / self.tau[x] @@ -123,7 +117,7 @@ class Material(object): 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 @@ -131,7 +125,7 @@ class Material(object): 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 @@ -152,7 +146,7 @@ class Material(object): class PeplinskiSoil(object): """Soil objects that are characterised according to a mixing model by Peplinski (http://dx.doi.org/10.1109/36.387598).""" - + def __init__(self, ID, sandfraction, clayfraction, bulkdensity, sandpartdensity, watervolfraction): """ Args: @@ -163,7 +157,7 @@ class PeplinskiSoil(object): sandpartdensity (float): Density of the sand particles in the soil (g/cm3). watervolfraction (float): Two numbers that specify a range for the volumetric water fraction of the soil. """ - + self.ID = ID self.S = sandfraction self.C = clayfraction @@ -174,28 +168,28 @@ class PeplinskiSoil(object): def calculate_debye_properties(self, nbins, G): """Calculates the real and imaginery part of a Debye model for the soil as well as a conductivity. It uses a semi-empirical model (http://dx.doi.org/10.1109/36.387598). - + Args: nbins (int): Number of bins to use to create the different materials. G (class): Grid class instance - holds essential parameters describing the model. """ - + # Debye model properties of water f = 1.3e9 w = 2 * np.pi * f erealw = Material.watereri + ((Material.waterdeltaer) / (1 + (w * Material.watertau)**2)) eimagw = w * Material.watertau * ((Material.waterdeltaer) / (1 + (w * Material.watertau)**2)) - + a = 0.65 # Experimentally derived constant es = (1.01 + 0.44 * self.rs)**2 - 0.062 b1 = 1.2748 - 0.519 * self.S - 0.152 * self.C b2 = 1.33797 - 0.603 * self.S - 0.166 * self.C - + # For frequencies in the range 0.3GHz to 1.3GHz sigf1 = 0.0467 + 0.2204 * self.rb - 0.411 * self.S + 0.6614 * self.C # For frequencies in the range 1.4GHz to 18GHz sigf2 = -1.645 + 1.939 * self.rb - 2.25622 * self.S + 1.594 * self.C - + # Generate a set of bins based on the given volumetric water fraction values mubins = np.linspace(self.mu[0], self.mu[1], nbins + 1) # Generate a range of volumetric water fraction values the mid-point of each bin to make materials from @@ -208,13 +202,13 @@ class PeplinskiSoil(object): er1 = (1 + (self.rb/self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a) - muiter[0]) ** (1/a) # Real part for frequencies in the range 0.3GHz to 1.3GHz er2 = 1.15 * er1 - 0.68 - + # Imaginary part for frequencies in the range 0.3GHz to 1.3GHz eri = er2 - (muiter[0]**(b2/a) * Material.waterdeltaer) - + # Effective conductivity sig = muiter[0]**(b2/a) * ((sigf1 * (self.rs - self.rb)) / (self.rs * muiter[0])) - + # Check to see if the material already exists before creating a new one requiredID = '|{:.4f}|'.format(float(muiter[0])) material = next((x for x in G.materials if x.ID == requiredID), None) @@ -231,7 +225,7 @@ class PeplinskiSoil(object): m.deltaer.append(er2 - m.er) m.tau.append(Material.watertau) G.materials.append(m) - + muiter.iternext() diff --git a/gprMax/sources.py b/gprMax/sources.py index ba4e06f4..c75d6101 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -20,11 +20,9 @@ from copy import deepcopy import numpy as np -from .constants import c, floattype -from .grid import Ix, Iy, Iz -from .utilities import round_value -from .materials import material_ids -from .exceptions import GeneralError +from gprMax.constants import c, floattype +from gprMax.grid import Ix, Iy, Iz +from gprMax.utilities import round_value class Source(object): @@ -45,18 +43,15 @@ class Source(object): class VoltageSource(Source): - """The voltage source can be a hard source if it's resistance is zero, - i.e. the time variation of the specified electric field component is - prescribed. If it's resistance is non-zero it behaves as a resistive - voltage source.""" - + """The voltage source can be a hard source if it's resistance is zero, i.e. the time variation of the specified electric field component is prescribed. If it's resistance is non-zero it behaves as a resistive voltage source.""" + def __init__(self): - super().__init__() + super(Source, self).__init__() self.resistance = None def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G): """Updates electric field values for a voltage source. - + Args: abstime (float): Absolute time. updatecoeffsE (memory view): numpy array of electric field update coefficients. @@ -64,7 +59,7 @@ class VoltageSource(Source): Ex, Ey, Ez (memory view): numpy array of electric field values. G (class): Grid class instance - holds essential parameters describing the model. """ - + if abstime >= self.start and abstime <= self.stop: # Set the time of the waveform evaluation to account for any delay in the start time = abstime - self.start @@ -72,7 +67,7 @@ class VoltageSource(Source): j = self.ycoord k = self.zcoord waveform = next(x for x in G.waveforms if x.ID == self.waveformID) - + if self.polarisation is 'x': if self.resistance != 0: componentID = 'E' + self.polarisation @@ -96,7 +91,7 @@ class VoltageSource(Source): def create_material(self, G): """Create a new material at the voltage source location that adds the voltage source conductivity to the underlying parameters. - + Args: G (class): Grid class instance - holds essential parameters describing the model. """ @@ -105,7 +100,7 @@ class VoltageSource(Source): i = self.xcoord j = self.ycoord k = self.zcoord - + componentID = 'E' + self.polarisation requirednumID = G.ID[G.IDlookup[componentID], i, j, k] material = next(x for x in G.materials if x.numID == requirednumID) @@ -113,7 +108,7 @@ class VoltageSource(Source): newmaterial.ID = material.ID + '+VoltageSource_' + str(self.resistance) newmaterial.numID = len(G.materials) newmaterial.average = False - + # Add conductivity of voltage source to underlying conductivity if self.polarisation == 'x': newmaterial.se += G.dx / (self.resistance * G.dy * G.dz) @@ -128,13 +123,13 @@ class VoltageSource(Source): class HertzianDipole(Source): """The Hertzian dipole is an additive source (electric current density).""" - + def __init__(self): - super().__init__() + super(Source, self).__init__() def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G): """Updates electric field values for a Hertzian dipole. - + Args: abstime (float): Absolute time. updatecoeffsE (memory view): numpy array of electric field update coefficients. @@ -142,7 +137,7 @@ class HertzianDipole(Source): Ex, Ey, Ez (memory view): numpy array of electric field values. G (class): Grid class instance - holds essential parameters describing the model. """ - + if abstime >= self.start and abstime <= self.stop: # Set the time of the waveform evaluation to account for any delay in the start time = abstime - self.start @@ -150,7 +145,7 @@ class HertzianDipole(Source): j = self.ycoord k = self.zcoord waveform = next(x for x in G.waveforms if x.ID == self.waveformID) - + if self.polarisation is 'x': componentID = 'E' + self.polarisation Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dy * G.dz)) @@ -166,13 +161,13 @@ class HertzianDipole(Source): class MagneticDipole(Source): """The magnetic dipole is an additive source (magnetic current density).""" - + def __init__(self): - super().__init__() + super(Source, self).__init__() def update_magnetic(self, abstime, updatecoeffsH, ID, Hx, Hy, Hz, G): """Updates magnetic field values for a magnetic dipole. - + Args: abstime (float): Absolute time. updatecoeffsH (memory view): numpy array of magnetic field update coefficients. @@ -180,7 +175,7 @@ class MagneticDipole(Source): Hx, Hy, Hz (memory view): numpy array of magnetic field values. G (class): Grid class instance - holds essential parameters describing the model. """ - + if abstime >= self.start and abstime <= self.stop: # Set the time of the waveform evaluation to account for any delay in the start time = abstime - self.start @@ -188,7 +183,7 @@ class MagneticDipole(Source): j = self.ycoord k = self.zcoord waveform = next(x for x in G.waveforms if x.ID == self.waveformID) - + if self.polarisation is 'x': Hx[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) @@ -200,60 +195,45 @@ class MagneticDipole(Source): class TransmissionLine(Source): - """The transmission line source is a one-dimensional transmission - line which is attached virtually to a grid cell. - """ - + """The transmission line source is a one-dimensional transmission line which is attached virtually to a grid cell.""" + def __init__(self, G): """ Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (class): Grid class instance - holds essential parameters describing the model. """ - - super().__init__() + + super(Source, self).__init__() self.resistance = None - + # Coefficients for ABC termination of end of the transmission line self.abcv0 = 0 self.abcv1 = 0 - # Spatial step of transmission line (based on magic time step for - # dispersionless behaviour) + # Spatial step of transmission line (based on magic time step for dispersionless behaviour) self.dl = c * G.dt - - # Number of nodes in the transmission line (initially a long line - # to calculate incident voltage and current); - # consider putting ABCs/PML at end + + # Number of nodes in the transmission line (initially a long line to calculate incident voltage and current); consider putting ABCs/PML at end self.nl = round_value(0.667 * G.iterations) - - # Nodal position of the one-way injector excitation in the - # transmission line + + # Nodal position of the one-way injector excitation in the transmission line self.srcpos = 5 - + # Nodal position of where line connects to antenna/main grid self.antpos = 10 - + self.voltage = np.zeros(self.nl, dtype=floattype) self.current = np.zeros(self.nl, dtype=floattype) self.Vinc = np.zeros(G.iterations, dtype=floattype) self.Iinc = np.zeros(G.iterations, dtype=floattype) - - self.type_str = 'TransmissionLine' - - def setID(self): - self.ID = self.type_str + '(' + str(self.xcoord) + ',' + str(self.ycoord) + ',' + str(self.zcoord) + ')' - + def calculate_incident_V_I(self, G): - """Calculates the incident voltage and current with a long length - transmission line not connected to the main grid from: - http://dx.doi.org/10.1002/mop.10415 - + """Calculates the incident voltage and current with a long length transmission line not connected to the main grid from: http://dx.doi.org/10.1002/mop.10415 + Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (class): Grid class instance - holds essential parameters describing the model. """ - + abstime = 0 for timestep in range(G.iterations): self.Vinc[timestep] = self.voltage[self.antpos] @@ -263,36 +243,33 @@ class TransmissionLine(Source): self.update_current(abstime, G) abstime += 0.5 * G.dt - # Shorten number of nodes in the transmission line before use - # with main grid + # Shorten number of nodes in the transmission line before use with main grid self.nl = self.antpos + 1 - + def update_abc(self, G): """Updates absorbing boundary condition at end of the transmission line. - + Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (class): Grid class instance - holds essential parameters describing the model. """ - + h = (c * G.dt - self.dl) / (c * G.dt + self.dl) - + self.voltage[0] = h * (self.voltage[1] - self.abcv0) + self.abcv1 self.abcv0 = self.voltage[0] self.abcv1 = self.voltage[1] def update_voltage(self, time, G): """Updates voltage values along the transmission line. - + Args: time (float): Absolute time. - G (class): Grid class instance - holds essential parameters - describing the model. + G (class): Grid class instance - holds essential parameters describing the model. """ - + # Update all the voltage values along the line self.voltage[1:self.nl] -= self.resistance * (c * G.dt / self.dl) * (self.current[1:self.nl] - self.current[0:self.nl - 1]) - + # Update the voltage at the position of the one-way injector excitation waveform = next(x for x in G.waveforms if x.ID == self.waveformID) self.voltage[self.srcpos] += (c * G.dt / self.dl) * waveform.amp * waveform.calculate_value(time - 0.5 * G.dt, G.dt) @@ -302,13 +279,12 @@ class TransmissionLine(Source): def update_current(self, time, G): """Updates current values along the transmission line. - + Args: time (float): Absolute time. - G (class): Grid class instance - holds essential parameters - describing the model. + G (class): Grid class instance - holds essential parameters describing the model. """ - + # Update all the current values along the line self.current[0:self.nl - 1] -= (1 / self.resistance) * (c * G.dt / self.dl) * (self.voltage[1:self.nl] - self.voltage[0:self.nl - 1]) @@ -317,25 +293,23 @@ class TransmissionLine(Source): self.current[self.srcpos - 1] += (c * G.dt / self.dl) * waveform.amp * waveform.calculate_value(time - 0.5 * G.dt, G.dt) * (1 / self.resistance) def update_electric(self, abstime, Ex, Ey, Ez, G): - """Updates electric field value in the main grid from voltage value - in the transmission line. - + """Updates electric field value in the main grid from voltage value in the transmission line. + Args: abstime (float): Absolute time. Ex, Ey, Ez (memory view): numpy array of electric field values. - G (class): Grid class instance - holds essential parameters - describing the model. + G (class): Grid class instance - holds essential parameters describing the model. """ + if abstime >= self.start and abstime <= self.stop: - # Set the time of the waveform evaluation to account - # for any delay in the start + # Set the time of the waveform evaluation to account for any delay in the start time = abstime - self.start i = self.xcoord j = self.ycoord k = self.zcoord - + self.update_voltage(time, G) - + if self.polarisation is 'x': Ex[i, j, k] = - self.voltage[self.antpos] / G.dx @@ -346,24 +320,21 @@ class TransmissionLine(Source): Ez[i, j, k] = - self.voltage[self.antpos] / G.dz def update_magnetic(self, abstime, Hx, Hy, Hz, G): - """Updates current value in transmission line from magnetic field - values in the main grid. - + """Updates current value in transmission line from magnetic field values in the main grid. + Args: abstime (float): Absolute time. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (class): Grid class instance - holds essential parameters - describing the model. + G (class): Grid class instance - holds essential parameters describing the model. """ - + if abstime >= self.start and abstime <= self.stop: - # Set the time of the waveform evaluation to account - # for any delay in the start + # Set the time of the waveform evaluation to account for any delay in the start time = abstime - self.start i = self.xcoord j = self.ycoord k = self.zcoord - + if self.polarisation is 'x': self.current[self.antpos] = Ix(i, j, k, G.Hy, G.Hz, G) @@ -375,15 +346,3 @@ class TransmissionLine(Source): self.update_current(time, G) - -class TEMTransmissionLine(TransmissionLine): - - def __init__(self, G): - super().__init__(G) - self.type_str = 'TEMTransmissionLine' - - def update_magnetic(self, abstime, Hx, Hy, Hz, G): - pass - - def update_electric(self, abstime, Ex, Ey, Ez, G): - pass