publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-04 11:36:52 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

Formatting cleanups.

浏览代码
这个提交包含在：
Craig Warren
2017-01-26 16:23:29 +00:00
父节点 1cb71fdb2c
当前提交 c03b5efa63
共有 25 个文件被更改，包括 322 次插入 和 328 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

26
docs/source/conf.py
查看文件

@@ -225,25 +225,25 @@ htmlhelp_basename = 'gprMaxdoc'
# -- Options for LaTeX output ---------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
'papersize': 'a4paper',
# The paper size ('letterpaper' or 'a4paper').
'papersize': 'a4paper',
# The font size ('10pt', '11pt' or '12pt').
#'pointsize': '10pt',
# The font size ('10pt', '11pt' or '12pt').
#'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#'preamble': '',
# Additional stuff for the LaTeX preamble.
#'preamble': '',
# Latex figure (float) alignment
#'figure_align': 'htbp',
# Latex figure (float) alignment
#'figure_align': 'htbp',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title,
# author, documentclass [howto, manual, or own class]).
latex_documents = [
(master_doc, 'gprMax.tex', 'gprMax User Guide',
'Craig Warren and Antonis Giannopoulos', 'manual'),
(master_doc, 'gprMax.tex', 'gprMax User Guide',
'Craig Warren and Antonis Giannopoulos', 'manual'),
]
# The name of an image file (relative to this directory) to place at the top of
@@ -286,9 +286,9 @@ man_pages = [
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
(master_doc, 'gprMax', 'gprMax User Guide',
author, 'gprMax', 'One line description of project.',
'Miscellaneous'),
(master_doc, 'gprMax', 'gprMax User Guide',
author, 'gprMax', 'One line description of project.',
'Miscellaneous'),
]
# Documents to append as an appendix to all manuals.

1
gprMax/exceptions.py
查看文件

@@ -23,6 +23,7 @@ init()
sys.tracebacklimit = None
class GeneralError(ValueError):
"""Handles general errors. Subclasses the ValueError class."""

4
gprMax/gprMax.py
查看文件

@@ -110,7 +110,7 @@ def run_main(args):
# Get information about host machine
hostinfo = get_host_info()
print('\nHost: {}; {} ({} cores); {} RAM; {}'.format(hostinfo['machineID'], hostinfo['cpuID'], hostinfo['cpucores'], human_size(hostinfo['ram'], a_kilobyte_is_1024_bytes=True), hostinfo['osversion']))
# Create a separate namespace that users can access in any Python code blocks in the input file
usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'input_directory': inputdirectory}
@@ -211,7 +211,7 @@ def run_benchmark_sim(args, inputfile, usernamespace):
# Get information about host machine
hostinfo = get_host_info()
machineIDlong = '; '.join([hostinfo['machineID'], hostinfo['cpuID'], hostinfo['osversion']])
# Number of threads to test - start from max physical CPU cores and divide in half until 1
minthreads = 1
maxthreads = hostinfo['cpucores']

4
gprMax/grid.py
查看文件

@@ -196,7 +196,7 @@ def dispersion_analysis(G):
# Shift powers so that frequency with maximum power is at zero decibels
power -= np.amax(power)
# Get frequency for max power
freqmaxpower = np.where(power[1::] == np.amax(power[1::]))[0][0]
@@ -253,7 +253,7 @@ def dispersion_analysis(G):
# Physical phase velocity error (percentage)
results['deltavp'] = (((vp * c) - c) / c) * 100
# Store rounded down value of grid sampling density
results['N'] = int(np.floor(results['N']))

10
gprMax/input_cmds_multiuse.py
查看文件

@@ -86,12 +86,12 @@ def process_multicmds(multicmds, G):
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
xcoord = G.calculate_coord('x', tmp[1])
ycoord = G.calculate_coord('y', tmp[2])
zcoord = G.calculate_coord('z', tmp[3])
resistance = float(tmp[4])
check_coordinates(xcoord, ycoord, zcoord)
if xcoord < G.pmlthickness['xminus'] or xcoord > G.nx - G.pmlthickness['xplus'] or ycoord < G.pmlthickness['yminus'] or ycoord > G.ny - G.pmlthickness['yplus'] or zcoord < G.pmlthickness['zminus'] or zcoord > G.nz - G.pmlthickness['zplus']:
print(Fore.RED + "WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.' + Style.RESET_ALL)
@@ -149,7 +149,7 @@ def process_multicmds(multicmds, G):
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
xcoord = G.calculate_coord('x', tmp[1])
ycoord = G.calculate_coord('y', tmp[2])
zcoord = G.calculate_coord('z', tmp[3])
@@ -222,7 +222,7 @@ def process_multicmds(multicmds, G):
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
xcoord = G.calculate_coord('x', tmp[1])
ycoord = G.calculate_coord('y', tmp[2])
zcoord = G.calculate_coord('z', tmp[3])
@@ -283,7 +283,7 @@ def process_multicmds(multicmds, G):
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
xcoord = G.calculate_coord('x', tmp[1])
ycoord = G.calculate_coord('y', tmp[2])
zcoord = G.calculate_coord('z', tmp[3])

6
gprMax/pml.py
查看文件

@@ -166,7 +166,7 @@ class PML(object):
self.nx = xf - xs
self.ny = yf - ys
self.nz = zf - zs
if self.direction[0] == 'x':
self.d = G.dx
self.thickness = self.nx
@@ -176,11 +176,11 @@ class PML(object):
elif self.direction[0] == 'z':
self.d = G.dz
self.thickness = self.nz
self.CFS = G.cfs
if not self.CFS:
self.CFS = [CFS()]
self.initialise_field_arrays()
def initialise_field_arrays(self):

2
gprMax/sources.py
查看文件

@@ -187,7 +187,7 @@ class MagneticDipole(Source):
Hx[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dy * G.dz))
elif self.polarisation == 'y':
Hy[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dy * G.dz))
elif self.polarisation == 'z':

28
setup.py
查看文件

@@ -36,12 +36,12 @@ import sys
# Importing _version__.py before building can cause issues.
with open('gprMax/_version.py', 'r') as fd:
version = re.search(r'^__version__\s*=\s*[\'"]([^\'"]*)[\'"]',
fd.read(), re.MULTILINE).group(1)
fd.read(), re.MULTILINE).group(1)
# Parse package name from init file. Importing __init__.py / gprMax will break as gprMax depends on compiled .pyx files.
with open('gprMax/__init__.py', 'r') as fd:
packagename = re.search(r'^__name__\s*=\s*[\'"]([^\'"]*)[\'"]',
fd.read(), re.MULTILINE).group(1)
fd.read(), re.MULTILINE).group(1)
packages = [packagename, 'tests', 'tools', 'user_libs']
@@ -104,7 +104,7 @@ if 'cleanall' in sys.argv:
# Set compiler options
# Windows
if sys.platform == 'win32':
compile_args = ['/O2', '/openmp', '/w'] # No static linking as no static version of OpenMP library.
compile_args = ['/O2', '/openmp', '/w'] # No static linking as no static version of OpenMP library.
linker_args = []
extra_objects = []
# Mac OS X - needs gcc (usually via HomeBrew) because the default compiler LLVM (clang) does not support OpenMP
@@ -116,7 +116,7 @@ elif sys.platform == 'darwin':
os.environ['CC'] = gccpath[-1].split(os.sep)[-1]
else:
raise('Cannot find gcc 4.x, 5.x or 6.x in /usr/local/bin. gprMax requires gcc to be installed - easily done through the Homebrew package manager (http://brew.sh). Note: gcc with OpenMP support, i.e. --without-multilib, must be installed')
compile_args = ['-O3', '-w', '-fopenmp', '-march=native'] # Sometimes worth testing with '-fstrict-aliasing', '-fno-common'
compile_args = ['-O3', '-w', '-fopenmp', '-march=native'] # Sometimes worth testing with '-fstrict-aliasing', '-fno-common'
linker_args = ['-fopenmp']
extra_objects = []
# Linux
@@ -139,7 +139,7 @@ for file in cythonfiles:
include_dirs=[np.get_include()],
extra_compile_args=compile_args,
extra_link_args=linker_args,
extra_objects = extra_objects)
extra_objects=extra_objects)
extensions.append(extension)
# Cythonize (build .c files)
@@ -162,15 +162,15 @@ setup(name=packagename,
description='Electromagnetic Modelling Software based on the Finite-Difference Time-Domain (FDTD) method',
license='GPLv3+',
classifiers=[
'Environment :: Console',
'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)',
'Operating System :: MacOS :: MacOS X',
'Operating System :: Microsoft :: Windows :: Windows 7',
'Operating System :: POSIX :: Linux',
'Programming Language :: Cython',
'Programming Language :: Python :: 3.4',
'Topic :: Scientific/Engineering'
],
'Environment :: Console',
'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)',
'Operating System :: MacOS :: MacOS X',
'Operating System :: Microsoft :: Windows :: Windows 7',
'Operating System :: POSIX :: Linux',
'Programming Language :: Cython',
'Programming Language :: Python :: 3.4',
'Topic :: Scientific/Engineering'
],
ext_modules=extensions,
packages=packages,
include_dirs=[np.get_include()])

29
tests/analytical_solutions.py
查看文件

@@ -6,13 +6,13 @@ from gprMax.waveforms import Waveform
def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
"""Analytical solution of a z-directed Hertzian dipole in free space with a Gaussian current waveform (http://dx.doi.org/10.1016/0021-9991(83)90103-1).
Args:
iterations (int): Number of time steps.
dt (float): Time step (seconds).
dxdydz (float): Tuple of spatial resolution (metres).
rx (float): Tuple of coordinates of receiver position relative to transmitter position (metres).
Returns:
fields (float): Array contain electric and magnetic field components.
"""
@@ -22,7 +22,7 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
w.type = 'gaussiandot'
w.amp = 1
w.freq = 1e9
# Waveform integral
wint = Waveform()
wint.type = 'gaussian'
@@ -43,7 +43,7 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
dx = dxdydz[0]
dy = dxdydz[1]
dz = dxdydz[2]
# Length of Hertzian dipole
dl = dz
@@ -103,35 +103,35 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
# Calculate fields
for timestep in range(iterations):
# Calculate values for waveform, I * dl (current multiplied by dipole length) to match gprMax behaviour
fint_Ex = wint.calculate_value((timestep * dt) - tau_Ex, dt) * dl
f_Ex = w.calculate_value((timestep * dt) - tau_Ex, dt) * dl
fdot_Ex = wdot.calculate_value((timestep * dt) - tau_Ex, dt) * dl
fint_Ey = wint.calculate_value((timestep * dt) - tau_Ey, dt) * dl
f_Ey= w.calculate_value((timestep * dt) - tau_Ey, dt) * dl
f_Ey = w.calculate_value((timestep * dt) - tau_Ey, dt) * dl
fdot_Ey = wdot.calculate_value((timestep * dt) - tau_Ey, dt) * dl
fint_Ez = wint.calculate_value((timestep * dt) - tau_Ez, dt) * dl
f_Ez = w.calculate_value((timestep * dt) - tau_Ez, dt) * dl
fdot_Ez = wdot.calculate_value((timestep * dt) - tau_Ez, dt) * dl
fint_Hx = wint.calculate_value((timestep * dt) - tau_Hx, dt) * dl
f_Hx = w.calculate_value((timestep * dt) - tau_Hx, dt) * dl
fdot_Hx = wdot.calculate_value((timestep * dt) - tau_Hx, dt) * dl
fint_Hy = wint.calculate_value((timestep * dt) - tau_Hy, dt) * dl
f_Hy= w.calculate_value((timestep * dt) - tau_Hy, dt) * dl
f_Hy = w.calculate_value((timestep * dt) - tau_Hy, dt) * dl
fdot_Hy = wdot.calculate_value((timestep * dt) - tau_Hy, dt) * dl
fint_Hz = wint.calculate_value((timestep * dt) - tau_Hz, dt) * dl
f_Hz = w.calculate_value((timestep * dt) - tau_Hz, dt) * dl
fdot_Hz = wdot.calculate_value((timestep * dt) - tau_Hz, dt) * dl
# Ex
fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * e0 * Er_x**5)) * (3 * (fint_Ex + (tau_Ex * f_Ex)) + (tau_Ex**2 * fdot_Ex))
# Ey
try:
tmp = Ey_y / Ey_x
@@ -156,4 +156,3 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
fields[timestep, 5] = 0
return fields

4
tests/benchmarking/plot_benchmark.py
查看文件

@@ -23,7 +23,7 @@ baseresult = np.load(args.baseresult)
# Get machine/CPU/OS details
try:
machineIDlong = str(baseresult['machineID'])
#machineIDlong = 'Dell PowerEdge R630; Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz; Linux (3.10.0-327.18.2.el7.x86_64)' # Use to manually describe machine
# machineIDlong = 'Dell PowerEdge R630; Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz; Linux (3.10.0-327.18.2.el7.x86_64)' # Use to manually describe machine
machineID = machineIDlong.split(';')[0]
except KeyError:
hostinfo = get_host_info()
@@ -112,5 +112,3 @@ ax.set_ylim(bottom=1, top=ax.get_ylim()[1] * 1.1)
fig.savefig(os.path.join(os.path.dirname(args.baseresult), machineID.replace(' ', '_') + '.png'), dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
plt.show()

10
tests/test_experimental.py
查看文件

@@ -66,15 +66,15 @@ modelmax = np.where(np.abs(model) == 1)[0][0]
# Real results
with open(args.realfile, 'r') as f:
real = np.loadtxt(f)
real[:,1] = real[:,1] / np.amax(np.abs(real[:,1]))
realmax = np.where(np.abs(real[:,1]) == 1)[0][0]
real[:, 1] = real[:, 1] / np.amax(np.abs(real[:, 1]))
realmax = np.where(np.abs(real[:, 1]) == 1)[0][0]
difftime = - (timemodel[modelmax] - real[realmax,0])
difftime = - (timemodel[modelmax] - real[realmax, 0])
# Plot modelled and real data
fig, ax = plt.subplots(num=args.modelfile + ' versus ' + args.realfile, figsize=(20, 10), facecolor='w', edgecolor='w')
ax.plot(timemodel + difftime, model, 'r', lw=2, label='Model')
ax.plot(real[:,0], real[:,1], 'r', ls='--', lw=2, label='Experiment')
ax.plot(real[:, 0], real[:, 1], 'r', ls='--', lw=2, label='Experiment')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_xlim([0, timemodel[-1]])
@@ -85,6 +85,6 @@ ax.grid()
# Save a PDF/PNG of the figure
savename = os.path.abspath(os.path.dirname(args.modelfile)) + os.sep + os.path.splitext(os.path.split(args.modelfile)[1])[0] + '_vs_' + os.path.splitext(os.path.split(args.realfile)[1])[0]
#fig.savefig(savename + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
#fig.savefig((savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
# fig.savefig((savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
plt.show()

77
tests/test_models_basic.py
查看文件

@@ -36,11 +36,10 @@ from gprMax.exceptions import GeneralError
from tests.analytical_solutions import hertzian_dipole_fs
"""Compare field outputs
Usage:
cd gprMax
python -m tests.test_basic_suite
"""
basepath = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'models_basic')
@@ -57,21 +56,21 @@ path = '/rxs/rx1/'
starttime = perf_counter()
for i, model in enumerate(testmodels):
testresults[model] = {}
# Run model
api(os.path.join(basepath, model + os.path.sep + model + '.in'))
# Special case for analytical comparison
if model == 'hertzian_dipole_fs_analytical':
# Get output for model file
filetest = h5py.File(os.path.join(basepath, model + os.path.sep + model + '.out'),'r')
filetest = h5py.File(os.path.join(basepath, model + os.path.sep + model + '.out'), 'r')
testresults[model]['Test version'] = filetest.attrs['gprMax']
# Get available field output component names
outputstest = list(filetest[path].keys())
# Arrays for storing time
floattype = filetest[path + outputstest[0]].dtype
timetest = np.zeros((filetest.attrs['Iterations']), dtype=floattype)
@@ -82,12 +81,12 @@ for i, model in enumerate(testmodels):
datatest = np.zeros((filetest.attrs['Iterations'], len(outputstest)), dtype=floattype)
print(datatest.shape)
for ID, name in enumerate(outputstest):
datatest[:,ID] = filetest[path + str(name)][:]
datatest[:, ID] = filetest[path + str(name)][:]
# Tx/Rx position to feed to analytical solution
rxpos = filetest[path].attrs['Position']
txpos = filetest['/srcs/src1/'].attrs['Position']
rxposrelative = ((rxpos[0] - txpos[0]), (rxpos[1] - txpos[1]), (rxpos[2] - txpos[2]))
rxposrelative = ((rxpos[0] - txpos[0]), (rxpos[1] - txpos[1]), (rxpos[2] - txpos[2]))
# Analytical solution of a dipole in free space
dataref = hertzian_dipole_fs(filetest.attrs['Iterations'], filetest.attrs['dt'], filetest.attrs['dx, dy, dz'], rxposrelative)
@@ -97,15 +96,15 @@ for i, model in enumerate(testmodels):
# Diffs
datadiffs = np.zeros(datatest.shape, dtype=floattype)
for i in range(len(outputstest)):
max = np.amax(np.abs(dataref[:,i]))
max = np.amax(np.abs(dataref[:, i]))
try:
datadiffs[:,i] = ((np.abs(dataref[:,i] - datatest[:,i])) / max) * 100
datadiffs[:, i] = ((np.abs(dataref[:, i] - datatest[:, i])) / max) * 100
except FloatingPointError:
print('FloatingPointError')
datadiffs[:,i] = 0
datadiffs[:, i] = 0
# Register test passed
threshold = 2 # Percent
threshold = 2 # Percent
if np.amax(np.amax(datadiffs)) < 2:
testresults[model]['Pass'] = True
else:
@@ -114,17 +113,17 @@ for i, model in enumerate(testmodels):
else:
# Get output for model and reference files
fileref = h5py.File(os.path.join(basepath, model + os.path.sep + model + '_ref.out'),'r')
filetest = h5py.File(os.path.join(basepath, model + os.path.sep + model + '.out'),'r')
fileref = h5py.File(os.path.join(basepath, model + os.path.sep + model + '_ref.out'), 'r')
filetest = h5py.File(os.path.join(basepath, model + os.path.sep + model + '.out'), 'r')
testresults[model]['Ref version'] = fileref.attrs['gprMax']
testresults[model]['Test version'] = filetest.attrs['gprMax']
# Get available field output component names
outputsref = list(fileref[path].keys())
outputstest = list(filetest[path].keys())
if outputsref != outputstest:
raise GeneralError('Field output components do not match reference solution')
# Check that type of float used to store fields matches
if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype:
raise GeneralError('Type of floating point number does not match reference solution')
@@ -147,8 +146,8 @@ for i, model in enumerate(testmodels):
dataref = np.zeros((fileref.attrs['Iterations'], len(outputsref)), dtype=floattype)
datatest = np.zeros((filetest.attrs['Iterations'], len(outputstest)), dtype=floattype)
for ID, name in enumerate(outputsref):
dataref[:,ID] = fileref[path + str(name)][:]
datatest[:,ID] = filetest[path + str(name)][:]
dataref[:, ID] = fileref[path + str(name)][:]
datatest[:, ID] = filetest[path + str(name)][:]
fileref.close()
filetest.close()
@@ -156,12 +155,12 @@ for i, model in enumerate(testmodels):
# Diffs
datadiffs = np.zeros(datatest.shape, dtype=floattype)
for i in range(len(outputstest)):
max = np.nanmax(np.abs(dataref[:,i]))
max = np.nanmax(np.abs(dataref[:, i]))
try:
datadiffs[:,i] = ((np.abs(dataref[:,i] - datatest[:,i])) / max) * 100
datadiffs[:, i] = ((np.abs(dataref[:, i] - datatest[:, i])) / max) * 100
except FloatingPointError:
print('FloatingPointError')
datadiffs[:,i] = 0
datadiffs[:, i] = 0
# Register test passed
if not np.any(datadiffs):
@@ -172,18 +171,18 @@ for i, model in enumerate(testmodels):
# Plot datasets
fig1, ((ex1, hx1), (ey1, hy1), (ez1, hz1)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num=model + '.in', figsize=(20, 10), facecolor='w', edgecolor='w')
ex1.plot(timetest, datatest[:,0],'r', lw=2, label=model)
ex1.plot(timeref, dataref[:,0],'g', lw=2, ls='--', label=model + '(Ref)')
ey1.plot(timetest, datatest[:,1],'r', lw=2, label=model)
ey1.plot(timeref, dataref[:,1],'g', lw=2, ls='--', label=model + '(Ref)')
ez1.plot(timetest, datatest[:,2],'r', lw=2, label=model)
ez1.plot(timeref, dataref[:,2],'g', lw=2, ls='--', label=model + '(Ref)')
hx1.plot(timetest, datatest[:,3],'r', lw=2, label=model)
hx1.plot(timeref, dataref[:,3],'g', lw=2, ls='--', label=model + '(Ref)')
hy1.plot(timetest, datatest[:,4],'r', lw=2, label=model)
hy1.plot(timeref, dataref[:,4],'g', lw=2, ls='--', label=model + '(Ref)')
hz1.plot(timetest, datatest[:,5],'r', lw=2, label=model)
hz1.plot(timeref, dataref[:,5],'g', lw=2, ls='--', label=model + '(Ref)')
ex1.plot(timetest, datatest[:, 0], 'r', lw=2, label=model)
ex1.plot(timeref, dataref[:, 0], 'g', lw=2, ls='--', label=model + '(Ref)')
ey1.plot(timetest, datatest[:, 1], 'r', lw=2, label=model)
ey1.plot(timeref, dataref[:, 1], 'g', lw=2, ls='--', label=model + '(Ref)')
ez1.plot(timetest, datatest[:, 2], 'r', lw=2, label=model)
ez1.plot(timeref, dataref[:, 2], 'g', lw=2, ls='--', label=model + '(Ref)')
hx1.plot(timetest, datatest[:, 3], 'r', lw=2, label=model)
hx1.plot(timeref, dataref[:, 3], 'g', lw=2, ls='--', label=model + '(Ref)')
hy1.plot(timetest, datatest[:, 4], 'r', lw=2, label=model)
hy1.plot(timeref, dataref[:, 4], 'g', lw=2, ls='--', label=model + '(Ref)')
hz1.plot(timetest, datatest[:, 5], 'r', lw=2, label=model)
hz1.plot(timeref, dataref[:, 5], 'g', lw=2, ls='--', label=model + '(Ref)')
ylabels = ['$E_x$, field strength [V/m]', '$H_x$, field strength [A/m]', '$E_y$, field strength [V/m]', '$H_y$, field strength [A/m]', '$E_z$, field strength [V/m]', '$H_z$, field strength [A/m]']
for i, ax in enumerate(fig1.axes):
ax.set_ylabel(ylabels[i])
@@ -193,12 +192,12 @@ for i, model in enumerate(testmodels):
# Plot diffs
fig2, ((ex2, hx2), (ey2, hy2), (ez2, hz2)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num='Diffs: ' + model + '.in', figsize=(20, 10), facecolor='w', edgecolor='w')
ex2.plot(timeref, datadiffs[:,0],'r', lw=2, label='Ex')
ey2.plot(timeref, datadiffs[:,1],'r', lw=2, label='Ey')
ez2.plot(timeref, datadiffs[:,2],'r', lw=2, label='Ez')
hx2.plot(timeref, datadiffs[:,3],'r', lw=2, label='Hx')
hy2.plot(timeref, datadiffs[:,4],'r', lw=2, label='Hy')
hz2.plot(timeref, datadiffs[:,5],'r', lw=2, label='Hz')
ex2.plot(timeref, datadiffs[:, 0], 'r', lw=2, label='Ex')
ey2.plot(timeref, datadiffs[:, 1], 'r', lw=2, label='Ey')
ez2.plot(timeref, datadiffs[:, 2], 'r', lw=2, label='Ez')
hx2.plot(timeref, datadiffs[:, 3], 'r', lw=2, label='Hx')
hy2.plot(timeref, datadiffs[:, 4], 'r', lw=2, label='Hy')
hz2.plot(timeref, datadiffs[:, 5], 'r', lw=2, label='Hz')
ylabels = ['$E_x$, difference [%]', '$H_x$, difference [%]', '$E_y$, difference [%]', '$H_y$, difference [%]', '$E_z$, difference [%]', '$H_z$, difference [%]']
for i, ax in enumerate(fig2.axes):
ax.set_ylabel(ylabels[i])

20
tools/Paraview macros/gprMax_info.py
查看文件

@@ -56,7 +56,7 @@ with open(model.FileName[0], 'r') as f:
if materials:
# Get range of data
mat_datarange = model.CellData.GetArray('Material').GetRange()
# Create threshold for materials (name and numeric value)
for x in range(0, int(mat_datarange[1]) + 1):
for y in range(len(materials)):
@@ -64,9 +64,9 @@ if materials:
threshold = Threshold(Input=model)
threshold.Scalars = 'Material'
threshold.ThresholdRange = [materials[y][0], materials[y][0]]
RenameSource(materials[y][1], threshold)
if materials[y][0] != 1:
# Show data in view
thresholddisplay = Show(threshold, renderview)
@@ -75,7 +75,7 @@ if materials:
if srcs_pml:
# Get ranges of data
srcs_pml_datarange = model.CellData.GetArray('Sources_PML').GetRange()
# Create threshold for sources/pml (name and numeric value)
for x in range(1, int(srcs_pml_datarange[1]) + 1):
for y in range(len(srcs_pml)):
@@ -83,20 +83,20 @@ if srcs_pml:
threshold = Threshold(Input=model)
threshold.Scalars = 'Sources_PML'
threshold.ThresholdRange = [srcs_pml[y][0], srcs_pml[y][0]]
RenameSource(srcs_pml[y][1], threshold)
# Show data in view
thresholddisplay = Show(threshold, renderview)
thresholddisplay.ColorArrayName = 'Sources_PML'
if srcs_pml[y][0] == 1:
thresholddisplay.Opacity = 0.5
if rxs:
# Get ranges of data
rxs_datarange = model.CellData.GetArray('Receivers').GetRange()
# Create threshold for sources/pml (name and numeric value)
for x in range(1, int(rxs_datarange[1]) + 1):
for y in range(len(rxs)):
@@ -104,9 +104,9 @@ if rxs:
threshold = Threshold(Input=model)
threshold.Scalars = 'Receivers'
threshold.ThresholdRange = [rxs[y][0], rxs[y][0]]
RenameSource(rxs[y][1], threshold)
# Show data in view
thresholddisplay = Show(threshold, renderview)
thresholddisplay.ColorArrayName = 'Receivers'

20
tools/inputfile_old2new.py
查看文件

@@ -51,7 +51,7 @@ hertziandipoles = []
lindex = 0
while(lindex < len(inputlines)):
if inputlines[lindex].startswith('#') and not inputlines[lindex].startswith('##'):
cmd = inputlines[lindex].split(':')
cmdname = cmd[0].lower()
@@ -66,7 +66,7 @@ while(lindex < len(inputlines)):
inputlines.insert(lindex, replacement)
lindex += 1
dx_dy_dz = (float(params[0]), float(params[1]), float(params[1]))
elif cmdname == '#dx_dy_dz':
dx_dy_dz = (float(params[0]), float(params[1]), float(params[2]))
lindex += 1
@@ -75,16 +75,16 @@ while(lindex < len(inputlines)):
lindex += 1
else:
lindex +=1
lindex += 1
lindex = 0
while(lindex < len(inputlines)):
if inputlines[lindex].startswith('#') and not inputlines[lindex].startswith('##'):
cmd = inputlines[lindex].split(':')
cmdname = cmd[0].lower()
params = cmd[1].split()
if cmdname == '#domain':
if model2D:
# Syntax of old command: #domain: x y
@@ -105,7 +105,7 @@ while(lindex < len(inputlines)):
else:
timewindow = int(params[0])
lindex += 1
elif cmdname == '#num_of_procs':
# Syntax of old command: #num_of_procs: nthreads
replacement = '#num_threads: {}'.format(params[0])
@@ -222,7 +222,7 @@ while(lindex < len(inputlines)):
replacement = '#cylinder: {} {} {} {} {} {} {} {}'.format(params[1], params[3], params[4], params[2], params[3], params[4], params[5], params[6])
elif params[0] == 'y':
replacement = '#cylinder: {} {} {} {} {} {} {} {}'.format(params[3], params[1], params[4], params[3], params[2], params[4], params[5], params[6])
elif params[0] =='z':
elif params[0] == 'z':
replacement = '#cylinder: {} {} {} {} {} {} {} {}'.format(params[3], params[4], params[1], params[3], params[4], params[2], params[5], params[6])
print("Command '{}', replaced with '{}'".format(inputlines[lindex], replacement))
@@ -237,7 +237,7 @@ while(lindex < len(inputlines)):
inputlines.pop(lindex)
inputlines.insert(lindex, replacement)
lindex += 1
elif cmdname == '#cylindrical_segment':
print("Command '{}' has been removed as it is no longer supported. You can create a cylindrical segment by using a #box to cut through a #cylinder.".format(inputlines[lindex]))
inputlines.pop(lindex)
@@ -312,7 +312,7 @@ while(lindex < len(inputlines)):
lindex += 1
else:
lindex +=1
lindex += 1
# Convert separate #line_source and associated #tx to #waveform and #hertzian_dipole
for source in linesources:
@@ -370,7 +370,7 @@ for source in voltagesources:
voltagesource = '#voltage_source: {} {} {} {} {} {} {} {}'.format(voltagesourcetx[1], voltagesourcetx[2], voltagesourcetx[3], voltagesourcetx[4], params[4], voltagesourcetx[5], voltagesourcetx[6], voltagesourcetx[7])
else:
voltagesource = '#voltage_source: {} {} {} {} {} {}'.format(voltagesourcetx[1], voltagesourcetx[2], voltagesourcetx[3], voltagesourcetx[4], params[4], voltagesourcetx[5])
print("Commands '{}' and '{}', replaced with '{}' and '{}'".format(source, tx, waveform, voltagesource))
inputlines.remove(source)
inputlines.remove(tx)

5
tools/outputfiles_merge.py
查看文件

@@ -44,7 +44,7 @@ fout = h5py.File(outputfile, 'w')
for model in range(modelruns):
fin = h5py.File(basefilename + str(model + 1) + '.out', 'r')
nrx = fin.attrs['nrx']
# Write properties for merged file on first iteration
if model == 0:
fout.attrs['Iterations'] = fin.attrs['Iterations']
@@ -63,7 +63,7 @@ for model in range(modelruns):
availableoutputs = list(fin[path].keys())
# For all receiver outputs
for output in availableoutputs:
fout[path + '/' + output][:,model] = fin[path + '/' + output][:]
fout[path + '/' + output][:, model] = fin[path + '/' + output][:]
fin.close()
@@ -74,4 +74,3 @@ if not check or check == 'y':
for model in range(modelruns):
file = basefilename + str(model + 1) + '.out'
os.remove(file)

49
tools/plot_Ascan.py
查看文件

@@ -31,16 +31,16 @@ from gprMax.receivers import Rx
def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
"""Plots electric and magnetic fields and currents from all receiver points in the given output file. Each receiver point is plotted in a new figure window.
Args:
filename (string): Filename (including path) of output file.
outputs (list): List of field/current components to plot.
fft (boolean): Plot FFT switch.
Returns:
plt (object): matplotlib plot object.
"""
# Open output file and read some attributes
f = h5py.File(filename, 'r')
nrx = f.attrs['nrx']
@@ -62,10 +62,10 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
for rx in range(1, nrx + 1):
path = '/rxs/rx' + str(rx) + '/'
availableoutputs = list(f[path].keys())
# If only a single output is required, create one subplot
if len(outputs) == 1:
# Check for polarity of output and if requested output is in file
if outputs[0][-1] == '-':
polarity = -1
@@ -75,10 +75,10 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
polarity = 1
outputtext = outputs[0]
output = outputs[0]
if output not in availableoutputs:
raise CmdInputError('{} output requested to plot, but the available output for receiver 1 is {}'.format(output, ', '.join(availableoutputs)))
outputdata = f[path + output][:] * polarity
# Plotting if FFT required
@@ -113,7 +113,7 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Power [dB]')
ax2.grid()
# Change colours and labels for magnetic field components or currents
if 'H' in outputs[0]:
plt.setp(line1, color='g')
@@ -127,17 +127,17 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
# plt.setp(ax1, ylabel=outputtext + ' current [A]')
# plt.setp(stemlines, 'color', 'b')
# plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
plt.show()
# Plotting if no FFT required
else:
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]', ylabel=outputtext + ' field strength [V/m]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
line = ax.plot(time, outputdata,'r', lw=2, label=outputtext)
line = ax.plot(time, outputdata, 'r', lw=2, label=outputtext)
ax.set_xlim([0, np.amax(time)])
#ax.set_ylim([-15, 20])
ax.grid()
if 'H' in output:
plt.setp(line, color='g')
plt.setp(ax, ylabel=outputtext + ', field strength [A/m]')
@@ -150,7 +150,7 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
gs = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3)
for output in outputs:
# Check for polarity of output and if requested output is in file
if output[-1] == 'm':
polarity = -1
@@ -159,41 +159,41 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
else:
polarity = 1
outputtext = output
# Check if requested output is in file
if output not in availableoutputs:
raise CmdInputError('Output(s) requested to plot: {}, but available output(s) for receiver {} in the file: {}'.format(', '.join(outputs), rx, ', '.join(availableoutputs)))
outputdata = f[path + output][:] * polarity
if output == 'Ex':
ax = plt.subplot(gs[0, 0])
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
ax.plot(time, outputdata, 'r', lw=2, label=outputtext)
ax.set_ylabel(outputtext + ', field strength [V/m]')
#ax.set_ylim([-15, 20])
elif output == 'Ey':
ax = plt.subplot(gs[1, 0])
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
ax.plot(time, outputdata, 'r', lw=2, label=outputtext)
ax.set_ylabel(outputtext + ', field strength [V/m]')
#ax.set_ylim([-15, 20])
elif output == 'Ez':
ax = plt.subplot(gs[2, 0])
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
ax.plot(time, outputdata, 'r', lw=2, label=outputtext)
ax.set_ylabel(outputtext + ', field strength [V/m]')
#ax.set_ylim([-15, 20])
elif output == 'Hx':
ax = plt.subplot(gs[0, 1])
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
ax.plot(time, outputdata, 'g', lw=2, label=outputtext)
ax.set_ylabel(outputtext + ', field strength [A/m]')
#ax.set_ylim([-0.03, 0.03])
elif output == 'Hy':
ax = plt.subplot(gs[1, 1])
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
ax.plot(time, outputdata, 'g', lw=2, label=outputtext)
ax.set_ylabel(outputtext + ', field strength [A/m]')
#ax.set_ylim([-0.03, 0.03])
elif output == 'Hz':
ax = plt.subplot(gs[2, 1])
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
ax.plot(time, outputdata, 'g', lw=2, label=outputtext)
ax.set_ylabel(outputtext + ', field strength [A/m]')
#ax.set_ylim([-0.03, 0.03])
# elif output == 'Ix':
@@ -228,6 +228,5 @@ if __name__ == "__main__":
parser.add_argument('-fft', action='store_true', help='plot FFT (single output must be specified)', default=False)
args = parser.parse_args()
plt = mpl_plot(args.outputfile, args.outputs, fft=args.fft)
plt.show()
plthandle = mpl_plot(args.outputfile, args.outputs, fft=args.fft)
plthandle.show()

15
tools/plot_Bscan.py
查看文件

@@ -29,12 +29,12 @@ from gprMax.exceptions import CmdInputError
def get_output_data(filename, rxnumber, rxcomponent):
"""Gets B-scan output data from a model.
Args:
filename (string): Filename (including path) of output file.
rxnumber (int): Receiver output number.
rxcomponent (str): Receiver output field/current component.
Returns:
outputdata (array): Array of A-scans, i.e. B-scan data.
dt (float): Temporal resolution of the model.
@@ -68,19 +68,19 @@ def get_output_data(filename, rxnumber, rxcomponent):
def mpl_plot(outputdata, dt, rxnumber, rxcomponent):
"""Creates a plot (with matplotlib) of the B-scan.
Args:
outputdata (array): Array of A-scans, i.e. B-scan data.
dt (float): Temporal resolution of the model.
rxnumber (int): Receiver output number.
rxcomponent (str): Receiver output field/current component.
Returns:
plt (object): matplotlib plot object.
"""
fig = plt.figure(num='rx' + str(rxnumber), figsize=(20, 10), facecolor='w', edgecolor='w')
plt.imshow(outputdata, extent=[0, outputdata.shape[1], outputdata.shape[0]*dt, 0], interpolation='nearest', aspect='auto', cmap='seismic', vmin=-np.amax(np.abs(outputdata)), vmax=np.amax(np.abs(outputdata)))
plt.imshow(outputdata, extent=[0, outputdata.shape[1], outputdata.shape[0] * dt, 0], interpolation='nearest', aspect='auto', cmap='seismic', vmin=-np.amax(np.abs(outputdata)), vmax=np.amax(np.abs(outputdata)))
plt.xlabel('Trace number')
plt.ylabel('Time [s]')
plt.grid()
@@ -118,7 +118,6 @@ if __name__ == "__main__":
for rx in range(1, nrx + 1):
outputdata, dt = get_output_data(args.outputfile, rx, args.rx_component)
plt = mpl_plot(outputdata, dt, rx, args.rx_component)
plt.show()
plthandle = mpl_plot(outputdata, dt, rx, args.rx_component)
plthandle.show()

59
tools/plot_antenna_params.py
查看文件

@@ -30,14 +30,14 @@ from gprMax.exceptions import CmdInputError
def calculate_antenna_params(filename, tltxnumber=1, tlrxnumber=None, rxnumber=None, rxcomponent=None):
"""Calculates antenna parameters - incident, reflected and total volatges and currents; s11, (s21) and input impedance.
Args:
filename (string): Filename (including path) of output file.
tltxnumber (int): Transmitter antenna - transmission line number
tlrxnumber (int): Receiver antenna - transmission line number
rxnumber (int): Receiver antenna - output number
rxcomponent (str): Receiver antenna - output electric field component
Returns:
antennaparams (dict): Antenna parameters.
"""
@@ -65,8 +65,8 @@ def calculate_antenna_params(filename, tltxnumber=1, tlrxnumber=None, rxnumber=N
Iinc = f[tltxpath + 'Iinc'][:]
# Total (incident + reflected) voltages/currents
Vtotal = f[tltxpath +'Vtotal'][:]
Itotal = f[tltxpath +'Itotal'][:]
Vtotal = f[tltxpath + 'Vtotal'][:]
Itotal = f[tltxpath + 'Itotal'][:]
# Reflected voltages/currents
Vref = Vtotal - Vinc
@@ -75,17 +75,17 @@ def calculate_antenna_params(filename, tltxnumber=1, tlrxnumber=None, rxnumber=N
# If a receiver antenna is used (with a transmission line or receiver), get received voltage for s21
if tlrxnumber:
tlrxpath = '/tls/tl' + str(tlrxnumber) + '/'
Vrec = f[tlrxpath +'Vtotal'][:]
Vrec = f[tlrxpath + 'Vtotal'][:]
elif rxnumber:
rxpath = '/rxs/rx' + str(rxnumber) + '/'
availableoutputs = list(f[rxpath].keys())
if rxcomponent not in availableoutputs:
raise CmdInputError('{} output requested, but the available output for receiver {} is {}'.format(rxcomponent, rxnumber, ', '.join(availableoutputs)))
rxpath += rxcomponent
# Received voltage
if rxcomponent == 'Ex':
Vrec = f[rxpath][:] * -1 * dxdydz[0]
@@ -120,12 +120,12 @@ def calculate_antenna_params(filename, tltxnumber=1, tlrxnumber=None, rxnumber=N
Vtotalp = 20 * np.log10(np.abs((np.fft.fft(Vtotal) * delaycorrection)))
Itotalp = 20 * np.log10(np.abs(np.fft.fft(Itotal)))
s11 = 20 * np.log10(s11)
# Create dictionary of antenna parameters
antennaparams = {'time': time, 'freqs': freqs, 'Vinc': Vinc, 'Vincp': Vincp, 'Iinc': Iinc, 'Iincp': Iincp,
'Vref': Vref, 'Vrefp': Vrefp, 'Iref': Iref, 'Irefp': Irefp,
'Vtotal': Vtotal, 'Vtotalp': Vtotalp, 'Itotal': Itotal, 'Itotalp': Itotalp,
's11': s11, 'zin': zin, 'yin': yin}
'Vref': Vref, 'Vrefp': Vrefp, 'Iref': Iref, 'Irefp': Irefp,
'Vtotal': Vtotal, 'Vtotalp': Vtotalp, 'Itotal': Itotal, 'Itotalp': Itotalp,
's11': s11, 'zin': zin, 'yin': yin}
if tlrxnumber or rxnumber:
s21 = 20 * np.log10(s21)
antennaparams['s21'] = s21
@@ -135,7 +135,7 @@ def calculate_antenna_params(filename, tltxnumber=1, tlrxnumber=None, rxnumber=N
def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref, Irefp, Vtotal, Vtotalp, Itotal, Itotalp, s11, zin, yin, s21=None):
"""Plots antenna parameters - incident, reflected and total volatges and currents; s11, (s21) and input impedance.
Args:
filename (string): Filename (including path) of output file.
time (array): Simulation time.
@@ -145,11 +145,11 @@ def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref,
Vtotal, Vtotalp, Itotal, Itotalp (array): Time and frequency domain representations of total voltage and current.
s11, s21 (array): s11 and, optionally, s21 parameters.
zin, yin (array): Input impedance and input admittance parameters.
Returns:
plt (object): matplotlib plot object.
"""
# Set plotting range
pltrangemin = 1
# To a certain drop from maximum power
@@ -253,16 +253,16 @@ def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref,
ax.set_ylabel('Power [dB]')
ax.grid()
## Plot reflected (reflected) voltage
# Plot reflected (reflected) voltage
#ax = plt.subplot(gs1[4, 0])
#ax.plot(time, Vref, 'r', lw=2, label='Vref')
#ax.set_title('Reflected voltage')
#ax.set_xlabel('Time [s]')
#ax.set_ylabel('Voltage [V]')
#ax.set_xlim([0, np.amax(time)])
#ax.grid()
# ax.grid()
#
## Plot frequency spectra of reflected voltage
# Plot frequency spectra of reflected voltage
#ax = plt.subplot(gs1[4, 1])
#markerline, stemlines, baseline = ax.stem(freqs[pltrange], Vrefp[pltrange], '-.')
#plt.setp(baseline, 'linewidth', 0)
@@ -272,18 +272,18 @@ def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref,
#ax.set_title('Reflected voltage')
#ax.set_xlabel('Frequency [Hz]')
#ax.set_ylabel('Power [dB]')
#ax.grid()
# ax.grid()
#
## Plot reflected (reflected) current
# Plot reflected (reflected) current
#ax = plt.subplot(gs1[5, 0])
#ax.plot(time, Iref, 'b', lw=2, label='Iref')
#ax.set_title('Reflected current')
#ax.set_xlabel('Time [s]')
#ax.set_ylabel('Current [A]')
#ax.set_xlim([0, np.amax(time)])
#ax.grid()
# ax.grid()
#
## Plot frequency spectra of reflected current
# Plot frequency spectra of reflected current
#ax = plt.subplot(gs1[5, 1])
#markerline, stemlines, baseline = ax.stem(freqs[pltrange], Irefp[pltrange], '-.')
#plt.setp(baseline, 'linewidth', 0)
@@ -293,7 +293,7 @@ def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref,
#ax.set_title('Reflected current')
#ax.set_xlabel('Frequency [Hz]')
#ax.set_ylabel('Power [dB]')
#ax.grid()
# ax.grid()
# Figure 2
# Plot frequency spectra of s11
@@ -356,7 +356,7 @@ def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref,
#ax.set_ylim([-300, 300])
ax.grid()
## Plot input admittance (magnitude)
# Plot input admittance (magnitude)
#ax = plt.subplot(gs2[2, 0])
#markerline, stemlines, baseline = ax.stem(freqs[pltrange], np.abs(yin[pltrange]), '-.')
#plt.setp(baseline, 'linewidth', 0)
@@ -368,9 +368,9 @@ def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref,
#ax.set_ylabel('Admittance [Siemens]')
##ax.set_xlim([0.88e9, 1.02e9])
##ax.set_ylim([0, 0.035])
#ax.grid()
# ax.grid()
#
## Plot input admittance (phase)
# Plot input admittance (phase)
#ax = plt.subplot(gs2[2, 1])
#markerline, stemlines, baseline = ax.stem(freqs[pltrange], np.angle(yin[pltrange], deg=True), '-.')
#plt.setp(baseline, 'linewidth', 0)
@@ -382,7 +382,7 @@ def mpl_plot(filename, time, freqs, Vinc, Vincp, Iinc, Iincp, Vref, Vrefp, Iref,
#ax.set_ylabel('Phase [degrees]')
##ax.set_xlim([0.88e9, 1.02e9])
##ax.set_ylim([-40, 100])
#ax.grid()
# ax.grid()
# Save a PDF/PNG of the figure
#fig1.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_tl_params.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
@@ -405,6 +405,5 @@ if __name__ == "__main__":
args = parser.parse_args()
antennaparams = calculate_antenna_params(args.outputfile, args.tltx_num, args.tlrx_num, args.rx_num, args.rx_component)
plt = mpl_plot(args.outputfile, **antennaparams)
plt.show()
plthandle = mpl_plot(args.outputfile, **antennaparams)
plthandle.show()

27
tools/plot_source_wave.py
查看文件

@@ -32,19 +32,19 @@ np.seterr(divide='ignore')
def check_timewindow(timewindow, dt):
"""Checks and sets time window and number of iterations.
Args:
timewindow (float): Time window.
dt (float): Time discretisation.
Returns:
timewindow (float): Time window.
iterations (int): Number of interations.
"""
# Time window could be a string, float or int, so convert to string then check
timewindow = str(timewindow)
try:
timewindow = int(timewindow)
iterations = timewindow
@@ -62,18 +62,18 @@ def check_timewindow(timewindow, dt):
def mpl_plot(w, timewindow, dt, iterations, fft=False):
"""Plots waveform and prints useful information about its properties.
Args:
w (class): Waveform class instance.
timewindow (float): Time window.
dt (float): Time discretisation.
iterations (int): Number of iterations.
fft (boolean): Plot FFT switch.
Returns:
plt (object): matplotlib plot object.
"""
time = np.linspace(0, 1, iterations)
time *= (iterations * dt)
waveform = np.zeros(len(time))
@@ -104,7 +104,7 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False):
# Calculate pulse width for gaussian
if w.type == 'gaussian':
powerdrop = -3 #dB
powerdrop = -3 # dB
start = np.where((10 * np.log10(waveform / np.amax(waveform))) > powerdrop)[0][0]
stop = np.where((10 * np.log10(waveform[start:] / np.amax(waveform))) < powerdrop)[0][0] + start
print('Pulse width at {:d}dB, i.e. FWHM: {:g} s'.format(powerdrop, time[stop] - time[start]))
@@ -130,7 +130,7 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False):
pltrange = np.s_[0:pltrange]
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num=w.type, figsize=(20, 10), facecolor='w', edgecolor='w')
# Plot waveform
ax1.plot(time, waveform, 'r', lw=2)
ax1.set_xlabel('Time [s]')
@@ -153,7 +153,7 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False):
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Amplitude')
[ax.grid() for ax in fig.axes] # Turn on grid
[ax.grid() for ax in fig.axes] # Turn on grid
# Save a PDF/PNG of the figure
#fig.savefig(os.path.dirname(os.path.abspath(__file__)) + os.sep + w.type + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
@@ -163,7 +163,7 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False):
if __name__ == "__main__":
# Parse command line arguments
parser = argparse.ArgumentParser(description='Plot built-in waveforms that can be used for sources.', usage='cd gprMax; python -m tools.plot_source_wave type amp freq timewindow dt')
parser.add_argument('type', help='type of waveform', choices=Waveform.types)
@@ -187,6 +187,5 @@ if __name__ == "__main__":
w.freq = args.freq
timewindow, iterations = check_timewindow(args.timewindow, args.dt)
plt = mpl_plot(w, timewindow, args.dt, iterations, args.fft)
plt.show()
plthandle = mpl_plot(w, timewindow, args.dt, iterations, args.fft)
plthandle.show()

12
user_libs/AustinManWoman/head_only_hdf5.py
查看文件

@@ -1,4 +1,6 @@
import argparse, os
import argparse
import os
import h5py
# Parse command line arguments
@@ -6,18 +8,18 @@ parser = argparse.ArgumentParser(description='Writes a HDF5 file of AustinMan or
parser.add_argument('filename', help='name and path to (HDF5) file containing AustinMan or AustinWoman model')
args = parser.parse_args()
# Read full body HDF5 file
# Read full body HDF5 file
f = h5py.File(args.filename, 'r')
dx_dy_dz = f.attrs['dx, dy, dz']
data = f['/data'][:,:,:]
data = f['/data'][:, :, :]
# Define head as last 1/8 of total body height
nzhead = 7 * int(data.shape[2] / 8)
print('Dimensions of head model: {:g} x {:g} x {:g} cells'.format(data.shape[0], data.shape[1], data.shape[2] - nzhead))
# Write HDF5 file
headfile = os.path.splitext(args.filename)[0] + '_head.h5'
f = h5py.File(headfile, 'w')
f.attrs['dx, dy, dz'] = dx_dy_dz
f['/data'] = data[:,:,nzhead:data.shape[2]]
f['/data'] = data[:, :, nzhead:data.shape[2]]

34
user_libs/antenna_patterns/initial_save.py
查看文件

@@ -36,37 +36,37 @@ epsr = 5
# Observation radii and angles
radii = np.linspace(0.1, 0.3, 20)
theta = np.linspace(3, 357, 60) * (180/np.pi)
theta = np.linspace(3, 357, 60) * (180 / np.pi)
# Scaling of time-domain field pattern values by material impedance
impscaling = False
# Centre frequency of modelled antenna
f = 1.5e9 # GSSI 1.5GHz antenna model
f = 1.5e9 # GSSI 1.5GHz antenna model
# Largest dimension of antenna transmitting element
D = 0.060 # GSSI 1.5GHz antenna model
D = 0.060 # GSSI 1.5GHz antenna model
# Traces to plot for sanity checking
traceno = np.s_[:] # All traces
traceno = np.s_[:] # All traces
########################################
# Critical angle and velocity
if epsr:
mr = 1
z1 = np.sqrt(mr/epsr) * z0
z1 = np.sqrt(mr / epsr) * z0
v1 = c / np.sqrt(epsr)
thetac = np.round(np.arcsin(v1/c) * (180/np.pi))
wavelength = v1/f
thetac = np.round(np.arcsin(v1 / c) * (180 / np.pi))
wavelength = v1 / f
# Print some useful information
print('Centre frequency: {} GHz'.format(f/1e9))
print('Centre frequency: {} GHz'.format(f / 1e9))
if epsr:
print('Critical angle for Er {} is {} degrees'.format(epsr, thetac))
print('Wavelength: {:.3f} m'.format(wavelength))
print('Observation distance(s) from {:.3f} m ({:.1f} wavelengths) to {:.3f} m ({:.1f} wavelengths)'.format(radii[0], radii[0]/wavelength, radii[-1], radii[-1]/wavelength))
print('Theoretical boundary between reactive & radiating near-field (0.62*sqrt((D^3/wavelength): {:.3f} m'.format(0.62 * np.sqrt((D**3)/wavelength)))
print('Theoretical boundary between radiating near-field & far-field (2*D^2/wavelength): {:.3f} m'.format((2 * D**2)/wavelength))
print('Observation distance(s) from {:.3f} m ({:.1f} wavelengths) to {:.3f} m ({:.1f} wavelengths)'.format(radii[0], radii[0] / wavelength, radii[-1], radii[-1] / wavelength))
print('Theoretical boundary between reactive & radiating near-field (0.62*sqrt((D^3/wavelength): {:.3f} m'.format(0.62 * np.sqrt((D**3) / wavelength)))
print('Theoretical boundary between radiating near-field & far-field (2*D^2/wavelength): {:.3f} m'.format((2 * D**2) / wavelength))
# Load text file with coordinates of pattern origin
origin = np.loadtxt(os.path.splitext(outputfile)[0] + '_rxsorigin.txt')
@@ -77,13 +77,13 @@ iterations = f.attrs['Iterations']
dt = f.attrs['dt']
nrx = f.attrs['nrx']
if antenna:
nrx = nrx - 1 # Ignore first receiver point with full antenna model
nrx = nrx - 1 # Ignore first receiver point with full antenna model
start = 2
else:
start = 1
time = np.arange(0, dt * iterations, dt)
time = time / 1E-9
fs = 1 / dt # Sampling frequency
fs = 1 / dt # Sampling frequency
# Initialise arrays to store fields
coords = np.zeros((nrx, 3), dtype=np.float32)
@@ -115,7 +115,7 @@ for rx in range(0, nrx):
Hz[:, rx] = f[path + 'Hz'][:]
f.close()
## Plot traces for sanity checking
# Plot traces for sanity checking
#fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(num=outputfile, nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), figsize=(20, 10), facecolor='w', edgecolor='w')
#ax1.plot(time, Ex[:, traceno],'r', lw=2)
#ax1.set_ylabel('$E_x$, field strength [V/m]')
@@ -129,12 +129,12 @@ f.close()
#ax4.set_ylabel('$H_y$, field strength [A/m]')
#ax6.plot(time, Hz[:, traceno],'b', lw=2)
#ax6.set_ylabel('$H_z$, field strength [A/m]')
## Turn on grid
# Turn on grid
#[ax.grid() for ax in fig.axes]
#plt.show()
# plt.show()
# Calculate fields for patterns
rxstart = 0 # Index for rx points
rxstart = 0 # Index for rx points
for radius in range(0, len(radii)):
index = 0
# Observation points

33
user_libs/antenna_patterns/plot_fields.py
查看文件

@@ -35,13 +35,13 @@ epsr = 5
# Observation radii and angles
radii = np.linspace(0.1, 0.3, 20)
theta = np.linspace(3, 357, 60)
theta = np.deg2rad(np.append(theta, theta[0])) # Append start value to close circle
theta = np.deg2rad(np.append(theta, theta[0])) # Append start value to close circle
# Centre frequency of modelled antenna
f = 1.5e9 # GSSI 1.5GHz antenna model
f = 1.5e9 # GSSI 1.5GHz antenna model
# Largest dimension of antenna transmitting element
D = 0.060 # GSSI 1.5GHz antenna model
D = 0.060 # GSSI 1.5GHz antenna model
# Minimum value for plotting energy and ring steps (dB)
min = -72
@@ -51,19 +51,19 @@ step = 12
# Critical angle and velocity
if epsr:
mr = 1
z1 = np.sqrt(mr/epsr) * z0
z1 = np.sqrt(mr / epsr) * z0
v1 = c / np.sqrt(epsr)
thetac = np.round(np.rad2deg(np.arcsin(v1/c)))
wavelength = v1/f
thetac = np.round(np.rad2deg(np.arcsin(v1 / c)))
wavelength = v1 / f
# Print some useful information
print('Centre frequency: {} GHz'.format(f/1e9))
print('Centre frequency: {} GHz'.format(f / 1e9))
if epsr:
print('Critical angle for Er {} is {} degrees'.format(epsr, thetac))
print('Wavelength: {:.3f} m'.format(wavelength))
print('Observation distance(s) from {:.3f} m ({:.1f} wavelengths) to {:.3f} m ({:.1f} wavelengths)'.format(radii[0], radii[0]/wavelength, radii[-1], radii[-1]/wavelength))
print('Theoretical boundary between reactive & radiating near-field (0.62*sqrt((D^3/wavelength): {:.3f} m'.format(0.62 * np.sqrt((D**3)/wavelength)))
print('Theoretical boundary between radiating near-field & far-field (2*D^2/wavelength): {:.3f} m'.format((2 * D**2)/wavelength))
print('Observation distance(s) from {:.3f} m ({:.1f} wavelengths) to {:.3f} m ({:.1f} wavelengths)'.format(radii[0], radii[0] / wavelength, radii[-1], radii[-1] / wavelength))
print('Theoretical boundary between reactive & radiating near-field (0.62*sqrt((D^3/wavelength): {:.3f} m'.format(0.62 * np.sqrt((D**3) / wavelength)))
print('Theoretical boundary between radiating near-field & far-field (2*D^2/wavelength): {:.3f} m'.format((2 * D**2) / wavelength))
# Setup figure
fig = plt.figure(num=args.numpyfile, figsize=(8, 8), facecolor='w', edgecolor='w')
@@ -81,15 +81,15 @@ ax.annotate('Ground', xy=(np.deg2rad(270), 0), xytext=(8, -15), textcoords='offs
# Plot patterns
for patt in range(0, len(radii)):
pattplot = np.append(patterns[patt, :], patterns[patt, 0]) # Append start value to close circle
pattplot = pattplot / np.max(np.max(patterns)) # Normalise, based on set of patterns
pattplot = np.append(patterns[patt, :], patterns[patt, 0]) # Append start value to close circle
pattplot = pattplot / np.max(np.max(patterns)) # Normalise, based on set of patterns
ax.plot(theta, 10 * np.log10(pattplot), label='{:.2f}m'.format(radii[patt]), marker='.', ms=6, lw=1.5)
# Add Hertzian dipole plot
#hertzplot1 = np.append(hertzian[0, :], hertzian[0, 0]) # Append start value to close circle
# hertzplot1 = np.append(hertzian[0, :], hertzian[0, 0]) # Append start value to close circle
#hertzplot1 = hertzplot1 / np.max(np.max(hertzian))
#ax.plot(theta, 10 * np.log10(hertzplot1), label='Inf. dipole, 0.1m', color='black', ls='-.', lw=3)
#hertzplot2 = np.append(hertzian[-1, :], hertzian[-1, 0]) # Append start value to close circle
# hertzplot2 = np.append(hertzian[-1, :], hertzian[-1, 0]) # Append start value to close circle
#hertzplot2 = hertzplot2 / np.max(np.max(hertzian))
#ax.plot(theta, 10 * np.log10(hertzplot2), label='Inf. dipole, 0.58m', color='black', ls='--', lw=3)
@@ -103,13 +103,13 @@ ax.set_rmax(0)
ax.set_rlabel_position(45)
ax.set_yticks(np.arange(min, step, step))
yticks = ax.get_yticks().tolist()
yticks[-1]='0 dB'
yticks[-1] = '0 dB'
ax.set_yticklabels(yticks)
# Grid and legend
ax.grid(True)
handles, existlabels = ax.get_legend_handles_labels()
leg = ax.legend([handles[0], handles[-1]], [existlabels[0], existlabels[-1]], ncol=2, loc=(0.27,-0.12), frameon=False) # Plot just first and last legend entries
leg = ax.legend([handles[0], handles[-1]], [existlabels[0], existlabels[-1]], ncol=2, loc=(0.27, -0.12), frameon=False) # Plot just first and last legend entries
#leg = ax.legend([handles[0], handles[-3], handles[-2], handles[-1]], [existlabels[0], existlabels[-3], existlabels[-2], existlabels[-1]], ncol=4, loc=(-0.13,-0.12), frameon=False)
[legobj.set_linewidth(2) for legobj in leg.legendHandles]
@@ -120,4 +120,3 @@ fig.savefig(savename, dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.
#fig.savefig(savename, dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
plt.show()

150
user_libs/antennas.py
查看文件

@@ -12,16 +12,17 @@ from gprMax.input_cmd_funcs import *
moduledirectory = os.path.dirname(os.path.abspath(__file__))
def antenna_like_GSSI_1500(x, y, z, resolution=0.001, rotate90=False, **kwargs):
"""Inserts a description of an antenna similar to the GSSI 1.5GHz antenna. Can be used with 1mm (default) or 2mm spatial resolution. The external dimensions of the antenna are 170x108x45mm. One output point is defined between the arms of the receiever bowtie. The bowties are aligned with the y axis so the output is the y component of the electric field.
Args:
x, y, z (float): Coordinates of a location in the model to insert the antenna. Coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction.
resolution (float): Spatial resolution for the antenna model.
rotate90 (bool): Rotate model 90 degrees CCW in xy plane.
kwargs (dict): Optional variables, e.g. can be fed from an optimisation process.
"""
# Antenna geometry properties
casesize = (0.170, 0.108, 0.043)
casethickness = 0.002
@@ -32,7 +33,7 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, rotate90=False, **kwargs):
bowtiebase = 0.022
bowtieheight = 0.014
patchheight = 0.015
# Set origin for rotation to geometric centre of antenna in x-y plane if required
if rotate90:
rotate90origin = (x, y)
@@ -49,25 +50,25 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, rotate90=False, **kwargs):
absorbersig = kwargs['absorbersig']
rxres = 50
else:
#excitationfreq = 1.5e9 # GHz
#sourceresistance = 50 # Ohms
# excitationfreq = 1.5e9 # GHz
# sourceresistance = 50 # Ohms
#absorberEr = 1.7
#absorbersig = 0.59
# Values from http://hdl.handle.net/1842/4074
excitationfreq = 1.71e9
#sourceresistance = 4
sourceresistance = 230 # Correction for old (< 123) GprMax3D bug
sourceresistance = 230 #  Correction for old (< 123) GprMax3D bug
absorberEr = 1.58
absorbersig = 0.428
rxres = 925 # Resistance at Rx bowtie
rxres = 925 # Resistance at Rx bowtie
x = x - (casesize[0] / 2)
y = y - (casesize[1] / 2)
# Coordinates of source excitation point in antenna
tx = x + 0.114, y + 0.053, z + skidthickness
if resolution == 0.001:
dx = 0.001
dy = 0.001
@@ -92,43 +93,43 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, rotate90=False, **kwargs):
# Plastic case
box(x, y, z + skidthickness, x + casesize[0], y + casesize[1], z + skidthickness + casesize[2], 'hdpe', rotate90origin=rotate90origin)
box(x + casethickness, y + casethickness, z + skidthickness, x + casesize[0] - casethickness, y + casesize[1] - casethickness, z + skidthickness + casesize[2] - casethickness, 'free_space', rotate90origin=rotate90origin)
# Metallic enclosure
box(x + 0.025, y + casethickness, z + skidthickness, x + casesize[0] - 0.025, y + casesize[1] - casethickness, z + skidthickness + 0.027, 'pec', rotate90origin=rotate90origin)
# Absorber material, and foam (modelled as PCB material) around edge of absorber
box(x + 0.025 + shieldthickness, y + casethickness + shieldthickness, z + skidthickness, x + 0.025 + shieldthickness + 0.057, y + casesize[1] - casethickness - shieldthickness, z + skidthickness + 0.027 - shieldthickness - 0.001, 'pcb', rotate90origin=rotate90origin)
box(x + 0.025 + shieldthickness + foamsurroundthickness, y + casethickness + shieldthickness + foamsurroundthickness, z + skidthickness, x + 0.025 + shieldthickness + 0.057 - foamsurroundthickness, y + casesize[1] - casethickness - shieldthickness - foamsurroundthickness, z + skidthickness + 0.027 - shieldthickness, 'absorber', rotate90origin=rotate90origin)
box(x + 0.086, y + casethickness + shieldthickness, z + skidthickness, x + 0.086 + 0.057, y + casesize[1] - casethickness - shieldthickness, z + skidthickness + 0.027 - shieldthickness - 0.001, 'pcb', rotate90origin=rotate90origin)
box(x + 0.086 + foamsurroundthickness, y + casethickness + shieldthickness + foamsurroundthickness, z + skidthickness, x + 0.086 + 0.057 - foamsurroundthickness, y + casesize[1] - casethickness - shieldthickness - foamsurroundthickness, z + skidthickness + 0.027 - shieldthickness, 'absorber', rotate90origin=rotate90origin)
# PCB
box(x + 0.025 + shieldthickness + foamsurroundthickness, y + casethickness + shieldthickness + foamsurroundthickness, z + skidthickness, x + 0.086 - shieldthickness - foamsurroundthickness, y + casesize[1] - casethickness - shieldthickness - foamsurroundthickness, z + skidthickness + pcbthickness, 'pcb', rotate90origin=rotate90origin)
box(x + 0.086 + foamsurroundthickness, y + casethickness + shieldthickness + foamsurroundthickness, z + skidthickness, x + 0.086 + 0.057 - foamsurroundthickness, y + casesize[1] - casethickness - shieldthickness - foamsurroundthickness, z + skidthickness + pcbthickness, 'pcb', rotate90origin=rotate90origin)
# PCB components
if resolution == 0.001:
# Rx & Tx bowties
a = 0
b = 0
while b < 13:
plate(x + 0.045 + a*dx, y + 0.039 + b*dx, z + skidthickness, x + 0.065 - a*dx, y + 0.039 + b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.045 + a*dx, y + 0.067 - b*dx, z + skidthickness, x + 0.065 - a*dx, y + 0.067 - b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a*dx, y + 0.039 + b*dx, z + skidthickness, x + 0.124 - a*dx, y + 0.039 + b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a*dx, y + 0.067 - b*dx, z + skidthickness, x + 0.124 - a*dx, y + 0.067 - b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.045 + a * dx, y + 0.039 + b * dx, z + skidthickness, x + 0.065 - a * dx, y + 0.039 + b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.045 + a * dx, y + 0.067 - b * dx, z + skidthickness, x + 0.065 - a * dx, y + 0.067 - b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a * dx, y + 0.039 + b * dx, z + skidthickness, x + 0.124 - a * dx, y + 0.039 + b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a * dx, y + 0.067 - b * dx, z + skidthickness, x + 0.124 - a * dx, y + 0.067 - b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
b += 1
if a == 2 or a == 4 or a == 7:
plate(x + 0.045 + a*dx, y + 0.039 + b*dx, z + skidthickness, x + 0.065 - a*dx, y + 0.039 + b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.045 + a*dx, y + 0.067 - b*dx, z + skidthickness, x + 0.065 - a*dx, y + 0.067 - b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a*dx, y + 0.039 + b*dx, z + skidthickness, x + 0.124 - a*dx, y + 0.039 + b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a*dx, y + 0.067 - b*dx, z + skidthickness, x + 0.124 - a*dx, y + 0.067 - b*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.045 + a * dx, y + 0.039 + b * dx, z + skidthickness, x + 0.065 - a * dx, y + 0.039 + b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.045 + a * dx, y + 0.067 - b * dx, z + skidthickness, x + 0.065 - a * dx, y + 0.067 - b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a * dx, y + 0.039 + b * dx, z + skidthickness, x + 0.124 - a * dx, y + 0.039 + b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.104 + a * dx, y + 0.067 - b * dx, z + skidthickness, x + 0.124 - a * dx, y + 0.067 - b * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
b += 1
a += 1
# Rx extension section (upper y)
plate(x + 0.044, y + 0.068, z + skidthickness, x + 0.044 + bowtiebase, y + 0.068 + patchheight, z + skidthickness, 'pec', rotate90origin=rotate90origin)
# Tx extension section (upper y)
plate(x + 0.103, y + 0.068, z + skidthickness, x + 0.103 + bowtiebase, y + 0.068 + patchheight, z + skidthickness, 'pec', rotate90origin=rotate90origin)
# Edges that represent wire between bowtie halves in 1mm model
edge(tx[0] - 0.059, tx[1] - dy, tx[2], tx[0] - 0.059, tx[1], tx[2], 'pec', rotate90origin=rotate90origin)
edge(tx[0] - 0.059, tx[1] + dy, tx[2], tx[0] - 0.059, tx[1] + 0.002, tx[2], 'pec', rotate90origin=rotate90origin)
@@ -136,12 +137,12 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, rotate90=False, **kwargs):
edge(tx[0], tx[1] + dz, tx[2], tx[0], tx[1] + 0.002, tx[2], 'pec', rotate90origin=rotate90origin)
elif resolution == 0.002:
# Rx & Tx bowties
for a in range(0,6):
plate(x + 0.044 + a*dx, y + 0.040 + a*dx, z + skidthickness, x + 0.066 - a*dx, y + 0.040 + a*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.044 + a*dx, y + 0.064 - a*dx, z + skidthickness, x + 0.066 - a*dx, y + 0.064 - a*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.103 + a*dx, y + 0.040 + a*dx, z + skidthickness, x + 0.125 - a*dx, y + 0.040 + a*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.103 + a*dx, y + 0.064 - a*dx, z + skidthickness, x + 0.125 - a*dx, y + 0.064 - a*dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
# Rx & Tx bowties
for a in range(0, 6):
plate(x + 0.044 + a * dx, y + 0.040 + a * dx, z + skidthickness, x + 0.066 - a * dx, y + 0.040 + a * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.044 + a * dx, y + 0.064 - a * dx, z + skidthickness, x + 0.066 - a * dx, y + 0.064 - a * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.103 + a * dx, y + 0.040 + a * dx, z + skidthickness, x + 0.125 - a * dx, y + 0.040 + a * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
plate(x + 0.103 + a * dx, y + 0.064 - a * dx, z + skidthickness, x + 0.125 - a * dx, y + 0.064 - a * dx + dy, z + skidthickness, 'pec', rotate90origin=rotate90origin)
# Rx extension section (upper y)
plate(x + 0.044, y + 0.066, z + skidthickness, x + 0.044 + bowtiebase, y + 0.066 + patchheight, z + skidthickness, 'pec', rotate90origin=rotate90origin)
# Tx extension section (upper y)
@@ -154,19 +155,19 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, rotate90=False, **kwargs):
# Skid
box(x, y, z, x + casesize[0], y + casesize[1], z + skidthickness, 'hdpe', rotate90origin=rotate90origin)
# Geometry views
#geometry_view(x - dx, y - dy, z - dz, x + casesize[0] + dx, y + casesize[1] + dy, z + skidthickness + casesize[2] + dz, dx, dy, dz, 'antenna_like_GSSI_1500')
#geometry_view(x, y, z, x + casesize[0], y + casesize[1], z + 0.010, dx, dy, dz, 'antenna_like_GSSI_1500_pcb', type='f')
# Excitation - custom pulse
#print('#excitation_file: {}'.format(os.path.join(moduledirectory, 'GSSIgausspulse1.txt')))
#print('#transmission_line: y {} {} {} {} GSSIgausspulse1'.format(tx[0], tx[1], tx[2], sourceresistance))
# print('#excitation_file: {}'.format(os.path.join(moduledirectory, 'GSSIgausspulse1.txt')))
# print('#transmission_line: y {} {} {} {} GSSIgausspulse1'.format(tx[0], tx[1], tx[2], sourceresistance))
# Excitation - Gaussian pulse
print('#waveform: gaussian 1 {} myGaussian'.format(excitationfreq))
transmission_line('y', tx[0], tx[1], tx[2], sourceresistance, 'myGaussian', dxdy=(resolution, resolution), rotate90origin=rotate90origin)
# Output point - receiver bowtie
if resolution == 0.001:
edge(tx[0] - 0.059, tx[1], tx[2], tx[0] - 0.059, tx[1] + dy, tx[2], 'rxres', rotate90origin=rotate90origin)
@@ -178,25 +179,25 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, rotate90=False, **kwargs):
def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
"""Inserts a description of an antenna similar to the MALA 1.2GHz antenna. Can be used with 1mm (default) or 2mm spatial resolution. The external dimensions of the antenna are 184x109x46mm. One output point is defined between the arms of the receiever bowtie. The bowties are aligned with the y axis so the output is the y component of the electric field.
Args:
x, y, z (float): Coordinates of a location in the model to insert the antenna. Coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction.
resolution (float): Spatial resolution for the antenna model.
rotate90 (bool): Rotate model 90 degrees CCW in xy plane.
kwargs (dict): Optional variables, e.g. can be fed from an optimisation process.
"""
# Antenna geometry properties
casesize = (0.184, 0.109, 0.040)
casethickness = 0.002
cavitysize = (0.062, 0.062, 0.037)
cavitythickness = 0.001
pcbthickness = 0.002
polypropylenethickness = 0.003;
hdpethickness = 0.003;
polypropylenethickness = 0.003
hdpethickness = 0.003
skidthickness = 0.006
bowtieheight = 0.025
# Set origin for rotation to geometric centre of antenna in x-y plane if required
if rotate90:
rotate90origin = (x, y)
@@ -204,7 +205,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
else:
rotate90origin = ()
output = 'Ey'
# Unknown properties
if kwargs:
excitationfreq = kwargs['excitationfreq']
@@ -217,13 +218,13 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
sourceresistance = 1000
absorberEr = 6.49
absorbersig = 0.252
x = x - (casesize[0] / 2)
y = y - (casesize[1] / 2)
# Coordinates of source excitation point in antenna
tx = x + 0.063, y + 0.052, z + skidthickness
if resolution == 0.001:
dx = 0.001
dy = 0.001
@@ -234,25 +235,25 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
dz = 0.002
cavitysize = (0.062, 0.062, 0.036)
cavitythickness = 0.002
polypropylenethickness = 0.002;
hdpethickness = 0.004;
polypropylenethickness = 0.002
hdpethickness = 0.004
bowtieheight = 0.024
tx = x + 0.062, y + 0.052, z + skidthickness
else:
raise CmdInputError('This antenna module can only be used with a spatial resolution of 1mm or 2mm')
# SMD resistors - 3 on each Tx & Rx bowtie arm
txres = 470 # Ohms
txrescellupper = txres / 3 # Resistor over 3 cells
txsigupper = ((1 / txrescellupper) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
txrescelllower = txres / 4 # Resistor over 4 cells
txsiglower = ((1 / txrescelllower) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
rxres = 150 # Ohms
rxrescellupper = rxres / 3 # Resistor over 3 cells
rxsigupper = ((1 / rxrescellupper) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
rxrescelllower = rxres / 4 # Resistor over 4 cells
rxsiglower = ((1 / rxrescelllower) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
txres = 470 # Ohms
txrescellupper = txres / 3 # Resistor over 3 cells
txsigupper = ((1 / txrescellupper) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
txrescelllower = txres / 4 # Resistor over 4 cells
txsiglower = ((1 / txrescelllower) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
rxres = 150 # Ohms
rxrescellupper = rxres / 3 # Resistor over 3 cells
rxsigupper = ((1 / rxrescellupper) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
rxrescelllower = rxres / 4 # Resistor over 4 cells
rxsiglower = ((1 / rxrescelllower) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor
# Material definitions
material(absorberEr, absorbersig, 1, 0, 'absorber')
material(3, 0, 1, 0, 'pcb')
@@ -262,17 +263,17 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
material(3, txsigupper, 1, 0, 'txresupper')
material(3, rxsiglower, 1, 0, 'rxreslower')
material(3, rxsigupper, 1, 0, 'rxresupper')
# Antenna geometry
# Shield - metallic enclosure
box(x, y, z + skidthickness, x + casesize[0], y + casesize[1], z + skidthickness + casesize[2], 'pec', rotate90origin=rotate90origin)
box(x + 0.020, y + casethickness, z + skidthickness, x + 0.100, y + casesize[1] - casethickness, z + skidthickness + casethickness, 'free_space', rotate90origin=rotate90origin)
box(x + 0.100, y + casethickness, z + skidthickness, x + casesize[0] - casethickness, y + casesize[1] - casethickness, z + skidthickness + casethickness, 'free_space', rotate90origin=rotate90origin)
# Absorber material
box(x + 0.020, y + casethickness, z + skidthickness, x + 0.100, y + casesize[1] - casethickness, z + skidthickness + casesize[2] - casethickness, 'absorber', rotate90origin=rotate90origin)
box(x + 0.100, y + casethickness, z + skidthickness, x + casesize[0] - casethickness, y + casesize[1] - casethickness, z + skidthickness + casesize[2] - casethickness, 'absorber', rotate90origin=rotate90origin)
# Shield - cylindrical sections
cylinder(x + 0.055, y + casesize[1] - 0.008, z + skidthickness, x + 0.055, y + casesize[1] - 0.008, z + skidthickness + casesize[2] - casethickness, 0.008, 'pec', rotate90origin=rotate90origin)
cylinder(x + 0.055, y + 0.008, z + skidthickness, x + 0.055, y + 0.008, z + skidthickness + casesize[2] - casethickness, 0.008, 'pec', rotate90origin=rotate90origin)
@@ -286,24 +287,24 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
box(x + 0.054, y + 0.014, z + skidthickness, x + 0.056, y + 0.016, z + skidthickness + casesize[2] - casethickness, 'free_space', rotate90origin=rotate90origin)
box(x + 0.146, y + casesize[1] - 0.016, z + skidthickness, x + 0.148, y + casesize[1] - 0.014, z + skidthickness + casesize[2] - casethickness, 'free_space', rotate90origin=rotate90origin)
box(x + 0.146, y + 0.014, z + skidthickness, x + 0.148, y + 0.016, z + skidthickness + casesize[2] - casethickness, 'free_space', rotate90origin=rotate90origin)
# PCB
box(x + 0.020, y + 0.018, z + skidthickness, x + casesize[0] - casethickness, y + casesize[1] - 0.018, z + skidthickness + pcbthickness, 'pcb', rotate90origin=rotate90origin)
# Shield - Tx & Rx cavities
box(x + 0.032, y + 0.022, z + skidthickness, x + 0.032 + cavitysize[0], y + 0.022 + cavitysize[1], z + skidthickness + cavitysize[2], 'pec', rotate90origin=rotate90origin)
box(x + 0.032 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness, x + 0.032 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, z + skidthickness + cavitysize[2], 'absorber', rotate90origin=rotate90origin)
box(x + 0.108, y + 0.022, z + skidthickness, x + 0.108 + cavitysize[0], y + 0.022 + cavitysize[1], z + skidthickness + cavitysize[2], 'pec', rotate90origin=rotate90origin)
box(x + 0.108 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness, x + 0.108 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, z + skidthickness + cavitysize[2], 'free_space', rotate90origin=rotate90origin)
# Shield - Tx & Rx cavities - joining strips
box(x + 0.032 + cavitysize[0], y + 0.022 + cavitysize[1] - 0.006, z + skidthickness + cavitysize[2] - casethickness, x + 0.108, y + 0.022 + cavitysize[1], z + skidthickness + cavitysize[2], 'pec', rotate90origin=rotate90origin)
box(x + 0.032 + cavitysize[0], y + 0.022, z + skidthickness + cavitysize[2] - casethickness, x + 0.108, y + 0.022 + 0.006, z + skidthickness + cavitysize[2], 'pec', rotate90origin=rotate90origin)
# PCB - replace bits chopped by TX & Rx cavities
box(x + 0.032 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness, x + 0.032 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, z + skidthickness + pcbthickness, 'pcb', rotate90origin=rotate90origin)
box(x + 0.108 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness, x + 0.108 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, z + skidthickness + pcbthickness, 'pcb', rotate90origin=rotate90origin)
# PCB components
# Tx bowtie
if resolution == 0.001:
@@ -314,7 +315,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
elif resolution == 0.002:
triangle(tx[0], tx[1], tx[2], tx[0] - 0.026, tx[1] - bowtieheight, tx[2], tx[0] + 0.026, tx[1] - bowtieheight, tx[2], 0, 'pec', rotate90origin=rotate90origin)
triangle(tx[0], tx[1] + 0.002, tx[2], tx[0] - 0.026, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + 0.026, tx[1] + bowtieheight + 0.002, tx[2], 0, 'pec', rotate90origin=rotate90origin)
# Rx bowtie
if resolution == 0.001:
triangle(tx[0] + 0.076, tx[1] - 0.001, tx[2], tx[0] + 0.076 - 0.026, tx[1] - bowtieheight - 0.001, tx[2], tx[0] + 0.076 + 0.026, tx[1] - bowtieheight - 0.001, tx[2], 0, 'pec', rotate90origin=rotate90origin)
@@ -324,7 +325,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
elif resolution == 0.002:
triangle(tx[0] + 0.076, tx[1], tx[2], tx[0] + 0.076 - 0.026, tx[1] - bowtieheight, tx[2], tx[0] + 0.076 + 0.026, tx[1] - bowtieheight, tx[2], 0, 'pec', rotate90origin=rotate90origin)
triangle(tx[0] + 0.076, tx[1] + 0.002, tx[2], tx[0] + 0.076 - 0.026, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + 0.076 + 0.026, tx[1] + bowtieheight + 0.002, tx[2], 0, 'pec', rotate90origin=rotate90origin)
# Tx surface mount resistors (lower y coordinate)
if resolution == 0.001:
edge(tx[0] - 0.023, tx[1] - bowtieheight - 0.004, tx[2], tx[0] - 0.023, tx[1] - bowtieheight - dy, tx[2], 'txreslower', rotate90origin=rotate90origin)
@@ -340,7 +341,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
edge(tx[0] + dx, tx[1] - bowtieheight - 0.004, tx[2], tx[0] + dx, tx[1] - bowtieheight, tx[2], 'txreslower', rotate90origin=rotate90origin)
edge(tx[0] + 0.020, tx[1] - bowtieheight - 0.004, tx[2], tx[0] + 0.020, tx[1] - bowtieheight, tx[2], 'txreslower', rotate90origin=rotate90origin)
edge(tx[0] + 0.020 + dx, tx[1] - bowtieheight - 0.004, tx[2], tx[0] + 0.020 + dx, tx[1] - bowtieheight, tx[2], 'txreslower', rotate90origin=rotate90origin)
# Tx surface mount resistors (upper y coordinate)
if resolution == 0.001:
edge(tx[0] - 0.023, tx[1] + bowtieheight + 0.002, tx[2], tx[0] - 0.023, tx[1] + bowtieheight + 0.006, tx[2], 'txresupper', rotate90origin=rotate90origin)
@@ -356,7 +357,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
edge(tx[0] + dx, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + dx, tx[1] + bowtieheight + 0.006, tx[2], 'txresupper', rotate90origin=rotate90origin)
edge(tx[0] + 0.020, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + 0.020, tx[1] + bowtieheight + 0.006, tx[2], 'txresupper', rotate90origin=rotate90origin)
edge(tx[0] + 0.020 + dx, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + 0.020 + dx, tx[1] + bowtieheight + 0.006, tx[2], 'txresupper', rotate90origin=rotate90origin)
# Rx surface mount resistors (lower y coordinate)
if resolution == 0.001:
edge(tx[0] - 0.023 + 0.076, tx[1] - bowtieheight - 0.004, tx[2], tx[0] - 0.023 + 0.076, tx[1] - bowtieheight - dy, tx[2], 'rxreslower', rotate90origin=rotate90origin)
@@ -372,7 +373,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
edge(tx[0] + dx + 0.076, tx[1] - bowtieheight - 0.004, tx[2], tx[0] + dx + 0.076, tx[1] - bowtieheight, tx[2], 'rxreslower', rotate90origin=rotate90origin)
edge(tx[0] + 0.020 + 0.076, tx[1] - bowtieheight - 0.004, tx[2], tx[0] + 0.020 + 0.076, tx[1] - bowtieheight, tx[2], 'rxreslower', rotate90origin=rotate90origin)
edge(tx[0] + 0.020 + dx + 0.076, tx[1] - bowtieheight - 0.004, tx[2], tx[0] + 0.020 + dx + 0.076, tx[1] - bowtieheight, tx[2], 'rxreslower', rotate90origin=rotate90origin)
# Rx surface mount resistors (upper y coordinate)
if resolution == 0.001:
edge(tx[0] - 0.023 + 0.076, tx[1] + bowtieheight + 0.002, tx[2], tx[0] - 0.023 + 0.076, tx[1] + bowtieheight + 0.006, tx[2], 'rxresupper', rotate90origin=rotate90origin)
@@ -388,7 +389,7 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
edge(tx[0] + dx + 0.076, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + dx + 0.076, tx[1] + bowtieheight + 0.006, tx[2], 'rxresupper', rotate90origin=rotate90origin)
edge(tx[0] + 0.020 + 0.076, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + 0.020 + 0.076, tx[1] + bowtieheight + 0.006, tx[2], 'rxresupper', rotate90origin=rotate90origin)
edge(tx[0] + 0.020 + dx + 0.076, tx[1] + bowtieheight + 0.002, tx[2], tx[0] + 0.020 + dx + 0.076, tx[1] + bowtieheight + 0.006, tx[2], 'rxresupper', rotate90origin=rotate90origin)
# Skid
box(x, y, z, x + casesize[0], y + casesize[1], z + polypropylenethickness, 'polypropylene', rotate90origin=rotate90origin)
box(x, y, z + polypropylenethickness, x + casesize[0], y + casesize[1], z + polypropylenethickness + hdpethickness, 'hdpe', rotate90origin=rotate90origin)
@@ -396,11 +397,10 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, rotate90=False, **kwargs):
# Geometry views
#geometry_view(x - dx, y - dy, z - dz, x + casesize[0] + dx, y + casesize[1] + dy, z + casesize[2] + skidthickness + dz, dx, dy, dz, 'antenna_like_MALA_1200')
#geometry_view(x, y, z, x + casesize[0], y + casesize[1], z + 0.010, dx, dy, dz, 'antenna_like_MALA_1200_pcb', type='f')
# Excitation
print('#waveform: gaussian 1.0 {} myGaussian'.format(excitationfreq))
voltage_source('y', tx[0], tx[1], tx[2], sourceresistance, 'myGaussian', dxdy=(resolution, resolution), rotate90origin=rotate90origin)
# Output point - receiver bowtie
rx(tx[0] + 0.076, tx[1], tx[2], identifier='rxbowtie', to_save=[output], polarisation='y', dxdy=(resolution, resolution), rotate90origin=rotate90origin)

21
user_libs/optimisation_taguchi/fitness_functions.py
查看文件

@@ -24,6 +24,7 @@ from gprMax.exceptions import GeneralError
The second argument is a dictionary which can contain any number of additional arguments, e.g. names (IDs) of outputs (rxs) from input file
"""
def min_max_value(filename, args):
"""Minimum value from a response.
@@ -78,8 +79,8 @@ def xcorr(filename, args):
raise GeneralError('Cannot load reference response from {}'.format(refrespfile))
with open(refrespfile, 'r') as f:
refdata = np.loadtxt(f)
reftime = refdata[:,0] * 1e-9
refresp = refdata[:,1]
reftime = refdata[:, 0] * 1e-9
refresp = refdata[:, 1]
# Load response from output file
f = h5py.File(filename, 'r')
@@ -131,8 +132,8 @@ def xcorr(filename, args):
#fig, ax = plt.subplots(subplot_kw=dict(xlabel='Iterations', ylabel='Voltage [V]'), figsize=(20, 10), facecolor='w', edgecolor='w')
#ax.plot(refresp,'r', lw=2, label='refresp')
#ax.plot(modelresp,'b', lw=2, label='modelresp')
#ax.grid()
#plt.show()
# ax.grid()
# plt.show()
# Calculate cross-correlation
xcorr = signal.correlate(refresp, modelresp)
@@ -224,7 +225,7 @@ def compactness(filename, args):
peak = np.amax([np.amax(outputdata), np.abs(np.amin(outputdata))])
# Get peaks and troughs in signal
delta = peak / 50 # Considered a peak/trough if it has the max/min value, and was preceded (to the left) by a value lower by delta
delta = peak / 50 # Considered a peak/trough if it has the max/min value, and was preceded (to the left) by a value lower by delta
maxtab, mintab = peakdet(outputdata, delta)
peaks = maxtab + mintab
peaks.sort()
@@ -271,8 +272,8 @@ def spectral_centroid(x, samplerate):
magnitudes = np.abs(np.fft.rfft(x))
length = len(x)
# Positive frequencies
freqs = np.abs(np.fft.fftfreq(length, 1.0/samplerate)[:length//2+1])
centroid = np.sum(magnitudes*freqs) / np.sum(magnitudes)
freqs = np.abs(np.fft.fftfreq(length, 1.0 / samplerate)[:length // 2 + 1])
centroid = np.sum(magnitudes * freqs) / np.sum(magnitudes)
return centroid
@@ -294,7 +295,7 @@ def zero_crossings(x):
return indexzeros
def peakdet(v, delta, x = None):
def peakdet(v, delta, x=None):
"""Detect peaks and troughs in a vector (adapted from MATLAB script at http://billauer.co.il/peakdet.html).
A point is considered a maximum peak if it has the maximal value, and was preceded (to the left) by a value lower by delta.
@@ -341,14 +342,14 @@ def peakdet(v, delta, x = None):
mnpos = x[i]
if lookformax:
if this < mx-delta:
if this < mx - delta:
if int(mxpos) != 0:
maxtab.append(int(mxpos))
mn = this
mnpos = x[i]
lookformax = False
else:
if this > mn+delta:
if this > mn + delta:
if int(mnpos) != 0:
mintab.append(int(mnpos))
mx = this

4
user_libs/optimisation_taguchi/plot_results.py
查看文件

@@ -32,7 +32,7 @@ for item in optparamsinit:
print('History of parameter values:')
for key, value in optparamshist.items():
print('\t{}: {}'.format(key, value))
# Plot the history of fitness values and each optimised parameter values for the optimisation
plot_optimisation_history(fitnessvalueshist, optparamshist, optparamsinit)