From 4689a323694b319a2297e3e675feff038ab0ff21 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 6 Jan 2022 16:14:58 +0000 Subject: [PATCH] Removed this tool. --- tools/plot_antenna_params_UtEl.py | 432 ------------------------------ 1 file changed, 432 deletions(-) delete mode 100644 tools/plot_antenna_params_UtEl.py diff --git a/tools/plot_antenna_params_UtEl.py b/tools/plot_antenna_params_UtEl.py deleted file mode 100644 index 23801f0a..00000000 --- a/tools/plot_antenna_params_UtEl.py +++ /dev/null @@ -1,432 +0,0 @@ -# Copyright (C) 2015-2022: The University of Edinburgh -# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley -# -# 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 logging -from pathlib import Path - -import h5py -import matplotlib.gridspec as gridspec -import matplotlib.pyplot as plt -import numpy as np -from scipy.io import loadmat - -logger = logging.getLogger(__name__) - - -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 - file = Path(filename) - f = h5py.File(file, '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(f'Time window: {np.amax(time):g} s ({iterations} iterations)') - print(f'Time step: {dt:g} s') - print(f'Frequency bin spacing: {df:g} Hz') - - # 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: - logger.exception(f"{rxcomponent} output requested, but the available output for receiver {rxnumber} is {', '.join(availableoutputs)}") - raise ValueError - - 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 - s11 = np.abs(np.fft.fft(Vref) / np.fft.fft(Vinc)) - if tlrxnumber or rxnumber: - s21 = np.abs(np.fft.fft(Vrec) / np.fft.fft(Vinc)) - - # Calculate input impedance - zin = (np.fft.fft(Vtotal) * delaycorrection) / np.fft.fft(Itotal) - - # Calculate input admittance - 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 > 3e9)[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(f's11 minimum: {np.amin(s11[pltrange]):g} dB at {freqs[s11minfreq + pltrangemin]:g} Hz') - print(f'At {freqs[s11minfreq + pltrangemin]:g} Hz...') - print(f'Input impedance: {np.abs(zin[s11minfreq + pltrangemin]):.1f}{zin[s11minfreq + pltrangemin].imag:+.1f}j Ohms') - # print(f'Input admittance (mag): {np.abs(yin[s11minfreq + pltrangemin]):g} S') - # print(f'Input admittance (phase): {np.angle(yin[s11minfreq + pltrangemin], deg=True):.1f} deg') - - # 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 - # Load measured s11 - s11_meas = loadmat('/Volumes/Users/cwarren/Dropbox (Northumbria University)/Research/GPR-antennas/Utsi Electronics/500MHz_Vivaldi_antenna/measured/s11_measured.mat') - - # Plot frequency spectra of s11 - fig2, ax = plt.subplots(num='Antenna parameters', figsize=(20, 12), facecolor='w', edgecolor='w') - gs2 = gridspec.GridSpec(1, 1, hspace=0.3) - ax = plt.subplot(gs2[0, 0]) - markerline, stemlines, baseline = ax.stem(freqs[pltrange], s11[pltrange], '-.', use_line_collection=True) - markerline2, stemlines2, baseline2 = ax.stem(s11_meas['f'], s11_meas['s11_t'], '-.', use_line_collection=True) - plt.setp(baseline, 'linewidth', 0) - plt.setp(stemlines, 'color', 'w') - plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g') - plt.setp(baseline2, 'linewidth', 0) - plt.setp(stemlines2, 'color', 'w') - plt.setp(markerline2, 'markerfacecolor', 'b', 'markeredgecolor', 'b') - ax.plot(freqs[pltrange], s11[pltrange], 'g', lw=2, label='Simulation (Debye, 50$\Omega$)') - ax.plot(s11_meas['f'], s11_meas['s11_t'], 'b', lw=2, label='Measured') - ax.set_title('s11') - ax.set_xlabel('Frequency [Hz]') - ax.set_ylabel('Power [dB]') - ax.set_xlim([0, 3e9]) - ax.set_ylim([-40, 0]) - ax.grid(which='both', axis='both', linestyle='-.') - ax.legend(fontsize=12, frameon=False) - - # 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(filename.with_suffix('.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()