From cb5d26b4336b702f317835f44d2381f9dcb7a06c Mon Sep 17 00:00:00 2001
From: JingziC <140467318+JingziC@users.noreply.github.com>
Date: Mon, 29 Jan 2024 17:10:01 +0900
Subject: [PATCH] Removing use_line_collection in plotting
use_line_collection is removed in matplotlib 3.8.0
---
tools/plot_Ascan.py | 474 ++++++++++----------
tools/plot_antenna_params.py | 848 +++++++++++++++++------------------
tools/plot_source_wave.py | 346 +++++++-------
3 files changed, 834 insertions(+), 834 deletions(-)
diff --git a/tools/plot_Ascan.py b/tools/plot_Ascan.py
index 7c5d2f29..9b1ffed1 100644
--- a/tools/plot_Ascan.py
+++ b/tools/plot_Ascan.py
@@ -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 .
-
-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 .
+
+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()
diff --git a/tools/plot_antenna_params.py b/tools/plot_antenna_params.py
index b54b1710..dffa6155 100644
--- a/tools/plot_antenna_params.py
+++ b/tools/plot_antenna_params.py
@@ -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 .
-
-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 .
+
+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()
diff --git a/tools/plot_source_wave.py b/tools/plot_source_wave.py
index 88710fd1..4f8abb0a 100644
--- a/tools/plot_source_wave.py
+++ b/tools/plot_source_wave.py
@@ -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 .
-
-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 .
+
+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()