Changed 'field' terminology to 'output'.

tools/plot_Ascan.py

@@ -20,19 +20,20 @@ import os, argparse
 import h5py
 import numpy as np
 import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
 
 from gprMax.exceptions import CmdInputError
 
-"""Plots electric and magnetic fields from any receiver points in the given output file. Each receiver point is plotted in a new figure window."""
+"""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."""
 
-# Fields that can be plotted
+# Outputs that can be plotted
-fieldslist = ['Ex', 'Hx', 'Ey', 'Hy', 'Ez', 'Hz']
+outputslist = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz']
 
 # Parse command line arguments
-parser = argparse.ArgumentParser(description='Plots electric and magnetic fields 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 = 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('--fields', help='list of fields to be plotted, i.e. Ex Ey Ez', default=fieldslist, nargs='+')
+parser.add_argument('--outputs', help='list of outputs to be plotted, i.e. Ex Ey Ez', default=outputslist, nargs='+')
-parser.add_argument('-fft', action='store_true', default=False, help='plot FFT (single field component must be specified)')
+parser.add_argument('-fft', action='store_true', default=False, help='plot FFT (single output must be specified)')
 args = parser.parse_args()
 
 # Open output file and read some attributes
@@ -44,28 +45,28 @@ iterations = f.attrs['Iterations']
 time = np.arange(0, dt * iterations, dt)
 time = time / 1e-9
 
-# Check for valid field names
+# Check for valid output names
-for field in args.fields:
+for output in args.outputs:
-    if field not in fieldslist:
+    if output not in outputslist:
-        raise CmdInputError('{} not allowed. Options are: Ex Ey Ez Hx Hy Hz'.format(field))
+        raise CmdInputError('{} not allowed. Options are: Ex Ey Ez Hx Hy Hz'.format(output))
 
-# Check for single field component when doing a FFT
+# Check for single output component when doing a FFT
 if args.fft:
-    if not len(args.fields) == 1:
+    if not len(args.outputs) == 1:
-        raise CmdInputError('A single field component must be specified when using the -fft option')
+        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) + '/'
     
-    # If only a single field is required, create one subplot
+    # If only a single output is required, create one subplot
-    if len(args.fields) == 1:
+    if len(args.outputs) == 1:
-        fielddata = f[path + args.fields[0]][:]
+        outputdata = f[path + args.outputs[0]][:]
         
         # Plotting if FFT required
         if args.fft:
             # Calculate magnitude of frequency spectra of waveform
-            power = 10 * np.log10(np.abs(np.fft.fft(fielddata))**2)
+            power = 10 * np.log10(np.abs(np.fft.fft(outputdata))**2)
             freqs = np.fft.fftfreq(power.size, d=dt)
 
             # Shift powers so that frequency with maximum power is at zero decibels
@@ -75,11 +76,11 @@ for rx in range(1, nrx + 1):
             pltrange = np.where((np.amax(power) - power) > 60)[0][0] + 1
             pltrange = np.s_[0:pltrange]
 
-            # Plot time history of field component
+            # 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, fielddata, 'r', lw=2, label=args.fields[0])
+            line1 = ax1.plot(time, outputdata, 'r', lw=2, label=args.outputs[0])
             ax1.set_xlabel('Time [ns]')
-            ax1.set_ylabel(args.fields[0] + ' field strength [V/m]')
+            ax1.set_ylabel(args.outputs[0] + ' field strength [V/m]')
             ax1.set_xlim([0, np.amax(time)])
             ax1.grid()
 
@@ -94,10 +95,10 @@ for rx in range(1, nrx + 1):
             ax2.grid()
             
             # Change colours and labels for magnetic field components
-            if 'H' in args.fields[0]:
+            if 'H' in args.outputs[0]:
                 plt.setp(line1, color='b')
                 plt.setp(line2, color='b')
-                plt.setp(ax1, ylabel=args.fields[0] + ' field strength [A/m]')
+                plt.setp(ax1, ylabel=args.outputs[0] + ' field strength [A/m]')
                 plt.setp(stemlines, 'color', 'b')
                 plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
             
@@ -105,38 +106,61 @@ for rx in range(1, nrx + 1):
     
         # Plotting if no FFT required
         else:
-            fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [ns]', ylabel=args.fields[0] + ' field strength [V/m]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
+            fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [ns]', ylabel=args.outputs[0] + ' field strength [V/m]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
-            line = ax.plot(time, fielddata,'r', lw=2, label=args.fields[0])
+            line = ax.plot(time, outputdata,'r', lw=2, label=args.outputs[0])
             ax.set_xlim([0, np.amax(time)])
             ax.grid()
             
-            if 'H' in args.fields[0]:
+            if 'H' in args.outputs[0]:
                 plt.setp(line, color='b')
-                plt.setp(ax, ylabel=args.fields[0] + ' field strength [A/m]')
+                plt.setp(ax, ylabel=args.outputs[0] + ', field strength [A/m]')
+            elif 'I' in args.outputs[0]:
+                plt.setp(line, color='b')
+                plt.setp(ax, ylabel=args.outputs[0] + ', current [A]')
 
-    # If multiple fields required, creat all six subplots and populate only the specified ones
+    # If multiple fields required, creat all nine subplots and populate only the specified ones
     else:
-        fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
+        fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [ns]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
-        for field in args.fields:
+        gs = gridspec.GridSpec(3, 3, hspace=0.3)
-            fielddata = f[path + field][:]
+#        fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
-            if field == 'Ex':
+        for output in args.outputs:
-                ax1.plot(time, fielddata,'r', lw=2, label=field)
+            outputdata = f[path + output][:]
-                ax1.set_ylabel(field + ', field strength [V/m]')
+            if output == 'Ex':
-            elif field == 'Ey':
+                ax = plt.subplot(gs[0, 0])
-                ax3.plot(time, fielddata,'r', lw=2, label=field)
+                ax.plot(time, outputdata,'r', lw=2, label=output)
-                ax3.set_ylabel(field + ', field strength [V/m]')
+                ax.set_ylabel(output + ', field strength [V/m]')
-            elif field == 'Ez':
+            elif output == 'Ey':
-                ax5.plot(time, fielddata,'r', lw=2, label=field)
+                ax = plt.subplot(gs[1, 0])
-                ax5.set_ylabel(field + ', field strength [V/m]')
+                ax.plot(time, outputdata,'r', lw=2, label=output)
-            elif field == 'Hx':
+                ax.set_ylabel(output + ', field strength [V/m]')
-                ax2.plot(time, fielddata,'b', lw=2, label=field)
+            elif output == 'Ez':
-                ax2.set_ylabel(field + ', field strength [A/m]')
+                ax = plt.subplot(gs[2, 0])
-            elif field == 'Hy':
+                ax.plot(time, outputdata,'r', lw=2, label=output)
-                ax4.plot(time, fielddata,'b', lw=2, label=field)
+                ax.set_ylabel(output + ', field strength [V/m]')
-                ax4.set_ylabel(field + ', field strength [A/m]')
+            elif output == 'Hx':
-            elif field == 'Hz':
+                ax = plt.subplot(gs[0, 1])
-                ax6.plot(time, fielddata,'b', lw=2, label=field)
+                ax.plot(time, outputdata,'b', lw=2, label=output)
-                ax6.set_ylabel(field + ', field strength [A/m]')
+                ax.set_ylabel(output + ', field strength [A/m]')
+            elif output == 'Hy':
+                ax = plt.subplot(gs[1, 1])
+                ax.plot(time, outputdata,'b', lw=2, label=output)
+                ax.set_ylabel(output + ', field strength [A/m]')
+            elif output == 'Hz':
+                ax = plt.subplot(gs[2, 1])
+                ax.plot(time, outputdata,'b', lw=2, label=output)
+                ax.set_ylabel(output + ', field strength [A/m]')
+            elif output == 'Ix':
+                ax = plt.subplot(gs[0, 2])
+                ax.plot(time, outputdata,'b', lw=2, label=output)
+                ax.set_ylabel(output + ', current [A]')
+            elif output == 'Iy':
+                ax = plt.subplot(gs[1, 2])
+                ax.plot(time, outputdata,'b', lw=2, label=output)
+                ax.set_ylabel(output + ', current [A]')
+            elif output == 'Iz':
+                ax = plt.subplot(gs[2, 2])
+                ax.plot(time, outputdata,'b', lw=2, label=output)
+                ax.set_ylabel(output + ', current [A]')
         for ax in fig.axes:
             ax.set_xlim([0, np.amax(time)])
             ax.grid()

tools/plot_Bscan.py

@@ -25,39 +25,41 @@ from gprMax.exceptions import CmdInputError
 
 """Plots a B-scan image."""
 
-# Fields that can be plotted
+# Outputs that can be plotted
-fieldslist = ['Ex', 'Hx', 'Ey', 'Hy', 'Ez', 'Hz']
+outputslist = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz']
 
 # Parse command line arguments
 parser = argparse.ArgumentParser(description='Plots B-scan.', usage='cd gprMax; python -m tools.plot_Bscan outputfile --field fieldcomponent')
 parser.add_argument('outputfile', help='name of output file including path')
-parser.add_argument('--field', help='name of field to be plotted, i.e. Ex Ey Ez')
+parser.add_argument('--output', help='name of output to be plotted, i.e. Ex Ey Ez')
 args = parser.parse_args()
 
-# Check for valid field name
+# Check for valid output name
-if args.field not in fieldslist:
+if args.output not in outputslist:
-    raise CmdInputError('{} not allowed. Options are: Ex Ey Ez Hx Hy Hz'.format(args.field))
+    raise CmdInputError('{} not allowed. Options are: Ex Ey Ez Hx Hy Hz Ix Iy Iz'.format(args.output))
 
 # Open output file and read some attributes
 f = h5py.File(args.outputfile, 'r')
 path = '/rxs/rx1'
-data = f[path + '/' + args.field]
+outputdata = f[path + '/' + args.output]
 
 # Check that there is more than one A-scan present
-if data.shape[1] == 1:
+if outputdata.shape[1] == 1:
     raise CmdInputError('{} contains only a single A-scan.'.format(args.outputfile))
 
 # Plot B-scan image
 fig = plt.figure(num=args.outputfile, figsize=(20, 10), facecolor='w', edgecolor='w')
-plt.imshow(data, extent=[0, data.shape[1], data.shape[0]*f.attrs['dt']*1e9, 0], interpolation='nearest', aspect='auto', cmap='seismic', vmin=-np.amax(np.abs(data)), vmax=np.amax(np.abs(data)))
+plt.imshow(outputdata, extent=[0, outputdata.shape[1], outputdata.shape[0]*f.attrs['dt']*1e9, 0], interpolation='nearest', aspect='auto', cmap='seismic', vmin=-np.amax(np.abs(outputdata)), vmax=np.amax(np.abs(outputdata)))
 plt.xlabel('Trace number')
 plt.ylabel('Time [ns]')
 plt.grid()
 cb = plt.colorbar()
-if 'E' in args.field:
+if 'E' in args.output:
     cb.set_label('Field strength [V/m]')
-elif 'H' in args.field:
+elif 'H' in args.output:
     cb.set_label('Field strength [A/m]')
+elif 'I' in args.output:
+    cb.set_label('Current [A]')
 
 plt.show()
 #fig.savefig(os.path.splitext(os.path.abspath(file))[0] + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)