diff --git a/gprMax/_version.py b/gprMax/_version.py index 6fde4770..f424ece2 100644 --- a/gprMax/_version.py +++ b/gprMax/_version.py @@ -1,4 +1,3 @@ # This is where the version number is set and read by setup.py and conf.py (for the docs) __version__ = '3.0.5' - diff --git a/gprMax/constants.py b/gprMax/constants.py index 16ec223b..70c76458 100644 --- a/gprMax/constants.py +++ b/gprMax/constants.py @@ -31,4 +31,4 @@ z0 = np.sqrt(m0 / e0) # Main field arrays use floats (floattype) and complex numbers (complextype) floattype = np.float32 -complextype = np.complex64 \ No newline at end of file +complextype = np.complex64 diff --git a/gprMax/exceptions.py b/gprMax/exceptions.py index 6b087a20..e6165032 100644 --- a/gprMax/exceptions.py +++ b/gprMax/exceptions.py @@ -22,4 +22,4 @@ class GeneralError(ValueError): class CmdInputError(ValueError): """Handles errors in user specified commands. Subclasses the ValueError class.""" - pass \ No newline at end of file + pass diff --git a/gprMax/fractals.py b/gprMax/fractals.py index 271cd708..9f36ce5e 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.py @@ -26,16 +26,16 @@ np.seterr(divide='raise') class FractalSurface(object): """Fractal surfaces.""" - + surfaceIDs = ['xminus', 'xplus', 'yminus', 'yplus', 'zminus', 'zplus'] - + def __init__(self, xs, xf, ys, yf, zs, zf, dimension): """ Args: xs, xf, ys, yf, zs, zf (float): Extent of the fractal surface (one pair of coordinates must be equal to correctly define a surface). dimension (float): Fractal dimension that controls the fractal distribution. """ - + self.ID = None self.surfaceID = None self.xs = xs @@ -62,7 +62,7 @@ class FractalSurface(object): Args: G (class): Grid class instance - holds essential parameters describing the model. """ - + if self.xs == self.xf: surfacedims = (self.ny + 1, self.nz + 1) d = G.dx @@ -74,14 +74,14 @@ class FractalSurface(object): d = G.dz self.fractalsurface = np.zeros(surfacedims, dtype=complextype) - + # Positional vector at centre of array, scaled by weighting v1 = np.array([self.weighting[0]*(surfacedims[0])/2, self.weighting[1]*(surfacedims[1])/2]) - + # 2D array of random numbers to be convolved with the fractal function R = np.random.RandomState(self.seed) A = R.randn(surfacedims[0], surfacedims[1]) - + # 2D FFT A = np.fft.fftn(A) @@ -109,14 +109,14 @@ class FractalSurface(object): class FractalVolume(object): """Fractal volumes.""" - + def __init__(self, xs, xf, ys, yf, zs, zf, dimension): """ Args: xs, xf, ys, yf, zs, zf (float): Extent of the fractal volume. dimension (float): Fractal dimension that controls the fractal distribution. """ - + self.ID = None self.operatingonID= None self.xs = xs @@ -135,26 +135,26 @@ class FractalVolume(object): self.weighting = (1, 1, 1) self.nbins = 0 self.fractalsurfaces = [] - + def generate_fractal_volume(self, G): """Generate a 3D volume with a fractal distribution. Args: G (class): Grid class instance - holds essential parameters describing the model. """ - + self.fractalvolume = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=complextype) - + # Positional vector at centre of array, scaled by weighting v1 = np.array([self.weighting[0]*(self.nx + 1)/2, self.weighting[1]*(self.ny + 1)/2, self.weighting[2]*(self.nz + 1)/2]) - + # 3D array of random numbers to be convolved with the fractal function R = np.random.RandomState(self.seed) A = R.randn(self.nx + 1, self.ny + 1, self.nz + 1) - + # 3D FFT A = np.fft.fftn(A) - + for i in range(self.nx + 1): for j in range(self.ny + 1): for k in range(self.nz + 1): @@ -166,7 +166,7 @@ class FractalVolume(object): except FloatingPointError: rr = 0.9 self.fractalvolume[i, j, k] = A[i, j, k] * 1/(rr**self.b) - + # Shift the zero frequency component to the centre of the spectrum self.fractalvolume = np.fft.ifftshift(self.fractalvolume) # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT @@ -179,7 +179,7 @@ class FractalVolume(object): def generate_volume_mask(self): """Generate a 3D volume to use as a mask for adding rough surfaces, water and grass/roots. Zero signifies the mask is not set, one signifies the mask is set.""" - + self.mask = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.int8) maskxs = self.originalxs - self.xs maskxf = (self.originalxf - self.originalxs) + maskxs + 1 @@ -192,7 +192,7 @@ class FractalVolume(object): class Grass(object): """Geometry information for blades of grass.""" - + def __init__(self, numblades): """ Args: @@ -202,7 +202,7 @@ class Grass(object): self.numblades = numblades self.geometryparams = np.zeros((self.numblades, 6), dtype=floattype) self.seed = None - + # Randomly defined parameters that will be used to calculate geometry self.R1 = np.random.RandomState(self.seed) self.R2 = np.random.RandomState(self.seed) @@ -210,7 +210,7 @@ class Grass(object): self.R4 = np.random.RandomState(self.seed) self.R5 = np.random.RandomState(self.seed) self.R6 = np.random.RandomState(self.seed) - + for i in range(self.numblades): self.geometryparams[i, 0] = 10 + 20 * self.R1.random_sample() self.geometryparams[i, 1] = 10 + 20 * self.R2.random_sample() @@ -245,22 +245,10 @@ class Grass(object): Returns: x, y (float): x and y coordinates of grass root. """ - + self.geometryparams[root, 4] += -1 + 2 * self.R5.random_sample() self.geometryparams[root, 5] += -1 + 2 * self.R6.random_sample() x = round(self.geometryparams[root, 4]) y = round(self.geometryparams[root, 5]) return x, y - - - - - - - - - - - - diff --git a/gprMax/grid.py b/gprMax/grid.py index 698ca4f0..670ae354 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -272,6 +272,3 @@ def Iz(x, y, z, Hx, Hy, G): else: Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z]) return Iz - - - diff --git a/gprMax/input_cmd_funcs.py b/gprMax/input_cmd_funcs.py index 1df597f8..015aa149 100644 --- a/gprMax/input_cmd_funcs.py +++ b/gprMax/input_cmd_funcs.py @@ -394,5 +394,3 @@ def rx_steps(dx=0, dy=0, dz=0): c = Coordinate(dx, dy, dz) command('rx_steps', str(c)) return c - - diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index e86a011b..fbde599c 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/input_cmds_file.py @@ -33,20 +33,20 @@ def process_python_include_code(inputfile, usernamespace): Returns: processedlines (list): Input commands after Python processing. """ - + 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:')): - + # String to hold Python code to be executed pythoncode = '' x += 1 @@ -65,10 +65,10 @@ def process_python_include_code(inputfile, 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 = [] pythonstdout = [] @@ -77,34 +77,34 @@ def process_python_include_code(inputfile, usernamespace): hashcmds.append(line + '\n') elif line: pythonstdout.append(line) - + # Add commands to a list processedlines.extend(hashcmds) - + # Print any generated output that is not commands if pythonstdout: print('Python messages (from stdout): {}\n'.format(pythonstdout)) - + # 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 @@ -112,7 +112,7 @@ def process_python_include_code(inputfile, usernamespace): processedlines.append(inputlines[x]) x += 1 - + return processedlines @@ -125,7 +125,7 @@ def write_processed_file(inputfile, modelrun, numbermodelruns, processedlines): 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: @@ -160,12 +160,12 @@ def check_cmd_names(processedlines, checkessential=True): # 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', '#rx', '#rx_array', '#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 @@ -180,25 +180,24 @@ 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)) return singlecmds, multiplecmds, geometry - diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index e817de47..2cbe5247 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -41,7 +41,7 @@ def process_geometrycmds(geometry, G): for object in geometry: tmp = object.split() - + if tmp[0] == '#geometry_objects_file:': if len(tmp) != 6: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly five parameters') @@ -51,28 +51,28 @@ def process_geometrycmds(geometry, G): 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)) - + # Read materials from file with open(matfile, 'r') as f: # Strip out any newline characters and comments that must begin with double hashes materials = [line.rstrip() + '\n' for line in f if(not line.startswith('##') and line.rstrip('\n'))] - + numexistmaterials = len(G.materials) - + # Check validity of command names singlecmds, multicmds, geometry = check_cmd_names(materials, checkessential=False) - + # Process parameters for commands that can occur multiple times in the model process_multicmds(multicmds, G) - + # 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'] @@ -80,28 +80,28 @@ def process_geometrycmds(geometry, G): 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') - + build_voxels_from_array(xs, ys, zs, numexistmaterials, data, G.solid, G.rigidE, G.rigidH, G.ID) - + if G.messages: print('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)) - - + + 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: @@ -116,12 +116,12 @@ def process_geometrycmds(geometry, G): 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: @@ -146,33 +146,33 @@ def process_geometrycmds(geometry, G): else: for k in range(zs, zf): build_edge_z(xs, ys, k, material.numID, G.rigidE, G.rigidH, G.ID) - + if G.messages: print('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: @@ -187,7 +187,7 @@ def process_geometrycmds(geometry, G): 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: @@ -196,14 +196,14 @@ def process_geometrycmds(geometry, G): 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] @@ -216,7 +216,7 @@ def process_geometrycmds(geometry, G): # Isotropic case if len(materials) == 1: numIDx = numIDy = numIDz = materials[0].numID - + # Uniaxial anisotropic case elif len(materials) == 2: numIDy = materials[0].numID @@ -231,12 +231,12 @@ def process_geometrycmds(geometry, G): # 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) @@ -251,7 +251,7 @@ def process_geometrycmds(geometry, G): 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) @@ -263,7 +263,7 @@ def process_geometrycmds(geometry, G): 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]] @@ -285,7 +285,7 @@ def process_geometrycmds(geometry, G): 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 @@ -296,7 +296,7 @@ def process_geometrycmds(geometry, G): 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: @@ -331,7 +331,7 @@ def process_geometrycmds(geometry, G): if len(materials) == 1: averaging = materials[0].average and averagetriangularprism numID = numIDx = numIDy = numIDz = materials[0].numID - + # Uniaxial anisotropic case elif len(materials) == 3: averaging = False @@ -350,7 +350,7 @@ def process_geometrycmds(geometry, G): 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: @@ -358,7 +358,7 @@ def process_geometrycmds(geometry, G): # 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 @@ -383,7 +383,7 @@ def process_geometrycmds(geometry, G): 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]] @@ -405,14 +405,14 @@ def process_geometrycmds(geometry, G): 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: @@ -434,12 +434,12 @@ def process_geometrycmds(geometry, G): 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].average and averagebox numID = numIDx = numIDy = numIDz = materials[0].numID - + # Uniaxial anisotropic case elif len(materials) == 3: averaging = False @@ -458,12 +458,12 @@ def process_geometrycmds(geometry, G): 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' @@ -475,7 +475,7 @@ def process_geometrycmds(geometry, G): 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]] @@ -497,7 +497,7 @@ def process_geometrycmds(geometry, G): 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 @@ -505,7 +505,7 @@ def process_geometrycmds(geometry, G): 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)) @@ -515,12 +515,12 @@ def process_geometrycmds(geometry, G): 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].average and averagecylinder numID = numIDx = numIDy = numIDz = materials[0].numID - + # Uniaxial anisotropic case elif len(materials) == 3: averaging = False @@ -539,7 +539,7 @@ def process_geometrycmds(geometry, G): 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) @@ -556,7 +556,7 @@ def process_geometrycmds(geometry, G): 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]] @@ -610,7 +610,7 @@ def process_geometrycmds(geometry, G): if len(materials) == 1: averaging = materials[0].average and averagecylindricalsector numID = numIDx = numIDy = numIDz = materials[0].numID - + # Uniaxial anisotropic case elif len(materials) == 3: averaging = False @@ -629,7 +629,7 @@ def process_geometrycmds(geometry, G): 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: @@ -637,7 +637,7 @@ def process_geometrycmds(geometry, G): # 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 @@ -657,7 +657,7 @@ def process_geometrycmds(geometry, G): 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 @@ -680,7 +680,7 @@ def process_geometrycmds(geometry, G): 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]] @@ -702,7 +702,7 @@ def process_geometrycmds(geometry, G): 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) @@ -715,12 +715,12 @@ def process_geometrycmds(geometry, G): 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].average and averagesphere numID = numIDx = numIDy = numIDz = materials[0].numID - + # Uniaxial anisotropic case elif len(materials) == 3: averaging = False @@ -739,7 +739,7 @@ def process_geometrycmds(geometry, G): 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) @@ -762,14 +762,14 @@ def process_geometrycmds(geometry, G): seed = int(tmp[14]) 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: @@ -794,7 +794,7 @@ def process_geometrycmds(geometry, G): 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) @@ -833,7 +833,7 @@ def process_geometrycmds(geometry, G): 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) @@ -842,7 +842,7 @@ def process_geometrycmds(geometry, G): 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: @@ -863,7 +863,7 @@ def process_geometrycmds(geometry, G): 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: @@ -915,10 +915,10 @@ def process_geometrycmds(geometry, G): if fractalrange[1] > G.nz: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot apply fractal surface 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.surfaceID = requestedsurface surface.fractalrange = fractalrange @@ -940,7 +940,7 @@ def process_geometrycmds(geometry, G): 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) @@ -950,7 +950,7 @@ def process_geometrycmds(geometry, G): 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: @@ -967,7 +967,7 @@ def process_geometrycmds(geometry, G): 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: @@ -982,7 +982,7 @@ def process_geometrycmds(geometry, G): 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') @@ -1043,8 +1043,8 @@ def process_geometrycmds(geometry, G): if G.messages: print('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') @@ -1054,7 +1054,7 @@ def process_geometrycmds(geometry, G): 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) @@ -1064,7 +1064,7 @@ def process_geometrycmds(geometry, G): 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: @@ -1083,7 +1083,7 @@ def process_geometrycmds(geometry, G): 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: @@ -1129,10 +1129,10 @@ def process_geometrycmds(geometry, G): 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 @@ -1182,7 +1182,7 @@ def process_geometrycmds(geometry, G): 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) @@ -1194,7 +1194,7 @@ def process_geometrycmds(geometry, G): if G.messages: print('{} 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 @@ -1239,9 +1239,9 @@ def process_geometrycmds(geometry, G): 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': @@ -1456,4 +1456,3 @@ def process_geometrycmds(geometry, G): data = volume.fractalvolume.astype('int16', order='C') build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0, data, G.solid, G.rigidE, G.rigidH, G.ID) - diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 9c3ab03d..39f3772d 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -404,7 +404,7 @@ def process_multicmds(multicmds, G): 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)) @@ -744,4 +744,3 @@ def process_multicmds(multicmds, G): 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) - diff --git a/gprMax/input_cmds_singleuse.py b/gprMax/input_cmds_singleuse.py index 71268c86..57a5f058 100644 --- a/gprMax/input_cmds_singleuse.py +++ b/gprMax/input_cmds_singleuse.py @@ -35,7 +35,7 @@ def process_singlecmds(singlecmds, G): 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' @@ -196,7 +196,7 @@ def process_singlecmds(singlecmds, G): if 2*G.pmlthickness[0] >= G.nx or 2*G.pmlthickness[1] >= G.ny or 2*G.pmlthickness[2] >= G.nz or 2*G.pmlthickness[3] >= G.nx or 2*G.pmlthickness[4] >= G.ny or 2*G.pmlthickness[5] >= G.nz: raise CmdInputError(cmd + ' has too many cells for the domain size') - + # src_steps cmd = '#src_steps' if singlecmds[cmd] != 'None': @@ -238,7 +238,7 @@ def process_singlecmds(singlecmds, G): # 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) @@ -252,9 +252,8 @@ def process_singlecmds(singlecmds, G): w.uservalues = waveformvalues[:] else: w.uservalues = waveformvalues[:,waveform] - + if G.messages: print('User waveform {} created.'.format(w.ID)) G.waveforms.append(w) - diff --git a/gprMax/materials.py b/gprMax/materials.py index 820a4aef..fa6efc68 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -23,32 +23,32 @@ from gprMax.constants import e0, m0, complextype 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): """ Args: numID (int): Numeric identifier of the material. ID (str): Name of the material. """ - + self.numID = numID self.ID = ID self.type = 'standard' @@ -60,20 +60,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 @@ -89,7 +89,7 @@ class Material(object): 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) @@ -98,7 +98,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] @@ -116,7 +116,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 @@ -124,7 +124,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 @@ -145,7 +145,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: @@ -156,7 +156,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 @@ -172,23 +172,23 @@ class PeplinskiSoil(object): 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 @@ -201,13 +201,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) @@ -224,16 +224,5 @@ class PeplinskiSoil(object): m.deltaer.append(er2 - m.er) m.tau.append(Material.watertau) G.materials.append(m) - + muiter.iternext() - - - - - - - - - - - diff --git a/gprMax/optimisation_taguchi.py b/gprMax/optimisation_taguchi.py index bc656eec..783f2e1e 100644 --- a/gprMax/optimisation_taguchi.py +++ b/gprMax/optimisation_taguchi.py @@ -40,17 +40,17 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): if numbermodelruns > 1: raise CmdInputError('When a Taguchi optimisation is being carried out the number of model runs argument is not required') - + inputfileparts = os.path.splitext(inputfile) # Default maximum number of iterations of optimisation to perform (used if the stopping criterion is not achieved) maxiterations = 20 - + # Process Taguchi code blocks in the input file; pass in ordered dictionary to hold parameters to optimise tmp = usernamespace.copy() tmp.update({'optparams': OrderedDict()}) taguchinamespace = taguchi_code_blocks(inputfile, tmp) - + # Extract dictionaries and variables containing initialisation parameters optparams = taguchinamespace['optparams'] fitness = taguchinamespace['fitness'] @@ -62,7 +62,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): # Dictionary to hold history of optmised values of parameters optparamshist = OrderedDict((key, list()) for key in optparams) - + # Import specified fitness function fitness_metric = getattr(importlib.import_module('user_libs.optimisation_taguchi.fitness_functions'), fitness['name']) @@ -75,7 +75,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): print('\tOutput name(s) from model: {}'.format(fitness['args']['outputs'])) print('\tFitness function {} with stopping criterion {:g}'.format(fitness['name'], fitness['stop'])) print('\tMaximum iterations: {:g}'.format(maxiterations)) - + # Initialise arrays and lists to store parameters required throughout optimisation # Lower, central, and upper values for each parameter levels = np.zeros((s, k), dtype=floattype) @@ -91,7 +91,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): # Reset number of model runs to number of experiments numbermodelruns = N usernamespace['number_model_runs'] = numbermodelruns - + # Fitness values for each experiment fitnessvalues = [] @@ -111,14 +111,14 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): os.remove(outputfile) print('\nTaguchi optimisation, iteration {}: {} initial experiments with fitness values {}.'.format(iteration + 1, numbermodelruns, fitnessvalues)) - + # Calculate optimal levels from fitness values by building a response table; update dictionary of parameters with optimal values optparams, levelsopt = calculate_optimal_levels(optparams, levels, levelsopt, fitnessvalues, OA, N, k) # Update dictionary with history of parameters with optimal values for key, value in optparams.items(): optparamshist[key].append(value[0]) - + # Run a confirmation experiment with optimal values numbermodelruns = 1 usernamespace['number_model_runs'] = numbermodelruns @@ -130,7 +130,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): # Rename confirmation experiment output file so that it is retained for each iteraction os.rename(outputfile, os.path.splitext(outputfile)[0] + '_final' + str(iteration + 1) + '.out') - + print('\nTaguchi optimisation, iteration {} completed. History of optimal parameter values {} and of fitness values {}'.format(iteration + 1, dict(optparamshist), fitnessvalueshist, 68*'*')) iteration += 1 @@ -167,14 +167,14 @@ def taguchi_code_blocks(inputfile, taguchinamespace): Returns: processedlines (list): Input commands after Python processing. """ - + 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'))] - + # Store length of dict taglength = len(taguchinamespace) - + x = 0 while(x < len(inputlines)): if(inputlines[x].startswith('#taguchi:')): @@ -187,19 +187,19 @@ def taguchi_code_blocks(inputfile, taguchinamespace): x += 1 if x == len(inputlines): raise CmdInputError('Cannot find the end of the Taguchi code block, i.e. missing #end_taguchi: command.') - + # Compile code for faster execution taguchicompiledcode = compile(taguchicode, '', 'exec') # Execute code block & make available only usernamespace exec(taguchicompiledcode, taguchinamespace) - + x += 1 # Check if any Taguchi code blocks were found if len(taguchinamespace) == taglength: raise CmdInputError('No #taguchi and #end_taguchi code blocks found.') - + return taguchinamespace @@ -217,20 +217,20 @@ def construct_OA(optparams): s (int): Number of levels in OA t (int): Strength of OA """ - + oadirectory = os.path.join(os.path.dirname(os.path.abspath(__file__)), os.pardir, 'user_libs', 'optimisation_taguchi') oadirectory = os.path.abspath(oadirectory) # Properties of the orthogonal array (OA) # Strength t = 2 - + # Number of levels s = 3 - + # Number of parameters to optimise k = len(optparams) - + # Load the appropriate OA if k <= 4: OA = np.load(os.path.join(oadirectory, 'OA_9_4_3_2.npy')) @@ -259,18 +259,18 @@ def construct_OA(optparams): else: # THIS CASE NEEDS FURTHER TESTING print('\nTaguchi optimisation, WARNING: Optimising more than 7 parameters is currently an experimental feature!') - + p = int(np.ceil(np.log(k * (s - 1) + 1) / np.log(s))) - + # Number of experiments N = s**p # Number of columns cols = int((N - 1) / (s - 1)) - + # Algorithm to construct OA from: http://ieeexplore.ieee.org/xpl/articleDetails.jsp?reload=true&arnumber=6812898 OA = np.zeros((N + 1, cols + 1), dtype=np.int8) - + # Construct basic columns for ii in range(1, p + 1): col = int((s**(ii - 1) - 1) / (s - 1) + 1) @@ -340,7 +340,7 @@ def calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, le else: levels[0, p] = levels[1, p] - levelsdiff[p] levels[2, p] = levels[1, p] + levelsdiff[p] - + # Update dictionary of parameters to optimise with lists of new values; clear dictionary first optparams = OrderedDict((key, list()) for key in optparams) @@ -399,11 +399,11 @@ def calculate_optimal_levels(optparams, levels, levelsopt, fitnessvalues, OA, N, # Calculate optimal level from table of responses optlevel = np.where(responses == np.amax(responses))[0] - + # If 2 experiments produce the same fitness value pick first level (this shouldn't happen if the fitness function is designed correctly) if len(optlevel) > 1: optlevel = optlevel[0] - + levelsopt[p] = optlevel # Update dictionary of parameters to optimise with lists of new values; clear dictionary first @@ -425,7 +425,7 @@ def plot_optimisation_history(fitnessvalueshist, optparamshist, optparamsinit): """ import matplotlib.pyplot as plt - + # Plot history of fitness values fig, ax = plt.subplots(subplot_kw=dict(xlabel='Iterations', ylabel='Fitness value'), num='History of fitness values', figsize=(20, 10), facecolor='w', edgecolor='w') iterations = np.arange(1, len(fitnessvalueshist) + 1) @@ -443,6 +443,3 @@ def plot_optimisation_history(fitnessvalueshist, optparamshist, optparamsinit): ax.grid() p += 1 plt.show() - - - \ No newline at end of file diff --git a/gprMax/pml.py b/gprMax/pml.py index a906818b..52ea5b36 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -25,11 +25,11 @@ from gprMax.pml_2order_update import * class CFSParameter(object): """Individual CFS parameter (e.g. alpha, kappa, or sigma).""" - + # Allowable scaling profiles and directions scalingprofiles = {'constant': 0, 'linear': 1, 'quadratic': 2, 'cubic': 3, 'quartic': 4} scalingdirections = ['forward', 'reverse'] - + def __init__(self, ID=None, scaling='polynomial', scalingprofile=None, scalingdirection='forward', min=0, max=0): """ Args: @@ -40,7 +40,7 @@ class CFSParameter(object): min (float): Minimum value for parameter. max (float): Maximum value for parameter. """ - + self.ID = ID self.scaling = scaling self.scalingprofile = scalingprofile @@ -51,7 +51,7 @@ class CFSParameter(object): class CFS(object): """CFS term for PML.""" - + def __init__(self): """ Args: @@ -59,11 +59,11 @@ class CFS(object): 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, direction, er, mr, G): """Calculates an optimum value for sigma max based on underlying material properties. @@ -73,14 +73,14 @@ class CFS(object): mr (float): Average permeability of underlying material. G (class): Grid class instance - holds essential parameters describing the model. """ - + if direction[0] == 'x': d = G.dx elif direction[0] == 'y': d = G.dy elif direction[0] == 'z': d = G.dz - + # 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)) @@ -97,7 +97,7 @@ class CFS(object): 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] @@ -114,14 +114,14 @@ class CFS(object): Evalues (float): numpy array holding profile value for electric PML update. Hvalues (float): numpy array holding profile value for magnetic PML update. """ - + 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': @@ -132,26 +132,26 @@ class CFS(object): elif parameter.ID == 'sigma': Evalues *= self.sigma.max Hvalues *= self.sigma.max - + if parameter.scalingdirection == 'reverse': Evalues = Evalues[::-1] Hvalues = Hvalues[::-1] - + return Evalues, Hvalues class PML(object): """PML - the implementation comes from the derivation in: http://dx.doi.org/10.1109/TAP.2011.2180344""" - + directions = {0: 'xminus', 1: 'yminus', 2: 'zminus', 3: 'xplus', 4: 'yplus', 5: 'zplus'} - + def __init__(self, G, direction=None, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0): """ Args: xs, xf, ys, yf, zs, zf (float): Extent of the PML volume. cfs (list): CFS class instances associated with the PML. """ - + self.direction = direction self.xs = xs self.xf = xf @@ -166,10 +166,10 @@ class PML(object): if not self.CFS: self.CFS = [CFS()] self.initialise_field_arrays() - + def initialise_field_arrays(self): """Initialise arrays to store fields in PML.""" - + # Subscript notation, e.g. 'EPhiyxz' means the electric field Phi vector, of which the # component being corrected is y, the stretching direction is x, and field derivative # is z direction. @@ -200,7 +200,7 @@ class PML(object): 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 + 1), dtype=floattype) self.ERB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.ERE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) @@ -209,7 +209,7 @@ class PML(object): self.HRB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.HRE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.HRF = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) - + for x, cfs in enumerate(self.CFS): if not cfs.sigma.max: cfs.calculate_sigmamax(self.direction, er, mr, G) @@ -223,7 +223,7 @@ class PML(object): 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) - + # Magnetic PML update coefficients tmp = (2*e0*Hkappa) + G.dt * (Halpha * Hkappa + Hsigma) self.HRA[x, :] = (2*e0 + G.dt*Halpha) / tmp @@ -236,7 +236,7 @@ def build_pmls(G): """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). """ - + if G.messages: print('') @@ -245,7 +245,7 @@ def build_pmls(G): sumer = 0 # Sum of relative permittivities in PML slab summr = 0 # Sum of relative permeabilities in PML slab pmldirection=PML.directions[index] - + if pmldirection[0] == 'x': if pmldirection == 'xminus': pml = PML(G, direction=pmldirection, xf=pmlthickness, yf=G.ny, zf=G.nz) @@ -260,7 +260,7 @@ def build_pmls(G): summr += material.mr averageer = sumer / (G.ny * G.nz) averagemr = summr / (G.ny * G.nz) - + elif pmldirection[0] == 'y': if pmldirection == 'yminus': pml = PML(G, direction=pmldirection, yf=pmlthickness, xf=G.nx, zf=G.nz) @@ -290,12 +290,12 @@ def build_pmls(G): summr += material.mr averageer = sumer / (G.nx * G.ny) averagemr = summr / (G.nx * G.ny) - + if G.messages and G.pmlthickness.count(pmlthickness) != len(G.pmlthickness): print('PML slab ({} direction) with {} cells created.'.format(pml.direction, pml.thickness)) pml.calculate_update_coeffs(averageer, averagemr, G) - + # Where all the thicknesses of all the PML slabs are equal if G.messages and G.pmlthickness.count(G.pmlthickness[0]) == len(G.pmlthickness): if G.pmlthickness[0] == 0: @@ -392,4 +392,3 @@ def update_magnetic_pml(G): elif pml.direction == 'zplus': update_pml_2order_hx_zplus(pml.xs, pml.xf, pml.ys, pml.yf, pml.zs, pml.zf, G.nthreads, G.updatecoeffsH, G.ID, G.Hx, G.Ey, pml.HPhixzy, pml.HRA, pml.HRB, pml.HRE, pml.HRF, G.dz) update_pml_2order_hy_zplus(pml.xs, pml.xf, pml.ys, pml.yf, pml.zs, pml.zf, G.nthreads, G.updatecoeffsH, G.ID, G.Hy, G.Ex, pml.HPhiyzx, pml.HRA, pml.HRB, pml.HRE, pml.HRF, G.dz) - diff --git a/gprMax/receivers.py b/gprMax/receivers.py index 29dc29e5..dc2f6a5d 100644 --- a/gprMax/receivers.py +++ b/gprMax/receivers.py @@ -23,11 +23,11 @@ from gprMax.grid import Ix, Iy, Iz class Rx(object): """Receiver output points.""" - + availableoutputs = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz'] - + def __init__(self): - + self.ID = None self.outputs = OrderedDict() self.xcoord = None @@ -69,4 +69,4 @@ def store_outputs(timestep, Ex, Ey, Ez, Hx, Hy, Hz, G): for tlindex, tl in enumerate(G.transmissionlines): tl.Vtotal[timestep] = tl.voltage[tl.antpos] - tl.Itotal[timestep] = tl.current[tl.antpos] \ No newline at end of file + tl.Itotal[timestep] = tl.current[tl.antpos] diff --git a/gprMax/snapshots.py b/gprMax/snapshots.py index 7dd7accb..3984ea21 100644 --- a/gprMax/snapshots.py +++ b/gprMax/snapshots.py @@ -29,7 +29,7 @@ from gprMax.utilities import round_value class Snapshot(object): """Snapshots of the electric and magnetic field values.""" - + # Set string for byte order if sys.byteorder == 'little': byteorder = 'LittleEndian' @@ -52,7 +52,7 @@ class Snapshot(object): time (int): Iteration number to take the snapshot on. filename (str): Filename to save to. """ - + self.xs = xs self.ys = ys self.zs = zs @@ -73,22 +73,22 @@ class Snapshot(object): numbermodelruns (int): Total number of model runs. G (class): Grid class instance - holds essential parameters describing the model. """ - + # No Python 3 support for VTK at time of writing (03/2015) self.vtk_nx = self.xf - self.xs self.vtk_ny = self.yf - self.ys self.vtk_nz = self.zf - self.zs - + # Create directory and construct filename from user-supplied name and model run number if numbermodelruns == 1: snapshotdir = os.path.join(G.inputdirectory, os.path.splitext(G.inputfilename)[0] + '_snaps') else: snapshotdir = os.path.join(G.inputdirectory, os.path.splitext(G.inputfilename)[0] + '_snaps' + str(modelrun)) - + if not os.path.exists(snapshotdir): os.mkdir(snapshotdir) self.filename = os.path.abspath(os.path.join(snapshotdir, self.basefilename + '.vti')) - + # Calculate number of cells according to requested sampling self.vtk_xscells = round_value(self.xs / self.dx) self.vtk_xfcells = round_value(self.xf / self.dx) @@ -148,4 +148,3 @@ class Snapshot(object): self.filehandle.write('\n\n'.encode('utf-8')) self.filehandle.close() - diff --git a/gprMax/sources.py b/gprMax/sources.py index 209c5e15..66b3e63d 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -44,7 +44,7 @@ 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.""" - + def __init__(self): super(Source, self).__init__() self.resistance = None @@ -59,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 @@ -67,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 @@ -100,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) @@ -108,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) @@ -123,7 +123,7 @@ class VoltageSource(Source): class HertzianDipole(Source): """The Hertzian dipole is an additive source (electric current density).""" - + def __init__(self): super(Source, self).__init__() @@ -137,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 @@ -145,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)) @@ -161,7 +161,7 @@ class HertzianDipole(Source): class MagneticDipole(Source): """The magnetic dipole is an additive source (magnetic current density).""" - + def __init__(self): super(Source, self).__init__() @@ -175,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 @@ -183,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)) @@ -196,46 +196,46 @@ class MagneticDipole(Source): class TransmissionLine(Source): """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. """ - + 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) 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 self.nl = round_value(0.667 * G.iterations) - + # 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.Vtotal = np.zeros(G.iterations, dtype=floattype) self.Itotal = np.zeros(G.iterations, dtype=floattype) - + 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 Args: 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] @@ -247,16 +247,16 @@ class TransmissionLine(Source): # 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. """ - + 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] @@ -268,10 +268,10 @@ class TransmissionLine(Source): time (float): Absolute time. 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) @@ -286,7 +286,7 @@ class TransmissionLine(Source): time (float): Absolute time. 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]) @@ -302,16 +302,16 @@ class TransmissionLine(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 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 @@ -329,14 +329,14 @@ class TransmissionLine(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 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) @@ -347,4 +347,3 @@ class TransmissionLine(Source): self.current[self.antpos] = Iz(i, j, k, G.Hx, G.Hy, G) self.update_current(time, G) - diff --git a/gprMax/utilities.py b/gprMax/utilities.py index 19d23b40..1bd9c469 100644 --- a/gprMax/utilities.py +++ b/gprMax/utilities.py @@ -26,7 +26,7 @@ def logo(version): Args: version (str): Version number. """ - + licenseinfo = """ Copyright (C) 2015-2016: The University of Edinburgh Authors: Craig Warren and Antonis Giannopoulos @@ -43,7 +43,7 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with gprMax. If not, see .""" - + gprMaxlogo = """ __ __ __ _ _ __ _ __| \/ | __ ___ __ @@ -52,10 +52,10 @@ along with gprMax. If not, see .""" \__, | .__/|_| |_| |_|\__,_/_/\_\\ |___/|_| """ - + width = 65 url = 'www.gprmax.com' - + print('\nElectromagnetic modelling software based on the Finite-Difference \nTime-Domain (FDTD) method') print('\n{} {} {}'.format('*'*round((width - len(url))/2), url, '*'*round((width - len(url))/2))) print('{}'.format(gprMaxlogo)) @@ -69,7 +69,7 @@ def update_progress(progress): Args: progress (float): Number between zero and one to signify progress. """ - + # Modify this to change the length of the progress bar barLength = 50 block = round_value(barLength * progress) @@ -88,11 +88,11 @@ def round_value(value, decimalplaces=0): Returns: rounded (int/float): Rounded value. """ - + # Rounds to nearest integer (half values are rounded downwards) if decimalplaces == 0: rounded = int(d.Decimal(value).quantize(d.Decimal('1'),rounding=d.ROUND_HALF_DOWN)) - + # Rounds down to nearest float represented by number of decimal places else: precision = '1.{places}'.format(places='0' * decimalplaces) @@ -116,7 +116,7 @@ def human_size(size, a_kilobyte_is_1024_bytes=True): 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 @@ -124,5 +124,3 @@ def human_size(size, a_kilobyte_is_1024_bytes=True): return '{0:.1f}{1}'.format(size, suffix) raise ValueError('Number is too large.') - - diff --git a/gprMax/waveforms.py b/gprMax/waveforms.py index 2e612acd..57e26b82 100644 --- a/gprMax/waveforms.py +++ b/gprMax/waveforms.py @@ -23,9 +23,9 @@ from gprMax.utilities import round_value class Waveform(object): """Definitions of waveform shapes that can be used with sources.""" - + types = ['gaussian', 'gaussiandot', 'gaussiandotnorm', 'gaussiandotdot', 'gaussiandotdotnorm', 'ricker', 'sine', 'contsine', 'impulse', 'user'] - + def __init__(self): self.ID = None self.type = None @@ -43,7 +43,7 @@ class Waveform(object): Returns: waveform (float): Calculated value for waveform. """ - + # Coefficients for certain waveforms if self.type == 'gaussian' or self.type == 'gaussiandot' or self.type == 'gaussiandotdot': chi = 1 / self.freq @@ -53,21 +53,21 @@ class Waveform(object): chi = np.sqrt(2) / self.freq zeta = np.pi**2 * self.freq**2 delay = time - chi - + # Waveforms if self.type == 'gaussian': waveform = np.exp(-zeta * delay**2) - + elif self.type == 'gaussiandot': waveform = -2 * zeta * delay * np.exp(-zeta * delay**2) - + elif self.type == 'gaussiandotnorm': normalise = np.sqrt(np.exp(1) / (2 * zeta)) waveform = -2 * zeta * delay * np.exp(-zeta * delay**2) * normalise - + elif self.type == 'gaussiandotdot': waveform = 2 * zeta * (2 * zeta * delay**2 - 1) * np.exp(-zeta * delay**2) - + elif self.type == 'gaussiandotdotnorm': normalise = 1 / (2 * zeta) waveform = 2 * zeta * (2 * zeta * delay**2 - 1) * np.exp(-zeta * delay**2) * normalise @@ -80,7 +80,7 @@ class Waveform(object): waveform = np.sin(2 * np.pi * self.freq * time) if time * self.freq > 1: waveform = 0 - + elif self.type == 'contsine': rampamp = 0.25 ramp = rampamp * time * self.freq @@ -94,7 +94,7 @@ class Waveform(object): waveform = 1 elif time >= dt: waveform = 0 - + elif self.type == 'user': index = round_value(time / dt) # Check to see if there are still user specified values and if not use zero @@ -102,5 +102,5 @@ class Waveform(object): waveform = 0 else: waveform = self.uservalues[index] - - return waveform \ No newline at end of file + + return waveform diff --git a/gprMax/writer_hdf5.py b/gprMax/writer_hdf5.py index af5b025e..2ce2ffa0 100644 --- a/gprMax/writer_hdf5.py +++ b/gprMax/writer_hdf5.py @@ -91,4 +91,3 @@ def write_hdf5(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G): f['/rxs/rx' + str(rxindex + 1) + '/Iy'] = rx.outputs['Iy'] if 'Iz' in rx.outputs: f['/rxs/rx' + str(rxindex + 1) + '/Iz'] = rx.outputs['Iz'] -