diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index 48dfe479..a55560b5 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -58,7 +58,7 @@ class GeometryView(object): self.dz = dz self.basefilename = filename self.fileext = fileext - + if self.fileext == '.vti': # Calculate number of cells according to requested sampling for geometry view self.vtk_xscells = round_value(self.xs / self.dx) @@ -71,7 +71,7 @@ class GeometryView(object): self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells self.datawritesize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + 2 * (int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)) - + elif self.fileext == '.vtp': self.vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1) self.vtk_numpoint_components = 3 diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index 3602889a..32eb7097 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -391,7 +391,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): if all(value == 0 for value in G.pmlthickness.values()): if G.messages: print('PML boundaries: switched off') - pass # If all the PMLs are switched off don't need to build anything + pass # If all the PMLs are switched off don't need to build anything else: if G.messages: if all(value == G.pmlthickness['xminus'] for value in G.pmlthickness.values()): @@ -482,7 +482,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): geometryview.write_vtk(modelrun, numbermodelruns, G, pbar) pbar.close() if G.geometryobjectswrite: - + for i, geometryobject in enumerate(G.geometryobjectswrite): pbar = tqdm(total=geometryobject.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry object file {} of {}, {}'.format(i + 1, len(G.geometryobjectswrite), os.path.split(geometryobject.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable) geometryobject.write_hdf5(G, pbar) diff --git a/gprMax/grid.py b/gprMax/grid.py index 4cfe56d3..b67f0bbb 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -78,12 +78,12 @@ class FDTDGrid(Grid): self.title = '' self.messages = True self.tqdmdisable = False - - # Threshold (dB) down from maximum power (0dB) of main frequency used to calculate highest frequency for disperion analysis + + # Threshold (dB) down from maximum power (0dB) of main frequency used to calculate highest frequency for disperion analysis self.highestfreqthres = 60 # Maximum allowable percentage physical phase-velocity phase error self.maxnumericaldisp = 2 - + self.nx = 0 self.ny = 0 self.nz = 0 @@ -157,7 +157,7 @@ def dispersion_analysis(G): # Physical phase velocity error (percentage); grid sampling density; material with maximum permittivity; maximum frequency of interest results = {'deltavp': False, 'N': False, 'material': False, 'maxfreq': False} - + # Find maximum frequency maxfreqs = [] for waveform in G.waveforms: @@ -221,19 +221,19 @@ def dispersion_analysis(G): # Minimum wavelength minwavelength = minvelocity / results['maxfreq'] - + # Maximum spatial step delta = max(G.dx, G.dy, G.dz) - + # Courant stability factor S = (c * G.dt) / delta - + # Grid sampling density results['N'] = minwavelength / delta - + # Numerical phase velocity vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N']))) - + # Physical phase velocity error (percentage) results['deltavp'] = (((vp * c) - c) / c) * 100 diff --git a/gprMax/input_cmd_funcs.py b/gprMax/input_cmd_funcs.py index f57153c0..ecb64bda 100644 --- a/gprMax/input_cmd_funcs.py +++ b/gprMax/input_cmd_funcs.py @@ -38,7 +38,7 @@ Coordinate_tuple = namedtuple('Coordinate', ['x', 'y', 'z']) class Coordinate(Coordinate_tuple): """Subclass of a namedtuple where __str__ outputs 'x y z'""" - + def __str__(self): return '{:g} {:g} {:g}'.format(self.x, self.y, self.z) @@ -53,7 +53,7 @@ def command(cmd, *parameters): Returns: s (str): the printed string """ - + # remove Nones filtered = filter(lambda x: x is not None, parameters) # convert to str @@ -76,11 +76,11 @@ def command(cmd, *parameters): def rotate90_point(x, y, rotate90origin=()): """Rotates a point 90 degrees CCW in the x-y plane. - + Args: x, y (float): Coordinates. rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - + Returns: xrot, yrot (float): Rotated coordinates. """ @@ -158,7 +158,7 @@ def domain(x, y, z): Returns: domain (Coordinate): Namedtuple of the extent of the domain. """ - + domain = Coordinate(x, y, z) command('domain', domain) @@ -174,7 +174,7 @@ def dx_dy_dz(x, y, z): Returns: dx_dy_dz (float): Tuple of the spatial resolutions. """ - + dx_dy_dz = Coordinate(x, y, z) command('dx_dy_dz', dx_dy_dz) @@ -190,7 +190,7 @@ def time_window(time_window): Returns: time_window (float): Duration of simulation. """ - + command('time_window', time_window) return time_window @@ -206,7 +206,7 @@ def material(permittivity, conductivity, permeability, magconductivity, name): magconductivity (float): Magnetic loss of the material. name (str): Material identifier. """ - + command('material', permittivity, conductivity, permeability, magconductivity, name) @@ -222,7 +222,7 @@ def geometry_view(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, type='n'): Returns: s, f, d (tuple): 3 namedtuple Coordinate for the start, finish coordinates and spatial discretisation """ - + s = Coordinate(xs, ys, zs) f = Coordinate(xf, yf, zf) d = Coordinate(dx, dy, dz) @@ -243,7 +243,7 @@ def snapshot(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, filename): Returns: s, f, d (tuple): 3 namedtuple Coordinate for the start, finish coordinates and spatial discretisation """ - + s = Coordinate(xs, ys, zs) f = Coordinate(xf, yf, zf) d = Coordinate(dx, dy, dz) @@ -269,14 +269,13 @@ def edge(xs, ys, zs, xf, yf, zf, material, rotate90origin=()): Returns: s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates """ - + if rotate90origin: if xs == xf: polarisation = 'y' else: polarisation = 'x ' xs, ys, xf, yf = rotate90_edge(xs, ys, xf, yf, polarisation, rotate90origin) - s = Coordinate(xs, ys, zs) f = Coordinate(xf, yf, zf) @@ -296,10 +295,10 @@ def plate(xs, ys, zs, xf, yf, zf, material, rotate90origin=()): Returns: s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates """ - + if rotate90origin: xs, ys, xf, yf = rotate90_plate(xs, ys, xf, yf, rotate90origin) - + s = Coordinate(xs, ys, zs) f = Coordinate(xf, yf, zf) command('plate', s, f, material) @@ -320,12 +319,12 @@ def triangle(x1, y1, z1, x2, y2, z2, x3, y3, z3, thickness, material, averaging= Returns: v1, v2, v3 (tuple): 3 namedtuple Coordinate for the vertices """ - + if rotate90origin: x1, y1 = rotate90_point(x1, y1, rotate90origin) x2, y2 = rotate90_point(x2, y2, rotate90origin) x3, y3 = rotate90_point(x3, y3, rotate90origin) - + v1 = Coordinate(x1, y1, z1) v2 = Coordinate(x2, y2, z2) v3 = Coordinate(x3, y3, z3) @@ -346,10 +345,10 @@ def box(xs, ys, zs, xf, yf, zf, material, averaging='', rotate90origin=()): Returns: s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates """ - + if rotate90origin: xs, ys, xf, yf = rotate90_plate(xs, ys, xf, yf, rotate90origin) - + s = Coordinate(xs, ys, zs) f = Coordinate(xf, yf, zf) command('box', s, f, material, averaging) @@ -369,7 +368,7 @@ def sphere(x, y, z, radius, material, averaging=''): Returns: c (tuple): namedtuple Coordinate for the center of the sphere """ - + c = Coordinate(x, y, z) command('sphere', c, radius, material, averaging) @@ -389,11 +388,11 @@ def cylinder(x1, y1, z1, x2, y2, z2, radius, material, averaging='', rotate90ori Returns: c1, c2 (tuple): 2 namedtuple Coordinate for the centres of the two faces of the cylinder. """ - + if rotate90origin: x1, y1 = rotate90_point(x1, y1, rotate90origin) x2, y2 = rotate90_point(x2, y2, rotate90origin) - + c1 = Coordinate(x1, y1, z1) c2 = Coordinate(x2, y2, z2) command('cylinder', c1, c2, radius, material, averaging) @@ -414,7 +413,7 @@ def cylindrical_sector(axis, ctr1, ctr2, t1, t2, radius, startingangle, sweptang material (str): Material identifier(s). averaging (str): Turn averaging on or off. """ - + command('cylindrical_sector', axis, ctr1, ctr2, t1, t2, radius, startingangle, sweptangle, material, averaging) @@ -427,7 +426,7 @@ def excitation_file(file1): Returns: file1 (str): filename """ - + command('excitation_file', file1) return file1 @@ -445,7 +444,7 @@ def waveform(shape, amplitude, frequency, identifier): Returns: identifier (str): is an identifier for the waveform used to assign it to a source. """ - + command('waveform', shape, amplitude, frequency, identifier) return identifier @@ -466,7 +465,7 @@ def hertzian_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None Returns: coordinates (tuple): namedtuple Coordinate of the source location """ - + if rotate90origin: if polarisation == 'x': xf = f1 + dxdy[0] @@ -476,7 +475,7 @@ def hertzian_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None xf = f1 yf = f2 + dxdy[1] newpolarisation = 'x' - + f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) polarisation = newpolarisation @@ -502,7 +501,7 @@ def magnetic_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None Returns: coordinates (tuple): namedtuple Coordinate of the source location """ - + if rotate90origin: if polarisation == 'x': xf = f1 + dxdy[0] @@ -512,7 +511,7 @@ def magnetic_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None xf = f1 yf = f2 + dxdy[1] newpolarisation = 'x' - + f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) polarisation = newpolarisation @@ -539,7 +538,7 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_ Returns: coordinates (tuple): namedtuple Coordinate of the source location """ - + if rotate90origin: if polarisation == 'x': xf = f1 + dxdy[0] @@ -549,7 +548,7 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_ xf = f1 yf = f2 + dxdy[1] newpolarisation = 'x' - + f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) polarisation = newpolarisation @@ -562,7 +561,7 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_ def transmission_line(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #transmission_line: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove] - + Args: polarisation (str): is the polarisation of the source and can be 'x', 'y', or 'z'. f1 f2 f3 (float): are the coordinates (x,y,z) of the source in the model. @@ -572,11 +571,11 @@ def transmission_line(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_remove (float): is a time to remove the source. dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only. rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - + Returns: coordinates (tuple): namedtuple Coordinate of the source location """ - + if rotate90origin: if polarisation == 'x': xf = f1 + dxdy[0] @@ -586,14 +585,14 @@ def transmission_line(polarisation, f1, f2, f3, resistance, identifier, t0=None, xf = f1 yf = f2 + dxdy[1] newpolarisation = 'x' - + f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) polarisation = newpolarisation - + c = Coordinate(f1, f2, f3) # since command ignores None, this is safe: command('transmission_line', polarisation, str(c), resistance, identifier, t0, t_remove) - + return c @@ -611,7 +610,7 @@ def rx(x, y, z, identifier=None, to_save=None, polarisation=None, dxdy=None, rot Returns: coordinates (tuple): namedtuple Coordinate of the receiver location """ - + if rotate90origin: if polarisation == 'x': xf = x + dxdy[0] @@ -620,7 +619,7 @@ def rx(x, y, z, identifier=None, to_save=None, polarisation=None, dxdy=None, rot xf = x yf = y + dxdy[1] x, y, xf, yf = rotate90_edge(x, y, xf, yf, polarisation, rotate90origin) - + c = Coordinate(x, y, z) command('rx', str(c), identifier, ' '.join(to_save)) diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index df22ae39..46645ab1 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/input_cmds_file.py @@ -173,7 +173,7 @@ def check_cmd_names(processedlines, checkessential=True): cmd = processedlines[lindex].split(':') cmdname = cmd[0] cmdparams = cmd[1] - + # Check if there is space between command name and parameters, i.e. check first character of parameter string if ' ' not in cmdparams[0]: raise CmdInputError('There must be a space between the command name and parameters in ' + processedlines[lindex]) diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index 4eec939b..91dc10b7 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -65,7 +65,7 @@ def process_geometrycmds(geometry, G): matstr = os.path.splitext(os.path.split(matfile)[1])[0] numexistmaterials = len(G.materials) - + # Read materials from file with open(matfile, 'r') as f: # Strip out any newline characters and comments that must begin with double hashes @@ -73,7 +73,7 @@ def process_geometrycmds(geometry, G): # Check validity of command names singlecmdsimport, multicmdsimport, geometryimport = check_cmd_names(materials, checkessential=False) - + # Process parameters for commands that can occur multiple times in the model process_multicmds(multicmdsimport, G) @@ -873,7 +873,7 @@ def process_geometrycmds(geometry, G): # Search and process any modifiers for the fractal box for object in geometry: tmp = object.split() - + if tmp[0] == '#add_surface_roughness:': if len(tmp) < 13: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least twelve parameters') diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index bae08a76..a44302ed 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -502,11 +502,11 @@ def process_multicmds(multicmds, G): m.se = se m.mr = float(tmp[2]) m.sm = float(tmp[3]) - + # Set material averaging to False if infinite conductivity, i.e. pec if m.se == float('inf'): m.averagable = False - + if G.messages: tqdm.write('Material {} with epsr={:g}, sig={:g} S/m; mur={:g}, sig*={:g} S/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm)) diff --git a/gprMax/materials.py b/gprMax/materials.py index 5fabbb1a..06fe2858 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -156,7 +156,7 @@ def process_materials(G): else: materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]', '\nDelta eps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]', 'gamma\n[Hz]', '\nmu_r', 'sigma*\n[S/m]', 'Dielectric\nsmoothable']] - for material in G.materials: + for material in G.materials: # Calculate update coefficients for material material.calculate_update_coeffsE(G) material.calculate_update_coeffsH(G) @@ -181,23 +181,23 @@ def process_materials(G): materialtext.append('{:g}'.format(material.se)) if Material.maxpoles > 0: if 'debye' in material.type: - materialtext.append(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer)) - materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau)) + materialtext.append(', '.join('{:g}'.format(deltaer for deltaer in material.deltaer))) + materialtext.append(', '.join('{:g}'.format(tau for tau in material.tau))) materialtext.append('') materialtext.append('') materialtext.append('') elif 'lorentz' in material.type: - materialtext.append(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer)) + materialtext.append(', '.join('{:g}'.format(deltaer for deltaer in material.deltaer))) materialtext.append('') - materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau)) - materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha)) + materialtext.append(', '.join('{:g}'.format(tau for tau in material.tau))) + materialtext.append(', '.join('{:g}'.format(alpha for alpha in material.alpha))) materialtext.append('') elif 'drude' in material.type: materialtext.append('') materialtext.append('') - materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau)) + materialtext.append(', '.join('{:g}'.format(tau for tau in material.tau))) materialtext.append('') - materialtext.append(', '.join('{:g}'.format(alpha) for tau in material.alpha)) + materialtext.append(', '.join('{:g}'.format(alpha for alpha in material.alpha))) else: materialtext.append('') materialtext.append('') diff --git a/gprMax/pml.py b/gprMax/pml.py index 934c6013..7c3efaaa 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -262,7 +262,7 @@ class PML(object): def build_pmls(G, pbar): """This function builds instances of the PML and calculates the initial parameters and coefficients including setting profile (based on underlying material er and mr from solid array). - + Args: G (class): Grid class instance - holds essential parameters describing the model. pbar (class): Progress bar class instance. @@ -317,6 +317,6 @@ def build_pmls(G, pbar): summr += material.mr averageer = sumer / (G.nx * G.ny) averagemr = summr / (G.nx * G.ny) - + pml.calculate_update_coeffs(averageer, averagemr, G) pbar.update() diff --git a/gprMax/sources.py b/gprMax/sources.py index d1e2220c..48608050 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -364,7 +364,7 @@ class PlaneWave(Source): """ super(Source, self).__init__() - + # Coordinates defining Huygen's surface self.xs = 0 self.xf = 0 @@ -372,54 +372,54 @@ class PlaneWave(Source): self.yf = 0 self.zs = 0 self.zf = 0 - + # Spherical coordinates defining incident unit wavevector (k) - self.theta = 0 # 0 <= theta <= 180 - self.phi = 0 # 0 <= phi <= 360 - + self.theta = 0 # 0 <= theta <= 180 + self.phi = 0 # 0 <= phi <= 360 + # Angle that incident electric field makes with k cross z - self.psi = 0 # 0 <= psi <= 360 - + self.psi = 0 # 0 <= psi <= 360 + def calculate_origin(self, G): """Calculate origin of TF/SF interface with incident wavefront.""" - + if self.theta >= 0 and self.theta <= 90: if self.phi >= 0 and self.phi <= 90: self.xcoordorigin = 0 self.ycoordorigin = 0 self.zcoordorigin = 0 - + elif self.phi > 90 and self.phi <= 180: self.xcoordorigin = G.nx self.ycoordorigin = 0 self.zcoordorigin = 0 - + elif self.phi > 180 and self.phi <= 270: self.xcoordorigin = G.nx self.ycoordorigin = G.ny self.zcoordorigin = 0 - + elif self.phi > 270 and self.phi <= 360: self.xcoordorigin = 0 self.ycoordorigin = G.ny self.zcoordorigin = 0 - + elif self.theta > 90 and self.theta <= 180: if self.phi >= 0 and self.phi <= 90: self.xcoordorigin = 0 self.ycoordorigin = 0 self.zcoordorigin = G.nz - + elif self.phi > 90 and self.phi <= 180: self.xcoordorigin = G.nx self.ycoordorigin = 0 self.zcoordorigin = G.nz - + elif self.phi > 180 and self.phi <= 270: self.xcoordorigin = G.nx self.ycoordorigin = G.ny self.zcoordorigin = G.nz - + elif self.phi > 270 and self.phi <= 360: self.xcoordorigin = 0 self.ycoordorigin = G.ny @@ -427,31 +427,20 @@ class PlaneWave(Source): def calculate_vector_components(self): """Calculate components of incident fields.""" - + self.theta = np.deg2rad(self.theta) self.phi = np.deg2rad(self.phi) self.psi = np.deg2rad(self.psi) # Components of incident unit wavevector - self.kx = np.sin(theta) * np.cos(phi) - self.ky = np.sin(theta) * np.sin(phi) - self.kz = np.cos(theta) + self.kx = np.sin(self.theta) * np.cos(self.phi) + self.ky = np.sin(self.theta) * np.sin(self.phi) + self.kz = np.cos(self.theta) # Components of incident field vectors - self.Exinc = np.cos(psi) * np.sin(phi) - np.sin(psi) * np.cos(theta) * np.cos(phi) - self.Eyinc = -np.cos(psi) * np.cos(phi) - np.sin(psi) * np.cos(theta) * np.sin(phi) - self.Ezinc = np.sin(psi) * np.sin(theta) - self.Hxinc = np.sin(psi) * np.sin(phi) + np.cos(psi) * np.cos(theta) * np.cos(phi) - self.Hyinc = -np.sin(psi) * np.cos(phi) + np.cos(psi) * np.cos(theta) * np.sin(phi) - self.Hzinc = -np.cos(psi) * np.sin(theta) - - - - - - - - - - - + self.Exinc = np.cos(self.psi) * np.sin(self.phi) - np.sin(self.psi) * np.cos(self.theta) * np.cos(self.phi) + self.Eyinc = -np.cos(self.psi) * np.cos(self.phi) - np.sin(self.psi) * np.cos(self.theta) * np.sin(self.phi) + self.Ezinc = np.sin(self.psi) * np.sin(self.theta) + self.Hxinc = np.sin(self.psi) * np.sin(self.phi) + np.cos(self.psi) * np.cos(self.theta) * np.cos(self.phi) + self.Hyinc = -np.sin(self.psi) * np.cos(self.phi) + np.cos(self.psi) * np.cos(self.theta) * np.sin(self.phi) + self.Hzinc = -np.cos(self.psi) * np.sin(self.theta)