diff --git a/docs/source/images/AustinMan.png b/docs/source/images/user_libs/AustinMan.png
similarity index 100%
rename from docs/source/images/AustinMan.png
rename to docs/source/images/user_libs/AustinMan.png
diff --git a/docs/source/images/AustinMan_head.png b/docs/source/images/user_libs/AustinMan_head.png
similarity index 100%
rename from docs/source/images/AustinMan_head.png
rename to docs/source/images/user_libs/AustinMan_head.png
diff --git a/docs/source/images/user_libs/antenna_bowtie_opt.png b/docs/source/images/user_libs/antenna_bowtie_opt.png
new file mode 100644
index 00000000..36b8467a
Binary files /dev/null and b/docs/source/images/user_libs/antenna_bowtie_opt.png differ
diff --git a/docs/source/images/taguchi_process.png b/docs/source/images/user_libs/taguchi_process.png
similarity index 100%
rename from docs/source/images/taguchi_process.png
rename to docs/source/images/user_libs/taguchi_process.png
diff --git a/docs/source/user_libs_austinman.rst b/docs/source/user_libs_austinman.rst
index d60b4d38..5d14ccad 100644
--- a/docs/source/user_libs_austinman.rst
+++ b/docs/source/user_libs_austinman.rst
@@ -22,7 +22,7 @@ Information
 
 AustinMan and AustinWoman (http://bit.ly/AustinMan) are open source electromagnetic voxel models of the human body, which are developed by the Computational Electromagnetics Group (http://www.ece.utexas.edu/research/areas/electromagnetics-acoustics) at The University of Texas at Austin (http://www.utexas.edu). The models are based on data from the National Library of Medicine’s Visible Human Project (https://www.nlm.nih.gov/research/visible/visible_human.html).
 
-.. figure:: images/AustinMan_head.png
+.. figure:: images/user_libs/AustinMan_head.png
     :width: 600 px
 
     FDTD geometry mesh showing the head of the AustinMan model (2x2x2mm^3).
@@ -68,7 +68,7 @@ To insert a 2x2x2mm^3 AustinMan with the lower left corner 40mm from the origin
 
 For further information on the `#geometry_objects_file` see the section on object contruction commands in the :ref:`Input commands section <commands>`.
 
-.. figure:: images/AustinMan.png
+.. figure:: images/user_libs/AustinMan.png
     :width: 300 px
 
     FDTD geometry mesh showing the AustinMan body model (2x2x2mm^3).
diff --git a/docs/source/user_libs_opt_taguchi.rst b/docs/source/user_libs_opt_taguchi.rst
index 972cdccb..d16ed7a7 100644
--- a/docs/source/user_libs_opt_taguchi.rst
+++ b/docs/source/user_libs_opt_taguchi.rst
@@ -22,7 +22,7 @@ Information
     #
     # Please use the attribution at http://dx.doi.org/10.1190/1.3548506
 
-The package features an optimisation technique based on Taguchi's method. It allows the user to define parameters in an input file and optimise their values based on a fitness function.
+The package features an optimisation technique based on Taguchi's method. It allows users to define parameters in an input file and optimise their values based on a fitness function, for example it can be used to optimise material properties or geometry in a simulation.
 
 
 Taguchi's method
@@ -44,11 +44,13 @@ Package overview
 
 .. code-block:: none
 
+    antenna_bowtie_opt.in
     OA_9_4_3_2.npy
     OA_18_7_3_2.npy
     optimisation_taguchi_fitness.py
     optimisation_taguchi_plot.py
 
+* ``antenna_bowtie_opt.in`` is a example model of a bowtie antenna where values of loading resistors are optimised.
 * ``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.
@@ -58,7 +60,7 @@ Implementation
 
 The process by which Taguchi's method optimises parameters is illustrated in the following figure.
 
-.. figure:: images/taguchi_process.png
+.. figure:: images/user_libs/taguchi_process.png
     :width: 300 px
 
     Process associated with Taguchi's method.
@@ -93,3 +95,8 @@ Example
 
 The following example demonstrates using the Taguchi optimisation process to optimise values of loading resistors used in a bowtie antenna. The bowtie design features 3 slots in each arm of the bowtie where loading resistors are placed, and a substrate with a perimittivity of 4.8 is used. The antenna is modelled in free space, and an output point (the electric field value) is specified at a distance of 60 mm from the feed of the bowtie.
 
+.. figure:: images/user_libs/antenna_bowtie_opt.png
+    :width: 600 px
+
+    FDTD geometry mesh showing bowtie antenna with slots and loading resistors.
+
diff --git a/user_libs/optimisation_taguchi/antenna_bowtie_opt.in b/user_libs/optimisation_taguchi/antenna_bowtie_opt.in
new file mode 100644
index 00000000..53f1bfdb
--- /dev/null
+++ b/user_libs/optimisation_taguchi/antenna_bowtie_opt.in
@@ -0,0 +1,93 @@
+#taguchi:
+optparams['resinner'] = [0.1, 5000]
+optparams['resmiddle'] = [0.1, 5000]
+optparams['resouter'] = [0.1, 5000]
+fitness = {'name': 'compactness', 'stop': 30, 'args': {'outputs': 'Ex60mm'}}
+maxiterations = 5
+#end_taguchi:
+
+#python:
+
+import numpy as np
+from gprMax.input_cmd_funcs import *
+
+title = 'antenna_bowtie_opt'
+print('#title: {}'.format(title))
+
+domain = domain(0.180, 0.120, 0.160)
+dxdydz = dx_dy_dz(0.001, 0.001, 0.001)
+timewindow = time_window(5e-9)
+fr4_dims = (0.120, 0.060, 0.002)
+bowtie_dims = (0.050, 0.040) # Length, height
+flare_angle = np.arctan((bowtie_dims[1]/2) / bowtie_dims[0])
+tx_pos = (domain[0]/2, domain[1]/2, 0.050)
+
+# Vertical slot positions, relative to feed position, i.e. txpos[0]
+vcut_pos = (0.014, 0.027, 0.038)
+
+# Loading resistor values
+res = np.array([optparams['resinner'], optparams['resmiddle'], optparams['resouter']])
+rescond = ((1 / res) * (dxdydz[1] / (dxdydz[0] * dxdydz[2]))) / 2 # Divide by number of parallel edges per resistor
+
+# Materials
+material(4.8, 0, 1, 0, 'fr4')
+for i in range(len(res)):
+    material(1, rescond[i], 1, 0, 'res' + str(i + 1))
+
+# Source excitation and type
+print('#waveform: gaussian 1 2e9 mypulse')
+print('#transmission_line: x {:g} {:g} {:g} 50 mypulse'.format(tx_pos[0], tx_pos[1], tx_pos[2]))
+
+# Output point - distance from tx_pos in z direction
+print('#rx: {:g} {:g} {:g} Ex60mm Ex'.format(tx_pos[0], tx_pos[1], tx_pos[2] + 0.060))
+
+# Bowtie - upper x half
+triangle(tx_pos[0], tx_pos[1], tx_pos[2], tx_pos[0] + bowtie_dims[0], tx_pos[1] - bowtie_dims[1]/2, tx_pos[2], tx_pos[0] + bowtie_dims[0], tx_pos[1] + bowtie_dims[1]/2, tx_pos[2], 0, 'pec')
+
+# Bowtie - upper x half - vertical cuts
+for i in range(len(vcut_pos)):
+	for j in range(int(bowtie_dims[1] / dxdydz[2])):
+		edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - bowtie_dims[1]/2 + j * dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - bowtie_dims[1]/2 + j * dxdydz[1], tx_pos[2], 'free_space')
+
+# Bowtie - upper x half - vertical cuts - loading
+for i in range(len(vcut_pos)):
+    gap = ((vcut_pos[i] * np.tan(flare_angle) * 2) - 4*dxdydz[1]) / 5
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (0.5 * gap), tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (0.5 * gap), tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+
+# Bowtie - lower x half
+triangle(tx_pos[0] + dxdydz[0], tx_pos[1], tx_pos[2], tx_pos[0] - bowtie_dims[0], tx_pos[1] - bowtie_dims[1]/2, tx_pos[2], tx_pos[0] - bowtie_dims[0], tx_pos[1] + bowtie_dims[1]/2, tx_pos[2], 0, 'pec')
+
+# Bowtie - lower x half - cuts for loading
+for i in range(len(vcut_pos)):
+	for j in range(int(bowtie_dims[1] / dxdydz[2])):
+		edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - bowtie_dims[1]/2  + j * dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - bowtie_dims[1]/2  + j * dxdydz[1], tx_pos[2], 'free_space')
+
+# Bowtie - lower x half - vertical cuts - loading
+for i in range(len(vcut_pos)):
+    gap = ((vcut_pos[i] * np.tan(flare_angle) * 2) - 4*dxdydz[1]) / 5
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (0.5 * gap), tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (0.5 * gap), tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+    edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
+
+# PCB
+box(tx_pos[0] - fr4_dims[0]/2, tx_pos[1] - fr4_dims[1]/2, tx_pos[2] - fr4_dims[2], tx_pos[0] + fr4_dims[0]/2, tx_pos[1] + fr4_dims[1]/2, tx_pos[2], 'fr4')
+
+# Detailed geometry view of PCB and bowtie
+#geometry_view(tx_pos[0] - fr4_dims[0]/2, tx_pos[1] - fr4_dims[1]/2, tx_pos[2] - fr4_dims[2], tx_pos[0] + fr4_dims[0]/2, tx_pos[1] + fr4_dims[1]/2, tx_pos[2] + dxdydz[2], dxdydz[0], dxdydz[1], dxdydz[2], title + '_tx', type='f')
+
+# Geometry view of entire domain
+#geometry_view(0, 0, 0, domain[0], domain[1], domain[2], dxdydz[0], dxdydz[1], dxdydz[2], title)
+
+#end_python:
\ No newline at end of file