Updated style for FFT plots; added command line argument to specify transmission line number.

tools/plot_s11.py

@@ -21,6 +21,7 @@ import h5py
 import numpy as np
 np.seterr(divide='ignore', invalid='ignore')
 import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
 
 from gprMax.exceptions import CmdInputError
 
@@ -29,6 +30,7 @@ from gprMax.exceptions import CmdInputError
 # Parse command line arguments
 parser = argparse.ArgumentParser(description='Plots the s11 scattering parameter (input port voltage reflection coefficient) from an output file containing a transmission line source.', usage='cd gprMax; python -m tools.plot_s11 outputfile')
 parser.add_argument('outputfile', help='name of output file including path')
+parser.add_argument('-tln', default=1, type=int, help='transmission line number')
 args = parser.parse_args()
 
 # Open output file and read some attributes
@@ -39,7 +41,7 @@ iterations = f.attrs['Iterations']
 time = np.arange(0, dt * iterations, dt)
 time = time / 1e-9
 
-path = '/tls/tl1/'
+path = '/tls/tl' + str(args.tln) + '/'
 Vinc = f[path + 'Vinc'][:]
 Vscat = f[path + 'Vscat'][:]
 Vtotal = f[path +'Vtotal'][:]
@@ -48,34 +50,40 @@ Vtotal = f[path +'Vtotal'][:]
 Vincp = np.abs(np.fft.fft(Vinc))**2
 freqs = np.fft.fftfreq(Vincp.size, d=dt)
 Vscatp = np.abs(np.fft.fft(Vscat))**2
-s11 = np.abs(Vscatp / Vincp)
+s11 = Vscatp / Vincp
 
 # Convert to decibels
 Vincp = 10 * np.log10(Vincp)
 Vscatp = 10 * np.log10(Vscatp)
 s11 = 10 * np.log10(s11)
 
-# Set plotting range to a frequency
+# Set plotting range to -60dB from maximum power
-pltrange = np.where(freqs > 2e9)[0][0]
+pltrange = np.where((np.amax(Vincp) - Vincp) > 60)[0][0] + 1
-pltrange = np.s_[1:pltrange]
+pltrange = np.s_[0:pltrange]
 
 # Plot incident voltage
-fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(nrows=3, ncols=2, num='Incident and scattered voltages', figsize=(20, 10), facecolor='w', edgecolor='w')
+plt.subplots(num='Transmission line voltages & s11 parameter', figsize=(20, 10), facecolor='w', edgecolor='w')
+gs = gridspec.GridSpec(3, 2)
+ax1 = plt.subplot(gs[0, 0])
 ax1.plot(time, Vinc, 'r', lw=2, label='Vinc')
 ax1.set_xlabel('Time [ns]')
-ax1.set_ylabel('Incident (field) voltage [V]')
+ax1.set_ylabel('Incident voltage [V]')
 ax1.set_xlim([0, np.amax(time)])
 ax1.grid()
 
 # Plot frequency spectra of incident voltage
+ax2 = plt.subplot(gs[0, 1])
 markerline, stemlines, baseline = ax2.stem(freqs[pltrange]/1e9, Vincp[pltrange], '-.')
+plt.setp(baseline, 'linewidth', 0)
 plt.setp(stemlines, 'color', 'r')
 plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
+ax2.plot(freqs[pltrange]/1e9, Vincp[pltrange], 'r', lw=2)
 ax2.set_xlabel('Frequency [GHz]')
-ax2.set_ylabel('Power [dB]')
+ax2.set_ylabel('Incident voltage spectra [dB]')
 ax2.grid()
 
-# Plot scattered voltage
+# Plot scattered (field) voltage
+ax3 = plt.subplot(gs[1, 0])
 ax3.plot(time, Vscat, 'r', lw=2, label='Vscat')
 ax3.set_xlabel('Time [ns]')
 ax3.set_ylabel('Scattered (field) voltage [V]')
@@ -83,20 +91,26 @@ ax3.set_xlim([0, np.amax(time)])
 ax3.grid()
 
 # Plot frequency spectra of scattered voltage
+ax4 = plt.subplot(gs[1, 1])
 markerline, stemlines, baseline = ax4.stem(freqs[pltrange]/1e9, Vscatp[pltrange], '-.')
+plt.setp(baseline, 'linewidth', 0)
 plt.setp(stemlines, 'color', 'r')
 plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
+ax4.plot(freqs[pltrange]/1e9, Vscatp[pltrange], 'r', lw=2)
 ax4.set_xlabel('Frequency [GHz]')
-ax4.set_ylabel('Power [dB]')
+ax4.set_ylabel('Scattered (field) voltage spectra [dB]')
 ax4.grid()
 
 # Plot frequency spectra of s11
-markerline, stemlines, baseline = ax6.stem(freqs[pltrange]/1e9, s11[pltrange], '-.')
+ax5 = plt.subplot(gs[2, 1])
+markerline, stemlines, baseline = ax5.stem(freqs[pltrange]/1e9, s11[pltrange], '-.')
+plt.setp(baseline, 'linewidth', 0)
 plt.setp(stemlines, 'color', 'r')
 plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
-ax6.set_xlabel('Frequency [GHz]')
+ax5.plot(freqs[pltrange]/1e9, s11[pltrange], 'r', lw=2)
-ax6.set_ylabel('Power [dB]')
+ax5.set_xlabel('Frequency [GHz]')
-ax6.grid()
+ax5.set_ylabel('s11 [dB]')
+ax5.grid()
 
 plt.show()
 f.close()