Whitespace (W292,W293,W391) tidy.

这个提交包含在:
Craig Warren
2016-07-21 16:08:39 +01:00
父节点 cd823901d6
当前提交 e378bd66f9
共有 19 个文件被更改,包括 275 次插入317 次删除

查看文件

@@ -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'

查看文件

@@ -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
complextype = np.complex64

查看文件

@@ -22,4 +22,4 @@ class GeneralError(ValueError):
class CmdInputError(ValueError):
"""Handles errors in user specified commands. Subclasses the ValueError class."""
pass
pass

查看文件

@@ -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

查看文件

@@ -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

查看文件

@@ -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

查看文件

@@ -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

查看文件

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

查看文件

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

查看文件

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

查看文件

@@ -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()

查看文件

@@ -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, '<string>', '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()

查看文件

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

查看文件

@@ -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]
tl.Itotal[timestep] = tl.current[tl.antpos]

查看文件

@@ -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</AppendedData>\n</VTKFile>'.encode('utf-8'))
self.filehandle.close()

查看文件

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

查看文件

@@ -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 <http://www.gnu.org/licenses/>."""
gprMaxlogo = """
__ __
__ _ _ __ _ __| \/ | __ ___ __
@@ -52,10 +52,10 @@ along with gprMax. If not, see <http://www.gnu.org/licenses/>."""
\__, | .__/|_| |_| |_|\__,_/_/\_\\
|___/|_|
"""
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.')

查看文件

@@ -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
return waveform

查看文件

@@ -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']