diff --git a/docs/source/user_libs_opt_taguchi.rst b/docs/source/user_libs_opt_taguchi.rst
index a1953933..922c3aa1 100644
--- a/docs/source/user_libs_opt_taguchi.rst
+++ b/docs/source/user_libs_opt_taguchi.rst
@@ -45,15 +45,15 @@ Package overview
 .. code-block:: none
 
     antenna_bowtie_opt.in
+    fitness_functions.py
     OA_9_4_3_2.npy
     OA_18_7_3_2.npy
-    optimisation_taguchi_fitness.py
-    optimisation_taguchi_plot.py
+    plot_results.py
 
 * ``antenna_bowtie_opt.in`` is a example model of a bowtie antenna where values of loading resistors are optimised.
+* ``fitness_functions.py`` is a module containing fitness functions. There are some pre-built ones but users should add their own here.
 * ``OA_9_4_3_2.npy`` and ``OA_18_7_3_2.npy`` are NumPy archives containing pre-built OAs from http://neilsloane.com/oadir/
-* ``optimisation_taguchi_fitness.py`` is a module containing fitness functions. There are some pre-built ones but users should add their own here.
-* ``optimisation_taguchi_plot.py`` is a module for plotting the results, such as parameter values and convergence history, from an optimisation process when it has completed.
+* ``plot_results.py`` is a module for plotting the results, such as parameter values and convergence history, from an optimisation process when it has completed.
 
 Implementation
 --------------
@@ -67,7 +67,7 @@ The process by which Taguchi's method optimises parameters is illustrated in the
 
 In stage 1a, one of the 2 pre-built OAs will automatically be chosen depending on the number of parameters to optimise. Currently, up to 7 independent parameters can be optimised, although a method to construct OAs of any size is under testing.
 
-In stage 1b, a fitness function is required to set a goal against which to compare results from the optimisation process. A number of pre-built fitness functions can be found in the ``optimisation_taguchi_fitness.py`` module, e.g. ``minvalue``, ``maxvalue`` and ``xcorr``. Users can also easily add their own fitness functions to this module. All fitness functions must take two arguments:
+In stage 1b, a fitness function is required to set a goal against which to compare results from the optimisation process. A number of pre-built fitness functions can be found in the ``fitness_functions.py`` module, e.g. ``minvalue``, ``maxvalue`` and ``xcorr``. Users can also easily add their own fitness functions to this module. All fitness functions must take two arguments:
 
 * ``filename`` a string containing the full path and filename of the output file
 * ``args`` a dictionary which can contain any number of additional arguments for the function, e.g. names (IDs) of outputs (rxs) from input file
@@ -114,11 +114,11 @@ The bowtie design features 3 vertical slots (y-direction) in each arm of the bow
     :language: none
     :linenos:
 
-The first part of the input file (lines 1-7) contains the parameters to optimise, their initial ranges, and fitness function information for the optimisation process. Three parameters representing the resistor values are defined with ranges between 0.1 :math:`\Omega` and 5 :math:`k\Omega`. A fitness function called ``maxvalue`` with a stopping criterion of 50V/m. The output point in the model that will be used in the optimisation is specified as having the name ``Ex60mm``. Finally, a limit of 5 iterations is placed on the optimisation process, i.e. it will stop after 5 iterations irrespectively of whether it has reached the target of 50V/m.
+The first part of the input file (lines 1-7) contains the parameters to optimise, their initial ranges, and fitness function information for the optimisation process. Three parameters representing the resistor values are defined with ranges between 0.1 :math:`\Omega` and 5 :math:`k\Omega`. A pre-built fitness function called ``maxvalue`` is specified with a stopping criterion of 50V/m. The output point in the model that will be used in the optimisation is specified as having the name ``Ex60mm``. Finally, a limit of 5 iterations is placed on the optimisation process, i.e. it will stop after 5 iterations irrespectively of whether it has reached the target of 50V/m.
 
 The next part of the input file (lines 9-93) contains the model. For the most part there is nothing special about the way the model is defined - a mixture of Python, NumPy and functional forms of the input commands (available by importing the module ``input_cmd_funcs``) are used. However, it is worth pointing out how the values of the parameters to optimise are accessed. On line 29 a NumPy array of the values of the resistors is created. The values are accessed using their names as keys to the ``optparams`` dictionary. On line 30 the values of the resistors are converted to conductivities, which are used to create new materials (line 34-35). The resistors are then built by applying the materials to cell edges (e.g. lines 55-62). The output point in the model in specifed with the name ``Ex60mm`` and as having only an ``Ex`` field output (line 42).
 
-The optimisation process is run on the model using the ``--opt-taguchi`` command line flag:
+The optimisation process is run on the model using the ``--opt-taguchi`` command line flag.
 
 .. code-block:: none
 
@@ -127,3 +127,8 @@ The optimisation process is run on the model using the ``--opt-taguchi`` command
 Results
 ^^^^^^^
 
+When the optimisation has completed a summary will be printed showing histories of the parameter values and the fitness metric. These values are also saved (pickled) to file and can be plotted using the ``plot_results.py`` module, for example:
+
+.. code-block:: none
+
+    python -m user_libs.optimisation_taguchi.plot_results antenna_bowtie_opt_hist.pickle
diff --git a/gprMax/optimisation_taguchi.py b/gprMax/optimisation_taguchi.py
index 64983ac6..cf5310f8 100644
--- a/gprMax/optimisation_taguchi.py
+++ b/gprMax/optimisation_taguchi.py
@@ -62,7 +62,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace):
     optparamshist = OrderedDict((key, list()) for key in optparams)
     
     # Import specified fitness function
-    fitness_metric = getattr(importlib.import_module('user_libs.optimisation_taguchi.optimisation_taguchi_fitness'), fitness['name'])
+    fitness_metric = getattr(importlib.import_module('user_libs.optimisation_taguchi.fitness_functions'), fitness['name'])
 
     # Select OA
     OA, N, cols, k, s, t = construct_OA(optparams)
diff --git a/tools/plot_antenna_params.py b/tools/plot_antenna_params.py
index beb9c2c7..d8b15772 100644
--- a/tools/plot_antenna_params.py
+++ b/tools/plot_antenna_params.py
@@ -16,15 +16,13 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-import os, argparse
+import argparse, os
 import h5py
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.gridspec as gridspec
 #import scipy.io as sio
 
-moduledirectory = os.path.dirname(os.path.abspath(__file__))
-
 """Plots antenna parameters (s11 parameter and input impedance and admittance) from an output file containing a transmission line source."""
 
 # Parse command line arguments
@@ -60,6 +58,7 @@ Vinc = f[path + 'Vinc'][:]
 Iinc = f[path + 'Iinc'][:]
 Vtotal = f[path +'Vtotal'][:]
 Itotal = f[path +'Itotal'][:]
+Vrec = f['/rxs/rx1/Ex'][:] * -1
 f.close()
 Vref = Vtotal - Vinc
 Iref = Itotal - Iinc
@@ -72,13 +71,14 @@ delaycorrection = np.exp(-1j * 2 * np.pi * freqs * (dt / 2))
 
 # Calculate s11
 s11 = np.abs(np.fft.fft(Vref) * delaycorrection) / np.abs(np.fft.fft(Vinc) * delaycorrection)
+s21 = np.abs(np.fft.fft(Vrec)) / np.abs(np.fft.fft(Vinc) * delaycorrection)
 
 # Calculate input impedance
 zin = (np.fft.fft(Vtotal) * delaycorrection) / np.fft.fft(Itotal)
 
 # Load MoM zin from MATLAB antenna toolbox
 #MoM = {}
-#sio.loadmat(moduledirectory + '/../tests/numerical/vs_MoM_MATLAB/antenna_bowtie_fs/antenna_bowtie_fs_MoM.mat', MoM)
+#sio.loadmat('/../tests/numerical/vs_MoM_MATLAB/antenna_bowtie_fs/antenna_bowtie_fs_MoM.mat', MoM)
 
 # Calculate input admittance
 yin = np.fft.fft(Itotal) / (np.fft.fft(Vtotal) * delaycorrection)
@@ -91,6 +91,7 @@ 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)
+s21 = 20 * np.log10(s21)
 
 # Set plotting range
 pltrangemin = 1
@@ -240,7 +241,7 @@ ax.grid()
 # Figure 2
 # Plot frequency spectra of s11
 fig2, ax = plt.subplots(num='Antenna parameters', figsize=(20, 12), facecolor='w', edgecolor='w')
-gs2 = gridspec.GridSpec(3, 2, hspace=0.5)
+gs2 = gridspec.GridSpec(2, 2, hspace=0.5)
 ax = plt.subplot(gs2[0, 0])
 markerline, stemlines, baseline = ax.stem(freqs[pltrange], s11[pltrange], '-.')
 plt.setp(baseline, 'linewidth', 0)
@@ -251,7 +252,21 @@ ax.set_title('s11')
 ax.set_xlabel('Frequency [Hz]')
 ax.set_ylabel('Power [dB]')
 #ax.set_xlim([0.88e9, 1.02e9])
-#ax.set_ylim([-20, 0])
+ax.set_ylim([-20, 0])
+ax.grid()
+
+# Plot frequency spectra of s21
+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()
 
 # Plot input resistance (real part of impedance)
@@ -266,7 +281,7 @@ 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, 350])
+ax.set_ylim([0, 300])
 ax.grid()
 
 # Plot input reactance (imaginery part of impedance)
@@ -280,36 +295,36 @@ 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([-1400, 200])
+ax.set_ylim([-200, 100])
 ax.grid()
 
-# 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()
-
-# 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()
+## 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()
+#
+## 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()
 
 # Figure 3 - Comparison of numerical modelling techniques
 #fig3, ax = plt.subplots(num='FDTD vs MoM', figsize=(20, 5), facecolor='w', edgecolor='w')
diff --git a/user_libs/optimisation_taguchi/antenna_bowtie_opt.in b/user_libs/optimisation_taguchi/antenna_bowtie_opt.in
index 53f1bfdb..1138a2ce 100644
--- a/user_libs/optimisation_taguchi/antenna_bowtie_opt.in
+++ b/user_libs/optimisation_taguchi/antenna_bowtie_opt.in
@@ -2,7 +2,7 @@
 optparams['resinner'] = [0.1, 5000]
 optparams['resmiddle'] = [0.1, 5000]
 optparams['resouter'] = [0.1, 5000]
-fitness = {'name': 'compactness', 'stop': 30, 'args': {'outputs': 'Ex60mm'}}
+fitness = {'name': 'maxvalue', 'stop': 50, 'args': {'outputs': 'Ex60mm'}}
 maxiterations = 5
 #end_taguchi:
 
diff --git a/user_libs/optimisation_taguchi/optimisation_taguchi_fitness.py b/user_libs/optimisation_taguchi/fitness_functions.py
similarity index 100%
rename from user_libs/optimisation_taguchi/optimisation_taguchi_fitness.py
rename to user_libs/optimisation_taguchi/fitness_functions.py
diff --git a/user_libs/optimisation_taguchi/optimisation_taguchi_plot.py b/user_libs/optimisation_taguchi/plot_results.py
similarity index 94%
rename from user_libs/optimisation_taguchi/optimisation_taguchi_plot.py
rename to user_libs/optimisation_taguchi/plot_results.py
index d8dfa422..5b4829be 100644
--- a/user_libs/optimisation_taguchi/optimisation_taguchi_plot.py
+++ b/user_libs/optimisation_taguchi/plot_results.py
@@ -12,7 +12,7 @@ from gprMax.optimisation_taguchi import plot_optimisation_history
 """Plots the results (pickled to file) from a Taguchi optimisation process."""
 
 # Parse command line arguments
-parser = argparse.ArgumentParser(description='Plots the results (pickled to file) from a Taguchi optimisation process.', usage='cd gprMax; python -m user_libs.optimisation_taguchi_plot picklefile')
+parser = argparse.ArgumentParser(description='Plots the results (pickled to file) from a Taguchi optimisation process.', usage='cd gprMax; python -m user_libs.optimisation_taguchi.plot_results picklefile')
 parser.add_argument('picklefile', help='name of file including path')
 args = parser.parse_args()