diff --git a/.gitignore b/.gitignore
index db06575cf9f291d1fb9fa253fb5146c065523dd9..da42a86becea35fec627a26e6ee74710b0e2a699 100644
--- a/.gitignore
+++ b/.gitignore
@@ -80,9 +80,9 @@ tests/testRiemannHLLC
 tests/testMatrixInversion
 
 theory/latex/swift.pdf
-theory/kernel/kernels.pdf
-theory/kernel/kernel_derivatives.pdf
-theory/kernel/kernel_definitions.pdf
+theory/SPH/kernel/kernels.pdf
+theory/SPH/kernel/kernel_derivatives.pdf
+theory/SPH/kernel/kernel_definitions.pdf
 theory/paper_pasc/pasc_paper.pdf
 
 m4/libtool.m4
diff --git a/examples/SedovBlast_3D/plotSolution.py b/examples/SedovBlast_3D/plotSolution.py
index 1eea372b08e084a37c001f9c53a61667277c765b..8eca94c98acd7e081cdfc2c785aa1d19d601e81e 100644
--- a/examples/SedovBlast_3D/plotSolution.py
+++ b/examples/SedovBlast_3D/plotSolution.py
@@ -18,7 +18,7 @@
  # 
  ##############################################################################
 
-# Computes the analytical solution of the 2D Sedov blast wave.
+# Computes the analytical solution of the 3D Sedov blast wave.
 # The script works for a given initial box and dumped energy and computes the solution at a later time t.
 
 # Parameters
diff --git a/theory/SPH/kernel/kernel_definitions.tex b/theory/SPH/kernel/kernel_definitions.tex
index 8999636109ffadcbf148ce3c1fbccdc44feafe65..65ef5c82ad928398a69b4528c7914829f146ad6e 100644
--- a/theory/SPH/kernel/kernel_definitions.tex
+++ b/theory/SPH/kernel/kernel_definitions.tex
@@ -218,7 +218,7 @@ All kernels available in \swift are shown on Fig.~\ref{fig:sph:kernels}.
 
 \begin{figure}
 \includegraphics[width=\columnwidth]{kernel_derivatives.pdf}
-\caption{The first and secon derivatives of the kernel functions
+\caption{The first and second derivatives of the kernel functions
   available in \swift for a mean inter-particle separation $\langle
   x\rangle=1.5$ and a resolution $\eta=1.2348$.  A Gaussian kernel
   with the same smoothing length is shown for comparison.}
@@ -229,14 +229,13 @@ All kernels available in \swift are shown on Fig.~\ref{fig:sph:kernels}.
 \section{Kernel Derivatives}
 
 The derivatives of the kernel function have relatively simple
-expressions and are shown on Fig.~\ref{fig:sph:kernel_derivatives}.
+expressions and are shown on Fig.~\ref{fig:sph:kernel_derivatives}:
 
 \begin{eqnarray*}
- \vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|} \\
+ \vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|}, \\
  \frac{\partial W(\vec{x},h)}{\partial h} &=&- \frac{1}{h^4}\left[3f\left(\frac{|\vec{x}|}{h}\right) + 
-\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right]
+\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right].
 \end{eqnarray*}
-
 Note that for all the kernels listed above, $f'(0) = 0$. 
 
 \end{document}
diff --git a/theory/SPH/kernel/kernels.py b/theory/SPH/kernel/kernels.py
index 184379e5eafbcd12a1a47560ee88e02066da3942..9f0b853c6e9fce7180361457373f745cc2b8e804 100644
--- a/theory/SPH/kernel/kernels.py
+++ b/theory/SPH/kernel/kernels.py
@@ -1,5 +1,21 @@
-#!/usr/bin/env python2                                  
-# -*- coding: utf-8 -*-            
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016  Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program 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 Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
 import matplotlib
 matplotlib.use("Agg")
 from pylab import *
@@ -22,9 +38,9 @@ params = {
     'text.usetex': True,
 'figure.figsize' : (4.15,4.15),
 'figure.subplot.left'    : 0.14,
-'figure.subplot.right'   : 0.99  ,
-'figure.subplot.bottom'  : 0.08  ,
-'figure.subplot.top'     : 0.99  ,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.06,
+'figure.subplot.top'     : 0.99,
 'figure.subplot.wspace'  : 0.  ,
 'figure.subplot.hspace'  : 0.  ,
 'lines.markersize' : 6,
@@ -139,6 +155,10 @@ arrow(H_WendlandC6, 0.12*maxY , 0., -0.12*maxY*0.9, fc='y', ec='y', length_inclu
 plot([h, h], [0., maxY], 'k:', linewidth=0.5)
 text(h, maxY*0.35, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom')
 
+# Show sigma
+plot([h/2, h/2], [0., maxY], 'k:', linewidth=0.5)
+text(h/2, maxY*0.05, "$\\sigma\\equiv h/2$", rotation=90, backgroundcolor='w', ha='center', va='bottom')
+
 # Show <x>
 plot([dx, dx], [0., maxY], 'k:', linewidth=0.5)
 text(dx, maxY*0.35, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom')
@@ -163,6 +183,9 @@ plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$")
 # Show h
 plot([h, h], [0., 1.], 'k:', linewidth=0.5)
 
+# Show sigma
+plot([h/2, h/2], [0., 1.], 'k:', linewidth=0.5)
+
 # Show <x>
 plot([dx, dx], [0., 1.], 'k:', linewidth=0.5)
 
@@ -272,13 +295,16 @@ maxY = d_Gaussian(h/2, h)
 # Show h
 plot([h, h], [2*maxY, 0.1], 'k:', linewidth=0.5)
 
+# Show sigma
+plot([h/2, h/2], [2*maxY, 0.1], 'k:', linewidth=0.5)
+
 # Show <x>
 plot([dx, dx], [2*maxY, 0.1], 'k:', linewidth=0.5)
 
 
 xlim(0., 2.5*h)
 gca().xaxis.set_ticklabels([])
-ylim(1.2*maxY, -0.1*maxY)
+ylim(1.1*maxY, -0.1*maxY)
 xlabel("$r$", labelpad=0)
 ylabel("$\\partial W(r,h)/\\partial r$", labelpad=0.5)
 legend(loc="lower right")
@@ -288,12 +314,14 @@ legend(loc="lower right")
 subplot(212)
 
 maxY = d2_Gaussian(h,h)
+plot([h/2, h/2], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
 plot([h, h], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
-text(h, -3.*maxY, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom')
-
 plot([dx, dx], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
+text(h/2, -3.*maxY, "$\\sigma\\equiv h/2$", rotation=90, backgroundcolor='w', ha='center', va='bottom')
+text(h, -3.*maxY, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom')
 text(dx, -3.*maxY, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom')
 
+
 plot([0, 2.5*h], [0., 0.], 'k--', linewidth=0.7)
 plot(xx, d2_Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$")
 plot(xx, d2W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$")
@@ -304,7 +332,7 @@ plot(xx, d2W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$")
 plot(xx, d2W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$")
 
 xlim(0., 2.5*h)
-ylim(-3.2*maxY, 1.4*maxY)
+ylim(-3.2*maxY, 1.3*maxY)
 xlabel("$r$", labelpad=0)
 ylabel("$\\partial^2 W(r,h)/\\partial r^2$", labelpad=0.5)
 
diff --git a/theory/SPH/kernel/run.sh b/theory/SPH/kernel/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..61b469025b705f298b367963700e99fb13f13bd4
--- /dev/null
+++ b/theory/SPH/kernel/run.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+python kernels.py
+pdflatex kernel_definitions.tex
+pdflatex kernel_definitions.tex