你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-03 19:26:50 +08:00
Removing use_line_collection in plotting
use_line_collection is removed in matplotlib 3.8.0
这个提交包含在：
JingziC
@@ -1,237 +1,237 @@
|
||||
# Copyright (C) 2015-2023: The University of Edinburgh
|
||||
# Authors: Craig Warren and Antonis Giannopoulos
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# 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/>.
|
||||
|
||||
import argparse
|
||||
import os
|
||||
import sys
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
|
||||
from gprMax.exceptions import CmdInputError
|
||||
from gprMax.receivers import Rx
|
||||
from gprMax.utilities import fft_power
|
||||
|
||||
|
||||
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']
|
||||
dt = f.attrs['dt']
|
||||
iterations = f.attrs['Iterations']
|
||||
time = np.linspace(0, (iterations - 1) * dt, num=iterations)
|
||||
|
||||
# Check there are any receivers
|
||||
if nrx == 0:
|
||||
raise CmdInputError('No receivers found in {}'.format(filename))
|
||||
|
||||
# Check for single output component when doing a FFT
|
||||
if fft:
|
||||
if not len(outputs) == 1:
|
||||
raise CmdInputError('A single output must be specified when using the -fft option')
|
||||
|
||||
# New plot for each receiver
|
||||
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
|
||||
outputtext = '-' + outputs[0][0:-1]
|
||||
output = outputs[0][0:-1]
|
||||
else:
|
||||
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
|
||||
if fft:
|
||||
# FFT
|
||||
freqs, power = fft_power(outputdata, dt)
|
||||
freqmaxpower = np.where(np.isclose(power, 0))[0][0]
|
||||
|
||||
# Set plotting range to -60dB from maximum power or 4 times
|
||||
# frequency at maximum power
|
||||
try:
|
||||
pltrange = np.where(power[freqmaxpower:] < -60)[0][0] + freqmaxpower + 1
|
||||
except:
|
||||
pltrange = freqmaxpower * 4
|
||||
|
||||
pltrange = np.s_[0:pltrange]
|
||||
|
||||
# Plot time history of output component
|
||||
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
line1 = ax1.plot(time, outputdata, 'r', lw=2, label=outputtext)
|
||||
ax1.set_xlabel('Time [s]')
|
||||
ax1.set_ylabel(outputtext + ' field strength [V/m]')
|
||||
ax1.set_xlim([0, np.amax(time)])
|
||||
ax1.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra
|
||||
markerline, stemlines, baseline = ax2.stem(freqs[pltrange], power[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
line2 = ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2)
|
||||
ax2.set_xlabel('Frequency [Hz]')
|
||||
ax2.set_ylabel('Power [dB]')
|
||||
ax2.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Change colours and labels for magnetic field components or currents
|
||||
if 'H' in outputs[0]:
|
||||
plt.setp(line1, color='g')
|
||||
plt.setp(line2, color='g')
|
||||
plt.setp(ax1, ylabel=outputtext + ' field strength [A/m]')
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
elif 'I' in outputs[0]:
|
||||
plt.setp(line1, color='b')
|
||||
plt.setp(line2, color='b')
|
||||
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)
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
# ax.set_ylim([-15, 20])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
if 'H' in output:
|
||||
plt.setp(line, color='g')
|
||||
plt.setp(ax, ylabel=outputtext + ', field strength [A/m]')
|
||||
elif 'I' in output:
|
||||
plt.setp(line, color='b')
|
||||
plt.setp(ax, ylabel=outputtext + ', current [A]')
|
||||
|
||||
# If multiple outputs required, create all nine subplots and populate only the specified ones
|
||||
else:
|
||||
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
if len(outputs) == 9:
|
||||
gs = gridspec.GridSpec(3, 3, hspace=0.3, wspace=0.3)
|
||||
else:
|
||||
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
|
||||
outputtext = '-' + output[0:-1]
|
||||
output = output[0:-1]
|
||||
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.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.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.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.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.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.set_ylabel(outputtext + ', field strength [A/m]')
|
||||
# ax.set_ylim([-0.03, 0.03])
|
||||
elif output == 'Ix':
|
||||
ax = plt.subplot(gs[0, 2])
|
||||
ax.plot(time, outputdata, 'b', lw=2, label=outputtext)
|
||||
ax.set_ylabel(outputtext + ', current [A]')
|
||||
elif output == 'Iy':
|
||||
ax = plt.subplot(gs[1, 2])
|
||||
ax.plot(time, outputdata, 'b', lw=2, label=outputtext)
|
||||
ax.set_ylabel(outputtext + ', current [A]')
|
||||
elif output == 'Iz':
|
||||
ax = plt.subplot(gs[2, 2])
|
||||
ax.plot(time, outputdata, 'b', lw=2, label=outputtext)
|
||||
ax.set_ylabel(outputtext + ', current [A]')
|
||||
for ax in fig.axes:
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Save a PDF/PNG of the figure
|
||||
# fig.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_rx' + str(rx) + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
# fig.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_rx' + str(rx) + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
|
||||
f.close()
|
||||
|
||||
return plt
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(description='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.', usage='cd gprMax; python -m tools.plot_Ascan outputfile')
|
||||
parser.add_argument('outputfile', help='name of output file including path')
|
||||
parser.add_argument('--outputs', help='outputs to be plotted', default=Rx.defaultoutputs, choices=['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz', 'Ex-', 'Ey-', 'Ez-', 'Hx-', 'Hy-', 'Hz-', 'Ix-', 'Iy-', 'Iz-'], nargs='+')
|
||||
parser.add_argument('-fft', action='store_true', help='plot FFT (single output must be specified)', default=False)
|
||||
args = parser.parse_args()
|
||||
|
||||
plthandle = mpl_plot(args.outputfile, args.outputs, fft=args.fft)
|
||||
plthandle.show()
|
||||
# Copyright (C) 2015-2023: The University of Edinburgh
|
||||
# Authors: Craig Warren and Antonis Giannopoulos
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# 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/>.
|
||||
|
||||
import argparse
|
||||
import os
|
||||
import sys
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
|
||||
from gprMax.exceptions import CmdInputError
|
||||
from gprMax.receivers import Rx
|
||||
from gprMax.utilities import fft_power
|
||||
|
||||
|
||||
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']
|
||||
dt = f.attrs['dt']
|
||||
iterations = f.attrs['Iterations']
|
||||
time = np.linspace(0, (iterations - 1) * dt, num=iterations)
|
||||
|
||||
# Check there are any receivers
|
||||
if nrx == 0:
|
||||
raise CmdInputError('No receivers found in {}'.format(filename))
|
||||
|
||||
# Check for single output component when doing a FFT
|
||||
if fft:
|
||||
if not len(outputs) == 1:
|
||||
raise CmdInputError('A single output must be specified when using the -fft option')
|
||||
|
||||
# New plot for each receiver
|
||||
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
|
||||
outputtext = '-' + outputs[0][0:-1]
|
||||
output = outputs[0][0:-1]
|
||||
else:
|
||||
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
|
||||
if fft:
|
||||
# FFT
|
||||
freqs, power = fft_power(outputdata, dt)
|
||||
freqmaxpower = np.where(np.isclose(power, 0))[0][0]
|
||||
|
||||
# Set plotting range to -60dB from maximum power or 4 times
|
||||
# frequency at maximum power
|
||||
try:
|
||||
pltrange = np.where(power[freqmaxpower:] < -60)[0][0] + freqmaxpower + 1
|
||||
except:
|
||||
pltrange = freqmaxpower * 4
|
||||
|
||||
pltrange = np.s_[0:pltrange]
|
||||
|
||||
# Plot time history of output component
|
||||
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
line1 = ax1.plot(time, outputdata, 'r', lw=2, label=outputtext)
|
||||
ax1.set_xlabel('Time [s]')
|
||||
ax1.set_ylabel(outputtext + ' field strength [V/m]')
|
||||
ax1.set_xlim([0, np.amax(time)])
|
||||
ax1.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra
|
||||
markerline, stemlines, baseline = ax2.stem(freqs[pltrange], power[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
line2 = ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2)
|
||||
ax2.set_xlabel('Frequency [Hz]')
|
||||
ax2.set_ylabel('Power [dB]')
|
||||
ax2.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Change colours and labels for magnetic field components or currents
|
||||
if 'H' in outputs[0]:
|
||||
plt.setp(line1, color='g')
|
||||
plt.setp(line2, color='g')
|
||||
plt.setp(ax1, ylabel=outputtext + ' field strength [A/m]')
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
elif 'I' in outputs[0]:
|
||||
plt.setp(line1, color='b')
|
||||
plt.setp(line2, color='b')
|
||||
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)
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
# ax.set_ylim([-15, 20])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
if 'H' in output:
|
||||
plt.setp(line, color='g')
|
||||
plt.setp(ax, ylabel=outputtext + ', field strength [A/m]')
|
||||
elif 'I' in output:
|
||||
plt.setp(line, color='b')
|
||||
plt.setp(ax, ylabel=outputtext + ', current [A]')
|
||||
|
||||
# If multiple outputs required, create all nine subplots and populate only the specified ones
|
||||
else:
|
||||
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
if len(outputs) == 9:
|
||||
gs = gridspec.GridSpec(3, 3, hspace=0.3, wspace=0.3)
|
||||
else:
|
||||
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
|
||||
outputtext = '-' + output[0:-1]
|
||||
output = output[0:-1]
|
||||
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.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.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.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.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.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.set_ylabel(outputtext + ', field strength [A/m]')
|
||||
# ax.set_ylim([-0.03, 0.03])
|
||||
elif output == 'Ix':
|
||||
ax = plt.subplot(gs[0, 2])
|
||||
ax.plot(time, outputdata, 'b', lw=2, label=outputtext)
|
||||
ax.set_ylabel(outputtext + ', current [A]')
|
||||
elif output == 'Iy':
|
||||
ax = plt.subplot(gs[1, 2])
|
||||
ax.plot(time, outputdata, 'b', lw=2, label=outputtext)
|
||||
ax.set_ylabel(outputtext + ', current [A]')
|
||||
elif output == 'Iz':
|
||||
ax = plt.subplot(gs[2, 2])
|
||||
ax.plot(time, outputdata, 'b', lw=2, label=outputtext)
|
||||
ax.set_ylabel(outputtext + ', current [A]')
|
||||
for ax in fig.axes:
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Save a PDF/PNG of the figure
|
||||
# fig.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_rx' + str(rx) + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
# fig.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_rx' + str(rx) + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
|
||||
f.close()
|
||||
|
||||
return plt
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(description='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.', usage='cd gprMax; python -m tools.plot_Ascan outputfile')
|
||||
parser.add_argument('outputfile', help='name of output file including path')
|
||||
parser.add_argument('--outputs', help='outputs to be plotted', default=Rx.defaultoutputs, choices=['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz', 'Ex-', 'Ey-', 'Ez-', 'Hx-', 'Hy-', 'Hz-', 'Ix-', 'Iy-', 'Iz-'], nargs='+')
|
||||
parser.add_argument('-fft', action='store_true', help='plot FFT (single output must be specified)', default=False)
|
||||
args = parser.parse_args()
|
||||
|
||||
plthandle = mpl_plot(args.outputfile, args.outputs, fft=args.fft)
|
||||
plthandle.show()
|
||||
|
||||
@@ -1,424 +1,424 @@
|
||||
# Copyright (C) 2015-2023: The University of Edinburgh
|
||||
# Authors: Craig Warren and Antonis Giannopoulos
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# 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/>.
|
||||
|
||||
import argparse
|
||||
import os
|
||||
import sys
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
|
||||
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.
|
||||
"""
|
||||
|
||||
# Open output file and read some attributes
|
||||
f = h5py.File(filename, 'r')
|
||||
dxdydz = f.attrs['dx_dy_dz']
|
||||
dt = f.attrs['dt']
|
||||
iterations = f.attrs['Iterations']
|
||||
|
||||
# Calculate time array and frequency bin spacing
|
||||
time = np.linspace(0, (iterations - 1) * dt, num=iterations)
|
||||
df = 1 / np.amax(time)
|
||||
|
||||
print('Time window: {:g} s ({} iterations)'.format(np.amax(time), iterations))
|
||||
print('Time step: {:g} s'.format(dt))
|
||||
print('Frequency bin spacing: {:g} Hz'.format(df))
|
||||
|
||||
# Read/calculate voltages and currents from transmitter antenna
|
||||
tltxpath = '/tls/tl' + str(tltxnumber) + '/'
|
||||
|
||||
# Incident voltages/currents
|
||||
Vinc = f[tltxpath + 'Vinc'][:]
|
||||
Iinc = f[tltxpath + 'Iinc'][:]
|
||||
|
||||
# Total (incident + reflected) voltages/currents
|
||||
Vtotal = f[tltxpath + 'Vtotal'][:]
|
||||
Itotal = f[tltxpath + 'Itotal'][:]
|
||||
|
||||
# Reflected voltages/currents
|
||||
Vref = Vtotal - Vinc
|
||||
Iref = Itotal - Iinc
|
||||
|
||||
# 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'][:]
|
||||
|
||||
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]
|
||||
elif rxcomponent == 'Ey':
|
||||
Vrec = f[rxpath][:] * -1 * dxdydz[1]
|
||||
elif rxcomponent == 'Ez':
|
||||
Vrec = f[rxpath][:] * -1 * dxdydz[2]
|
||||
f.close()
|
||||
|
||||
# Frequency bins
|
||||
freqs = np.fft.fftfreq(Vinc.size, d=dt)
|
||||
|
||||
# Delay correction - current lags voltage, so delay voltage to match current timestep
|
||||
delaycorrection = np.exp(1j * 2 * np.pi * freqs * (dt / 2))
|
||||
|
||||
# Calculate s11 and (optionally) s21
|
||||
with np.errstate(divide='ignore'):
|
||||
s11 = np.abs(np.fft.fft(Vref) / np.fft.fft(Vinc))
|
||||
if tlrxnumber or rxnumber:
|
||||
with np.errstate(divide='ignore'):
|
||||
s21 = np.abs(np.fft.fft(Vrec) / np.fft.fft(Vinc))
|
||||
|
||||
# Calculate input impedance
|
||||
with np.errstate(divide='ignore'):
|
||||
zin = (np.fft.fft(Vtotal) * delaycorrection) / np.fft.fft(Itotal)
|
||||
|
||||
# Calculate input admittance
|
||||
with np.errstate(divide='ignore'):
|
||||
yin = np.fft.fft(Itotal) / (np.fft.fft(Vtotal) * delaycorrection)
|
||||
|
||||
# Convert to decibels (ignore warning from taking a log of any zero values)
|
||||
with np.errstate(divide='ignore'):
|
||||
Vincp = 20 * np.log10(np.abs((np.fft.fft(Vinc) * delaycorrection)))
|
||||
Iincp = 20 * np.log10(np.abs(np.fft.fft(Iinc)))
|
||||
Vrefp = 20 * np.log10(np.abs((np.fft.fft(Vref) * delaycorrection)))
|
||||
Irefp = 20 * np.log10(np.abs(np.fft.fft(Iref)))
|
||||
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)
|
||||
|
||||
# Replace any NaNs or Infs from zero division
|
||||
Vincp[np.invert(np.isfinite(Vincp))] = 0
|
||||
Iincp[np.invert(np.isfinite(Iincp))] = 0
|
||||
Vrefp[np.invert(np.isfinite(Vrefp))] = 0
|
||||
Irefp[np.invert(np.isfinite(Irefp))] = 0
|
||||
Vtotalp[np.invert(np.isfinite(Vtotalp))] = 0
|
||||
Itotalp[np.invert(np.isfinite(Itotalp))] = 0
|
||||
s11[np.invert(np.isfinite(s11))] = 0
|
||||
|
||||
# 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}
|
||||
if tlrxnumber or rxnumber:
|
||||
with np.errstate(divide='ignore'): # Ignore warning from taking a log of any zero values
|
||||
s21 = 20 * np.log10(s21)
|
||||
s21[np.invert(np.isfinite(s21))] = 0
|
||||
antennaparams['s21'] = s21
|
||||
|
||||
return antennaparams
|
||||
|
||||
|
||||
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.
|
||||
freq (array): Frequencies for FFTs.
|
||||
Vinc, Vincp, Iinc, Iincp (array): Time and frequency domain representations of incident voltage and current.
|
||||
Vref, Vrefp, Iref, Irefp (array): Time and frequency domain representations of reflected voltage and current.
|
||||
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
|
||||
pltrangemax = np.where((np.amax(Vincp[1::]) - Vincp[1::]) > 60)[0][0] + 1
|
||||
# To a maximum frequency
|
||||
# pltrangemax = np.where(freqs > 6e9)[0][0]
|
||||
pltrange = np.s_[pltrangemin:pltrangemax]
|
||||
|
||||
# Print some useful values from s11, and input impedance
|
||||
s11minfreq = np.where(s11[pltrange] == np.amin(s11[pltrange]))[0][0]
|
||||
print('s11 minimum: {:g} dB at {:g} Hz'.format(np.amin(s11[pltrange]), freqs[s11minfreq + pltrangemin]))
|
||||
print('At {:g} Hz...'.format(freqs[s11minfreq + pltrangemin]))
|
||||
print('Input impedance: {:.1f}{:+.1f}j Ohms'.format(np.abs(zin[s11minfreq + pltrangemin]), zin[s11minfreq + pltrangemin].imag))
|
||||
# print('Input admittance (mag): {:g} S'.format(np.abs(yin[s11minfreq + pltrangemin])))
|
||||
# print('Input admittance (phase): {:.1f} deg'.format(np.angle(yin[s11minfreq + pltrangemin], deg=True)))
|
||||
|
||||
# Figure 1
|
||||
# Plot incident voltage
|
||||
fig1, ax = plt.subplots(num='Transmitter transmission line parameters', figsize=(20, 12), facecolor='w', edgecolor='w')
|
||||
gs1 = gridspec.GridSpec(4, 2, hspace=0.7)
|
||||
ax = plt.subplot(gs1[0, 0])
|
||||
ax.plot(time, Vinc, 'r', lw=2, label='Vinc')
|
||||
ax.set_title('Incident voltage')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Voltage [V]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of incident voltage
|
||||
ax = plt.subplot(gs1[0, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Vincp[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax.plot(freqs[pltrange], Vincp[pltrange], 'r', lw=2)
|
||||
ax.set_title('Incident voltage')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot incident current
|
||||
ax = plt.subplot(gs1[1, 0])
|
||||
ax.plot(time, Iinc, 'b', lw=2, label='Vinc')
|
||||
ax.set_title('Incident current')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Current [A]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of incident current
|
||||
ax = plt.subplot(gs1[1, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Iincp[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'b')
|
||||
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
ax.plot(freqs[pltrange], Iincp[pltrange], 'b', lw=2)
|
||||
ax.set_title('Incident current')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot total voltage
|
||||
ax = plt.subplot(gs1[2, 0])
|
||||
ax.plot(time, Vtotal, 'r', lw=2, label='Vinc')
|
||||
ax.set_title('Total (incident + reflected) voltage')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Voltage [V]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of total voltage
|
||||
ax = plt.subplot(gs1[2, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Vtotalp[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax.plot(freqs[pltrange], Vtotalp[pltrange], 'r', lw=2)
|
||||
ax.set_title('Total (incident + reflected) voltage')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot total current
|
||||
ax = plt.subplot(gs1[3, 0])
|
||||
ax.plot(time, Itotal, 'b', lw=2, label='Vinc')
|
||||
ax.set_title('Total (incident + reflected) current')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Current [A]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of total current
|
||||
ax = plt.subplot(gs1[3, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Itotalp[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'b')
|
||||
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
ax.plot(freqs[pltrange], Itotalp[pltrange], 'b', lw=2)
|
||||
ax.set_title('Total (incident + reflected) current')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of reflected voltage
|
||||
# ax = plt.subplot(gs1[4, 1])
|
||||
# markerline, stemlines, baseline = ax.stem(freqs[pltrange], Vrefp[pltrange], '-.', use_line_collection=True)
|
||||
# plt.setp(baseline, 'linewidth', 0)
|
||||
# plt.setp(stemlines, 'color', 'r')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
# ax.plot(freqs[pltrange], Vrefp[pltrange], 'r', lw=2)
|
||||
# ax.set_title('Reflected voltage')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Power [dB]')
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of reflected current
|
||||
# ax = plt.subplot(gs1[5, 1])
|
||||
# markerline, stemlines, baseline = ax.stem(freqs[pltrange], Irefp[pltrange], '-.', use_line_collection=True)
|
||||
# plt.setp(baseline, 'linewidth', 0)
|
||||
# plt.setp(stemlines, 'color', 'b')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
# ax.plot(freqs[pltrange], Irefp[pltrange], 'b', lw=2)
|
||||
# ax.set_title('Reflected current')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Power [dB]')
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Figure 2
|
||||
# Plot frequency spectra of s11
|
||||
fig2, ax = plt.subplots(num='Antenna parameters', figsize=(20, 12), facecolor='w', edgecolor='w')
|
||||
gs2 = gridspec.GridSpec(2, 2, hspace=0.3)
|
||||
ax = plt.subplot(gs2[0, 0])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], s11[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], s11[pltrange], 'g', lw=2)
|
||||
ax.set_title('s11')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
# ax.set_xlim([0, 5e9])
|
||||
# ax.set_ylim([-25, 0])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of s21
|
||||
if s21 is not None:
|
||||
ax = plt.subplot(gs2[0, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], s21[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], s21[pltrange], 'g', lw=2)
|
||||
ax.set_title('s21')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([-25, 50])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot input resistance (real part of impedance)
|
||||
ax = plt.subplot(gs2[1, 0])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], zin[pltrange].real, '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], zin[pltrange].real, 'g', lw=2)
|
||||
ax.set_title('Input impedance (resistive)')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Resistance [Ohms]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
ax.set_ylim(bottom=0)
|
||||
# ax.set_ylim([0, 300])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot input reactance (imaginery part of impedance)
|
||||
ax = plt.subplot(gs2[1, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], zin[pltrange].imag, '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], zin[pltrange].imag, 'g', lw=2)
|
||||
ax.set_title('Input impedance (reactive)')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Reactance [Ohms]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([-300, 300])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot input admittance (magnitude)
|
||||
# ax = plt.subplot(gs2[2, 0])
|
||||
# markerline, stemlines, baseline = ax.stem(freqs[pltrange], np.abs(yin[pltrange]), '-.', use_line_collection=True)
|
||||
# plt.setp(baseline, 'linewidth', 0)
|
||||
# plt.setp(stemlines, 'color', 'g')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
# ax.plot(freqs[pltrange], np.abs(yin[pltrange]), 'g', lw=2)
|
||||
# ax.set_title('Input admittance (magnitude)')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Admittance [Siemens]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([0, 0.035])
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot input admittance (phase)
|
||||
# ax = plt.subplot(gs2[2, 1])
|
||||
# markerline, stemlines, baseline = ax.stem(freqs[pltrange], np.angle(yin[pltrange], deg=True), '-.', use_line_collection=True)
|
||||
# plt.setp(baseline, 'linewidth', 0)
|
||||
# plt.setp(stemlines, 'color', 'g')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
# ax.plot(freqs[pltrange], np.angle(yin[pltrange], deg=True), 'g', lw=2)
|
||||
# ax.set_title('Input admittance (phase)')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Phase [degrees]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([-40, 100])
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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)
|
||||
# fig2.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_ant_params.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
# fig1.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_tl_params.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
# fig2.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_ant_params.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
|
||||
return plt
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(description='Plots antenna parameters (s11, s21 parameters and input impedance) from an output file containing a transmission line source.', usage='cd gprMax; python -m tools.plot_antenna_params outputfile')
|
||||
parser.add_argument('outputfile', help='name of output file including path')
|
||||
parser.add_argument('--tltx-num', default=1, type=int, help='transmitter antenna - transmission line number')
|
||||
parser.add_argument('--tlrx-num', type=int, help='receiver antenna - transmission line number')
|
||||
parser.add_argument('--rx-num', type=int, help='receiver antenna - output number')
|
||||
parser.add_argument('--rx-component', type=str, help='receiver antenna - output electric field component', choices=['Ex', 'Ey', 'Ez'])
|
||||
args = parser.parse_args()
|
||||
|
||||
antennaparams = calculate_antenna_params(args.outputfile, args.tltx_num, args.tlrx_num, args.rx_num, args.rx_component)
|
||||
plthandle = mpl_plot(args.outputfile, **antennaparams)
|
||||
plthandle.show()
|
||||
# Copyright (C) 2015-2023: The University of Edinburgh
|
||||
# Authors: Craig Warren and Antonis Giannopoulos
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# 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/>.
|
||||
|
||||
import argparse
|
||||
import os
|
||||
import sys
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
|
||||
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.
|
||||
"""
|
||||
|
||||
# Open output file and read some attributes
|
||||
f = h5py.File(filename, 'r')
|
||||
dxdydz = f.attrs['dx_dy_dz']
|
||||
dt = f.attrs['dt']
|
||||
iterations = f.attrs['Iterations']
|
||||
|
||||
# Calculate time array and frequency bin spacing
|
||||
time = np.linspace(0, (iterations - 1) * dt, num=iterations)
|
||||
df = 1 / np.amax(time)
|
||||
|
||||
print('Time window: {:g} s ({} iterations)'.format(np.amax(time), iterations))
|
||||
print('Time step: {:g} s'.format(dt))
|
||||
print('Frequency bin spacing: {:g} Hz'.format(df))
|
||||
|
||||
# Read/calculate voltages and currents from transmitter antenna
|
||||
tltxpath = '/tls/tl' + str(tltxnumber) + '/'
|
||||
|
||||
# Incident voltages/currents
|
||||
Vinc = f[tltxpath + 'Vinc'][:]
|
||||
Iinc = f[tltxpath + 'Iinc'][:]
|
||||
|
||||
# Total (incident + reflected) voltages/currents
|
||||
Vtotal = f[tltxpath + 'Vtotal'][:]
|
||||
Itotal = f[tltxpath + 'Itotal'][:]
|
||||
|
||||
# Reflected voltages/currents
|
||||
Vref = Vtotal - Vinc
|
||||
Iref = Itotal - Iinc
|
||||
|
||||
# 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'][:]
|
||||
|
||||
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]
|
||||
elif rxcomponent == 'Ey':
|
||||
Vrec = f[rxpath][:] * -1 * dxdydz[1]
|
||||
elif rxcomponent == 'Ez':
|
||||
Vrec = f[rxpath][:] * -1 * dxdydz[2]
|
||||
f.close()
|
||||
|
||||
# Frequency bins
|
||||
freqs = np.fft.fftfreq(Vinc.size, d=dt)
|
||||
|
||||
# Delay correction - current lags voltage, so delay voltage to match current timestep
|
||||
delaycorrection = np.exp(1j * 2 * np.pi * freqs * (dt / 2))
|
||||
|
||||
# Calculate s11 and (optionally) s21
|
||||
with np.errstate(divide='ignore'):
|
||||
s11 = np.abs(np.fft.fft(Vref) / np.fft.fft(Vinc))
|
||||
if tlrxnumber or rxnumber:
|
||||
with np.errstate(divide='ignore'):
|
||||
s21 = np.abs(np.fft.fft(Vrec) / np.fft.fft(Vinc))
|
||||
|
||||
# Calculate input impedance
|
||||
with np.errstate(divide='ignore'):
|
||||
zin = (np.fft.fft(Vtotal) * delaycorrection) / np.fft.fft(Itotal)
|
||||
|
||||
# Calculate input admittance
|
||||
with np.errstate(divide='ignore'):
|
||||
yin = np.fft.fft(Itotal) / (np.fft.fft(Vtotal) * delaycorrection)
|
||||
|
||||
# Convert to decibels (ignore warning from taking a log of any zero values)
|
||||
with np.errstate(divide='ignore'):
|
||||
Vincp = 20 * np.log10(np.abs((np.fft.fft(Vinc) * delaycorrection)))
|
||||
Iincp = 20 * np.log10(np.abs(np.fft.fft(Iinc)))
|
||||
Vrefp = 20 * np.log10(np.abs((np.fft.fft(Vref) * delaycorrection)))
|
||||
Irefp = 20 * np.log10(np.abs(np.fft.fft(Iref)))
|
||||
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)
|
||||
|
||||
# Replace any NaNs or Infs from zero division
|
||||
Vincp[np.invert(np.isfinite(Vincp))] = 0
|
||||
Iincp[np.invert(np.isfinite(Iincp))] = 0
|
||||
Vrefp[np.invert(np.isfinite(Vrefp))] = 0
|
||||
Irefp[np.invert(np.isfinite(Irefp))] = 0
|
||||
Vtotalp[np.invert(np.isfinite(Vtotalp))] = 0
|
||||
Itotalp[np.invert(np.isfinite(Itotalp))] = 0
|
||||
s11[np.invert(np.isfinite(s11))] = 0
|
||||
|
||||
# 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}
|
||||
if tlrxnumber or rxnumber:
|
||||
with np.errstate(divide='ignore'): # Ignore warning from taking a log of any zero values
|
||||
s21 = 20 * np.log10(s21)
|
||||
s21[np.invert(np.isfinite(s21))] = 0
|
||||
antennaparams['s21'] = s21
|
||||
|
||||
return antennaparams
|
||||
|
||||
|
||||
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.
|
||||
freq (array): Frequencies for FFTs.
|
||||
Vinc, Vincp, Iinc, Iincp (array): Time and frequency domain representations of incident voltage and current.
|
||||
Vref, Vrefp, Iref, Irefp (array): Time and frequency domain representations of reflected voltage and current.
|
||||
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
|
||||
pltrangemax = np.where((np.amax(Vincp[1::]) - Vincp[1::]) > 60)[0][0] + 1
|
||||
# To a maximum frequency
|
||||
# pltrangemax = np.where(freqs > 6e9)[0][0]
|
||||
pltrange = np.s_[pltrangemin:pltrangemax]
|
||||
|
||||
# Print some useful values from s11, and input impedance
|
||||
s11minfreq = np.where(s11[pltrange] == np.amin(s11[pltrange]))[0][0]
|
||||
print('s11 minimum: {:g} dB at {:g} Hz'.format(np.amin(s11[pltrange]), freqs[s11minfreq + pltrangemin]))
|
||||
print('At {:g} Hz...'.format(freqs[s11minfreq + pltrangemin]))
|
||||
print('Input impedance: {:.1f}{:+.1f}j Ohms'.format(np.abs(zin[s11minfreq + pltrangemin]), zin[s11minfreq + pltrangemin].imag))
|
||||
# print('Input admittance (mag): {:g} S'.format(np.abs(yin[s11minfreq + pltrangemin])))
|
||||
# print('Input admittance (phase): {:.1f} deg'.format(np.angle(yin[s11minfreq + pltrangemin], deg=True)))
|
||||
|
||||
# Figure 1
|
||||
# Plot incident voltage
|
||||
fig1, ax = plt.subplots(num='Transmitter transmission line parameters', figsize=(20, 12), facecolor='w', edgecolor='w')
|
||||
gs1 = gridspec.GridSpec(4, 2, hspace=0.7)
|
||||
ax = plt.subplot(gs1[0, 0])
|
||||
ax.plot(time, Vinc, 'r', lw=2, label='Vinc')
|
||||
ax.set_title('Incident voltage')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Voltage [V]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of incident voltage
|
||||
ax = plt.subplot(gs1[0, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Vincp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax.plot(freqs[pltrange], Vincp[pltrange], 'r', lw=2)
|
||||
ax.set_title('Incident voltage')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot incident current
|
||||
ax = plt.subplot(gs1[1, 0])
|
||||
ax.plot(time, Iinc, 'b', lw=2, label='Vinc')
|
||||
ax.set_title('Incident current')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Current [A]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of incident current
|
||||
ax = plt.subplot(gs1[1, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Iincp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'b')
|
||||
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
ax.plot(freqs[pltrange], Iincp[pltrange], 'b', lw=2)
|
||||
ax.set_title('Incident current')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot total voltage
|
||||
ax = plt.subplot(gs1[2, 0])
|
||||
ax.plot(time, Vtotal, 'r', lw=2, label='Vinc')
|
||||
ax.set_title('Total (incident + reflected) voltage')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Voltage [V]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of total voltage
|
||||
ax = plt.subplot(gs1[2, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Vtotalp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax.plot(freqs[pltrange], Vtotalp[pltrange], 'r', lw=2)
|
||||
ax.set_title('Total (incident + reflected) voltage')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot total current
|
||||
ax = plt.subplot(gs1[3, 0])
|
||||
ax.plot(time, Itotal, 'b', lw=2, label='Vinc')
|
||||
ax.set_title('Total (incident + reflected) current')
|
||||
ax.set_xlabel('Time [s]')
|
||||
ax.set_ylabel('Current [A]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of total current
|
||||
ax = plt.subplot(gs1[3, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], Itotalp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'b')
|
||||
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
ax.plot(freqs[pltrange], Itotalp[pltrange], 'b', lw=2)
|
||||
ax.set_title('Total (incident + reflected) current')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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)
|
||||
# plt.setp(stemlines, 'color', 'r')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
# ax.plot(freqs[pltrange], Vrefp[pltrange], 'r', lw=2)
|
||||
# ax.set_title('Reflected voltage')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Power [dB]')
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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)
|
||||
# plt.setp(stemlines, 'color', 'b')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
# ax.plot(freqs[pltrange], Irefp[pltrange], 'b', lw=2)
|
||||
# ax.set_title('Reflected current')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Power [dB]')
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Figure 2
|
||||
# Plot frequency spectra of s11
|
||||
fig2, ax = plt.subplots(num='Antenna parameters', figsize=(20, 12), facecolor='w', edgecolor='w')
|
||||
gs2 = gridspec.GridSpec(2, 2, hspace=0.3)
|
||||
ax = plt.subplot(gs2[0, 0])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], s11[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], s11[pltrange], 'g', lw=2)
|
||||
ax.set_title('s11')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
# ax.set_xlim([0, 5e9])
|
||||
# ax.set_ylim([-25, 0])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot frequency spectra of s21
|
||||
if s21 is not None:
|
||||
ax = plt.subplot(gs2[0, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], s21[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], s21[pltrange], 'g', lw=2)
|
||||
ax.set_title('s21')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([-25, 50])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot input resistance (real part of impedance)
|
||||
ax = plt.subplot(gs2[1, 0])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], zin[pltrange].real, '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], zin[pltrange].real, 'g', lw=2)
|
||||
ax.set_title('Input impedance (resistive)')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Resistance [Ohms]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
ax.set_ylim(bottom=0)
|
||||
# ax.set_ylim([0, 300])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# Plot input reactance (imaginery part of impedance)
|
||||
ax = plt.subplot(gs2[1, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange], zin[pltrange].imag, '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange], zin[pltrange].imag, 'g', lw=2)
|
||||
ax.set_title('Input impedance (reactive)')
|
||||
ax.set_xlabel('Frequency [Hz]')
|
||||
ax.set_ylabel('Reactance [Ohms]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([-300, 300])
|
||||
ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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)
|
||||
# plt.setp(stemlines, 'color', 'g')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
# ax.plot(freqs[pltrange], np.abs(yin[pltrange]), 'g', lw=2)
|
||||
# ax.set_title('Input admittance (magnitude)')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Admittance [Siemens]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([0, 0.035])
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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)
|
||||
# plt.setp(stemlines, 'color', 'g')
|
||||
# plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
# ax.plot(freqs[pltrange], np.angle(yin[pltrange], deg=True), 'g', lw=2)
|
||||
# ax.set_title('Input admittance (phase)')
|
||||
# ax.set_xlabel('Frequency [Hz]')
|
||||
# ax.set_ylabel('Phase [degrees]')
|
||||
# ax.set_xlim([0.88e9, 1.02e9])
|
||||
# ax.set_ylim([-40, 100])
|
||||
# ax.grid(which='both', axis='both', linestyle='-.')
|
||||
|
||||
# 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)
|
||||
# fig2.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_ant_params.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
# fig1.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_tl_params.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
# fig2.savefig(os.path.splitext(os.path.abspath(filename))[0] + '_ant_params.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
|
||||
return plt
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(description='Plots antenna parameters (s11, s21 parameters and input impedance) from an output file containing a transmission line source.', usage='cd gprMax; python -m tools.plot_antenna_params outputfile')
|
||||
parser.add_argument('outputfile', help='name of output file including path')
|
||||
parser.add_argument('--tltx-num', default=1, type=int, help='transmitter antenna - transmission line number')
|
||||
parser.add_argument('--tlrx-num', type=int, help='receiver antenna - transmission line number')
|
||||
parser.add_argument('--rx-num', type=int, help='receiver antenna - output number')
|
||||
parser.add_argument('--rx-component', type=str, help='receiver antenna - output electric field component', choices=['Ex', 'Ey', 'Ez'])
|
||||
args = parser.parse_args()
|
||||
|
||||
antennaparams = calculate_antenna_params(args.outputfile, args.tltx_num, args.tlrx_num, args.rx_num, args.rx_component)
|
||||
plthandle = mpl_plot(args.outputfile, **antennaparams)
|
||||
plthandle.show()
|
||||
|
||||
@@ -1,173 +1,173 @@
|
||||
# Copyright (C) 2015-2023: The University of Edinburgh
|
||||
# Authors: Craig Warren and Antonis Giannopoulos
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# 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/>.
|
||||
|
||||
import argparse
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from gprMax.exceptions import CmdInputError
|
||||
from gprMax.utilities import fft_power
|
||||
from gprMax.utilities import round_value
|
||||
from gprMax.waveforms import Waveform
|
||||
|
||||
|
||||
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
|
||||
timewindow = (timewindow - 1) * dt
|
||||
|
||||
except:
|
||||
timewindow = float(timewindow)
|
||||
if timewindow > 0:
|
||||
iterations = round_value((timewindow / dt)) + 1
|
||||
else:
|
||||
raise CmdInputError('Time window must have a value greater than zero')
|
||||
|
||||
return timewindow, iterations
|
||||
|
||||
|
||||
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, (iterations - 1) * dt, num=iterations)
|
||||
waveform = np.zeros(len(time))
|
||||
timeiter = np.nditer(time, flags=['c_index'])
|
||||
|
||||
while not timeiter.finished:
|
||||
waveform[timeiter.index] = w.calculate_value(timeiter[0], dt)
|
||||
timeiter.iternext()
|
||||
|
||||
print('Waveform characteristics...')
|
||||
print('Type: {}'.format(w.type))
|
||||
print('Maximum (absolute) amplitude: {:g}'.format(np.max(np.abs(waveform))))
|
||||
|
||||
if w.freq and not w.type == 'gaussian':
|
||||
print('Centre frequency: {:g} Hz'.format(w.freq))
|
||||
|
||||
if w.type == 'gaussian' or w.type == 'gaussiandot' or w.type == 'gaussiandotnorm' or w.type == 'gaussianprime' or w.type == 'gaussiandoubleprime':
|
||||
delay = 1 / w.freq
|
||||
print('Time to centre of pulse: {:g} s'.format(delay))
|
||||
elif w.type == 'gaussiandotdot' or w.type == 'gaussiandotdotnorm' or w.type == 'ricker':
|
||||
delay = np.sqrt(2) / w.freq
|
||||
print('Time to centre of pulse: {:g} s'.format(delay))
|
||||
|
||||
print('Time window: {:g} s ({} iterations)'.format(timewindow, iterations))
|
||||
print('Time step: {:g} s'.format(dt))
|
||||
|
||||
if fft:
|
||||
# FFT
|
||||
freqs, power = fft_power(waveform, dt)
|
||||
|
||||
# Set plotting range to 4 times frequency at max power of waveform or
|
||||
# 4 times the centre frequency
|
||||
freqmaxpower = np.where(np.isclose(power, 0))[0][0]
|
||||
if freqs[freqmaxpower] > w.freq:
|
||||
pltrange = np.where(freqs > 4 * freqs[freqmaxpower])[0][0]
|
||||
else:
|
||||
pltrange = np.where(freqs > 4 * w.freq)[0][0]
|
||||
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]')
|
||||
ax1.set_ylabel('Amplitude')
|
||||
|
||||
# Plot frequency spectra
|
||||
markerline, stemlines, baseline = ax2.stem(freqs[pltrange], power[pltrange], '-.', use_line_collection=True)
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2)
|
||||
ax2.set_xlabel('Frequency [Hz]')
|
||||
ax2.set_ylabel('Power [dB]')
|
||||
|
||||
else:
|
||||
fig, ax1 = plt.subplots(num=w.type, figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
|
||||
# Plot waveform
|
||||
ax1.plot(time, waveform, 'r', lw=2)
|
||||
ax1.set_xlabel('Time [s]')
|
||||
ax1.set_ylabel('Amplitude')
|
||||
|
||||
[ax.grid(which='both', axis='both', linestyle='-.') 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)
|
||||
# fig.savefig(os.path.dirname(os.path.abspath(__file__)) + os.sep + w.type + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
|
||||
return plt
|
||||
|
||||
|
||||
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)
|
||||
parser.add_argument('amp', type=float, help='amplitude of waveform')
|
||||
parser.add_argument('freq', type=float, help='centre frequency of waveform')
|
||||
parser.add_argument('timewindow', help='time window to view waveform')
|
||||
parser.add_argument('dt', type=float, help='time step to view waveform')
|
||||
parser.add_argument('-fft', action='store_true', help='plot FFT of waveform', default=False)
|
||||
args = parser.parse_args()
|
||||
|
||||
# Check waveform parameters
|
||||
if args.type.lower() not in Waveform.types:
|
||||
raise CmdInputError('The waveform must have one of the following types {}'.format(', '.join(Waveform.types)))
|
||||
if args.freq <= 0:
|
||||
raise CmdInputError('The waveform requires an excitation frequency value of greater than zero')
|
||||
|
||||
# Create waveform instance
|
||||
w = Waveform()
|
||||
w.type = args.type
|
||||
w.amp = args.amp
|
||||
w.freq = args.freq
|
||||
|
||||
timewindow, iterations = check_timewindow(args.timewindow, args.dt)
|
||||
plthandle = mpl_plot(w, timewindow, args.dt, iterations, args.fft)
|
||||
plthandle.show()
|
||||
# Copyright (C) 2015-2023: The University of Edinburgh
|
||||
# Authors: Craig Warren and Antonis Giannopoulos
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# 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/>.
|
||||
|
||||
import argparse
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from gprMax.exceptions import CmdInputError
|
||||
from gprMax.utilities import fft_power
|
||||
from gprMax.utilities import round_value
|
||||
from gprMax.waveforms import Waveform
|
||||
|
||||
|
||||
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
|
||||
timewindow = (timewindow - 1) * dt
|
||||
|
||||
except:
|
||||
timewindow = float(timewindow)
|
||||
if timewindow > 0:
|
||||
iterations = round_value((timewindow / dt)) + 1
|
||||
else:
|
||||
raise CmdInputError('Time window must have a value greater than zero')
|
||||
|
||||
return timewindow, iterations
|
||||
|
||||
|
||||
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, (iterations - 1) * dt, num=iterations)
|
||||
waveform = np.zeros(len(time))
|
||||
timeiter = np.nditer(time, flags=['c_index'])
|
||||
|
||||
while not timeiter.finished:
|
||||
waveform[timeiter.index] = w.calculate_value(timeiter[0], dt)
|
||||
timeiter.iternext()
|
||||
|
||||
print('Waveform characteristics...')
|
||||
print('Type: {}'.format(w.type))
|
||||
print('Maximum (absolute) amplitude: {:g}'.format(np.max(np.abs(waveform))))
|
||||
|
||||
if w.freq and not w.type == 'gaussian':
|
||||
print('Centre frequency: {:g} Hz'.format(w.freq))
|
||||
|
||||
if w.type == 'gaussian' or w.type == 'gaussiandot' or w.type == 'gaussiandotnorm' or w.type == 'gaussianprime' or w.type == 'gaussiandoubleprime':
|
||||
delay = 1 / w.freq
|
||||
print('Time to centre of pulse: {:g} s'.format(delay))
|
||||
elif w.type == 'gaussiandotdot' or w.type == 'gaussiandotdotnorm' or w.type == 'ricker':
|
||||
delay = np.sqrt(2) / w.freq
|
||||
print('Time to centre of pulse: {:g} s'.format(delay))
|
||||
|
||||
print('Time window: {:g} s ({} iterations)'.format(timewindow, iterations))
|
||||
print('Time step: {:g} s'.format(dt))
|
||||
|
||||
if fft:
|
||||
# FFT
|
||||
freqs, power = fft_power(waveform, dt)
|
||||
|
||||
# Set plotting range to 4 times frequency at max power of waveform or
|
||||
# 4 times the centre frequency
|
||||
freqmaxpower = np.where(np.isclose(power, 0))[0][0]
|
||||
if freqs[freqmaxpower] > w.freq:
|
||||
pltrange = np.where(freqs > 4 * freqs[freqmaxpower])[0][0]
|
||||
else:
|
||||
pltrange = np.where(freqs > 4 * w.freq)[0][0]
|
||||
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]')
|
||||
ax1.set_ylabel('Amplitude')
|
||||
|
||||
# Plot frequency spectra
|
||||
markerline, stemlines, baseline = ax2.stem(freqs[pltrange], power[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2)
|
||||
ax2.set_xlabel('Frequency [Hz]')
|
||||
ax2.set_ylabel('Power [dB]')
|
||||
|
||||
else:
|
||||
fig, ax1 = plt.subplots(num=w.type, figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
|
||||
# Plot waveform
|
||||
ax1.plot(time, waveform, 'r', lw=2)
|
||||
ax1.set_xlabel('Time [s]')
|
||||
ax1.set_ylabel('Amplitude')
|
||||
|
||||
[ax.grid(which='both', axis='both', linestyle='-.') 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)
|
||||
# fig.savefig(os.path.dirname(os.path.abspath(__file__)) + os.sep + w.type + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
|
||||
return plt
|
||||
|
||||
|
||||
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)
|
||||
parser.add_argument('amp', type=float, help='amplitude of waveform')
|
||||
parser.add_argument('freq', type=float, help='centre frequency of waveform')
|
||||
parser.add_argument('timewindow', help='time window to view waveform')
|
||||
parser.add_argument('dt', type=float, help='time step to view waveform')
|
||||
parser.add_argument('-fft', action='store_true', help='plot FFT of waveform', default=False)
|
||||
args = parser.parse_args()
|
||||
|
||||
# Check waveform parameters
|
||||
if args.type.lower() not in Waveform.types:
|
||||
raise CmdInputError('The waveform must have one of the following types {}'.format(', '.join(Waveform.types)))
|
||||
if args.freq <= 0:
|
||||
raise CmdInputError('The waveform requires an excitation frequency value of greater than zero')
|
||||
|
||||
# Create waveform instance
|
||||
w = Waveform()
|
||||
w.type = args.type
|
||||
w.amp = args.amp
|
||||
w.freq = args.freq
|
||||
|
||||
timewindow, iterations = check_timewindow(args.timewindow, args.dt)
|
||||
plthandle = mpl_plot(w, timewindow, args.dt, iterations, args.fft)
|
||||
plthandle.show()
|
||||
|
||||
在新工单中引用
屏蔽一个用户