diff --git a/.gitignore b/.gitignore
index da42a86becea35fec627a26e6ee74710b0e2a699..d8b6425dcb10a99f25201efffa2e1c10f0bca079 100644
--- a/.gitignore
+++ b/.gitignore
@@ -83,6 +83,7 @@ theory/latex/swift.pdf
theory/SPH/kernel/kernels.pdf
theory/SPH/kernel/kernel_derivatives.pdf
theory/SPH/kernel/kernel_definitions.pdf
+theory/SPH/flavours/sph_flavours.pdf
theory/paper_pasc/pasc_paper.pdf
m4/libtool.m4
diff --git a/theory/SPH/flavours/bibliography.bib b/theory/SPH/flavours/bibliography.bib
new file mode 100644
index 0000000000000000000000000000000000000000..365de80f3d112704516e29398d01828700377a3d
--- /dev/null
+++ b/theory/SPH/flavours/bibliography.bib
@@ -0,0 +1,82 @@
+@ARTICLE{Price2012,
+ author = {{Price}, D.~J.},
+ title = "{Smoothed particle hydrodynamics and magnetohydrodynamics}",
+ journal = {Journal of Computational Physics},
+archivePrefix = "arXiv",
+ eprint = {1012.1885},
+ primaryClass = "astro-ph.IM",
+ year = 2012,
+ month = feb,
+ volume = 231,
+ pages = {759-794},
+ doi = {10.1016/j.jcp.2010.12.011},
+ adsurl = {http://adsabs.harvard.edu/abs/2012JCoPh.231..759P},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+@ARTICLE{Springel2005,
+ author = {{Springel}, V.},
+ title = "{The cosmological simulation code GADGET-2}",
+ journal = {\mnras},
+ eprint = {astro-ph/0505010},
+ keywords = {methods: numerical, galaxies: interactions, dark matter},
+ year = 2005,
+ month = dec,
+ volume = 364,
+ pages = {1105-1134},
+ doi = {10.1111/j.1365-2966.2005.09655.x},
+ adsurl = {http://adsabs.harvard.edu/abs/2005MNRAS.364.1105S},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+
+@ARTICLE{Hopkins2013,
+ author = {{Hopkins}, P.~F.},
+ title = "{A general class of Lagrangian smoothed particle hydrodynamics methods and implications for fluid mixing problems}",
+ journal = {\mnras},
+archivePrefix = "arXiv",
+ eprint = {1206.5006},
+ primaryClass = "astro-ph.IM",
+ keywords = {hydrodynamics, instabilities, turbulence, methods: numerical, cosmology: theory},
+ year = 2013,
+ month = feb,
+ volume = 428,
+ pages = {2840-2856},
+ doi = {10.1093/mnras/sts210},
+ adsurl = {http://adsabs.harvard.edu/abs/2013MNRAS.428.2840H},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Springel2002,
+ author = {{Springel}, V. and {Hernquist}, L.},
+ title = "{Cosmological smoothed particle hydrodynamics simulations: the entropy equation}",
+ journal = {\mnras},
+ eprint = {astro-ph/0111016},
+ keywords = {methods: numerical, galaxies: evolution, galaxies: starburst, methods: numerical, galaxies: evolution, galaxies: starburst},
+ year = 2002,
+ month = jul,
+ volume = 333,
+ pages = {649-664},
+ doi = {10.1046/j.1365-8711.2002.05445.x},
+ adsurl = {http://adsabs.harvard.edu/abs/2002MNRAS.333..649S},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Balsara1995,
+ author = {{Balsara}, D.~S.},
+ title = "{von Neumann stability analysis of smooth particle hydrodynamics--suggestions for optimal algorithms}",
+ journal = {Journal of Computational Physics},
+ year = 1995,
+ volume = 121,
+ pages = {357-372},
+ doi = {10.1016/S0021-9991(95)90221-X},
+ adsurl = {http://adsabs.harvard.edu/abs/1995JCoPh.121..357B},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+
+
+
diff --git a/theory/SPH/flavours/run.sh b/theory/SPH/flavours/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b331e5d508f3ac900069dc8e99417796bf442231
--- /dev/null
+++ b/theory/SPH/flavours/run.sh
@@ -0,0 +1,5 @@
+#!/bin/bash
+pdflatex sph_flavours.tex
+bibtex sph_flavours.aux
+pdflatex sph_flavours.tex
+pdflatex sph_flavours.tex
diff --git a/theory/SPH/flavours/sph_flavours.tex b/theory/SPH/flavours/sph_flavours.tex
new file mode 100644
index 0000000000000000000000000000000000000000..26420ca85cd4c251c1ab372398b4d30f6b0606de
--- /dev/null
+++ b/theory/SPH/flavours/sph_flavours.tex
@@ -0,0 +1,295 @@
+\documentclass[fleqn, usenatbib, useAMS, a4paper]{mnras}
+\usepackage{graphicx}
+\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
+\usepackage{times}
+\usepackage[super]{nth}
+
+\newcommand{\swift}{{\sc swift}\xspace}
+\newcommand{\gadget}{{\sc gadget}\xspace}
+\newcommand{\dd}[2]{\frac{\partial #1}{\partial #2}}
+\renewcommand{\vec}[1]{{\mathbf{#1}}}
+\newcommand{\Wij}{\overline{\nabla_xW_{ij}}}
+
+%\setlength{\mathindent}{0pt}
+
+%opening
+\title{SPH flavours in SWIFT}
+\author{Matthieu Schaller}
+
+\begin{document}
+
+\maketitle
+
+For vector quantities, we define $\vec{a}_{ij} \equiv (\vec{a}_i -
+\vec{a}_j)$. Barred quantities ($\bar a_{ij}$) correspond to the
+average over two particles of a quantity: $\bar a_{ij} \equiv
+\frac{1}{2}(a_i + a_j)$. To simplify notations, we also define the
+vector quantity $\Wij \equiv \frac{1}{2}\left(W(\vec{x}_{ij}, h_i) +
+\nabla_x W(\vec{x}_{ij},h_j)\right)$.
+
+\subsection{Minimal SPH}
+\label{sec:sph:minimal}
+
+This is the simplest fully-conservative version of SPH using the
+internal energy $u$ as a thermal variable that can be written
+down. Its implementation in \swift should be understood as a text-book
+example and template for more advanced implementation. A full
+derivation of the equations can be found in the review of
+\cite{Price2012}. Our implementation follows their equations (27), (43),
+(44), (45), (101), (103) and (104) with $\beta=3$ and $\alpha_u=0$. We
+summarize the equations here.
+
+\subsubsection{Density and other fluid properties (\nth{1} neighbour loop)}
+
+For a set of particles $i$ with positions $\vec{x}_i$ with velocities
+$\vec{v}_i$, masses $m_i$, sthermal energy per unit mass $u_i$ and
+smoothing length $h_i$, we compute the density for each particle:
+
+\begin{equation}
+ \rho_i \equiv \rho(\vec{x}_i) = \sum_j m_j W(\vec{x}_{ij}, h_i),
+ \label{eq:sph:minimal:rho}
+\end{equation}
+and the derivative of its density with respect to $h$:
+
+\begin{equation}
+ \label{eq:sph:minimal:rho_dh}
+ \rho_{\partial h_i} \equiv \dd{\rho}{h}(\vec{x}_i) = \sum_j m_j \dd{W}{h}(\vec{x}_{ij}
+ , h_i).
+\end{equation}
+The gradient terms (``h-terms'') can then be computed from the density
+and its derivative:
+
+\begin{equation}
+ f_i \equiv \left(1 + \frac{h_i}{3\rho_i}\rho_{\partial h_i}
+ \right)^{-1}.
+ \label{eq:sph:minimal:f_i}
+\end{equation}
+Using the pre-defined equation of state, the pressure $P_i$ and the sound
+speed $c_i$ at the location of particle $i$ can now be computed from
+$\rho_i$ and $u_i$:
+
+\begin{align}
+ P_i &= P_{\rm eos}(\rho_i, u_i), \label{eq:sph:minimal:P}\\
+ c_i &= c_{\rm eos}(\rho_i, u_i). \label{eq:sph:minimal:c}
+\end{align}
+
+\subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
+
+We can then proceed with the second loop over
+neighbours. The signal velocity $v_{{\rm sig},ij}$ between two particles is given by
+
+\begin{align}
+ \mu_{ij} &=
+ \begin{cases}
+ \frac{\vec{v}_{ij} \cdot \vec{x}_{ij}}{|\vec{x}_{ij}|} & \rm{if}~
+ \vec{v}_{ij} \cdot \vec{x}_{ij} < 0,\\
+ 0 &\rm{otherwise}, \\
+ \end{cases}\nonumber\\
+ v_{{\rm sig},ij} &= c_i + c_j - 3\mu_{ij}. \label{eq:sph:minimal:v_sig}
+\end{align}
+We also use these two quantities for the calculation of the viscosity term:
+
+\begin{equation}
+\nu_{ij} = -\frac{1}{2}\frac{\alpha \mu_{ij} v_{{\rm
+ sig},ij}}{\bar\rho_{ij}}
+ \label{eq:sph:minimal:nu_ij}
+\end{equation}
+The fluid accelerations are then given by
+
+\begin{align}
+ \frac{d\vec{v}_i}{dt} = -\sum_j m_j &\left[\frac{f_iP_i}{\rho_i^2}
+ \nabla_x W(\vec{x}_{ij}, h_i) \nonumber
+ +\frac{f_jP_j}{\rho_j^2}\nabla_x W(\vec{x}_{ij},h_j)\right. \\
+ &+ \left. \bigg.\nu_{ij} \Wij \right], \label{eq:sph:minimal:dv_dt}
+\end{align}
+and the change in internal energy,
+
+\begin{align}
+ \frac{du_i}{dt} = \sum_j m_j &\left[\frac{f_iP_i}{\rho_i^2} \vec{v}_{ij}
+ \cdot \nabla_x W(\vec{x}_{ij}, h_i) \right. \label{eq:sph:minimal:du_dt}\\
+ &+\left. \frac{1}{2}\nu_{ij}\vec{v}_{ij}\cdot\Big. \Wij\right], \nonumber
+\end{align}
+where in both cases the first line corresponds to the standard SPH
+term and the second line to the viscuous accelerations.
+
+We also compute an estimator of the change in smoothing length to be
+used in the prediction step. This is an estimate of the local
+divergence of the velocity field compatible with the accelerations
+computed above:
+
+\begin{equation}
+ \frac{dh_i}{dt} = -\frac{1}{3}h_i \sum_j \frac{m_j}{\rho_j}
+ \vec{v}_{ij}\cdot \nabla_x W(\vec{x}_{ij}, h_i).
+ \label{eq:sph:minimal:dh_dt}
+\end{equation}
+and update the signal velocity of the particles:
+
+\begin{equation}
+ v_{{\rm sig},i} = \max_j \left( v_{{\rm sig},ij} \right).
+ \label{eq:sph:minimal:v_sig_update}
+\end{equation}
+All the quantities required for time integration have now been obtained.
+
+\subsubsection{Time integration}
+
+For each particle, we compute a time-step given by the CFL condition:
+
+\begin{equation}
+ \Delta t = 2 C_{\rm CFL} \frac{H_i}{v_{{\rm sig},i}},
+ \label{eq:sph:minimal:dt}
+\end{equation}
+where $C_{\rm CFL}$ is a free dimensionless parameter and $H_i = \gamma h_i$ is the
+kernel support size. Particles can then be ``kicked'' forward in time:
+\begin{align}
+ \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t \label{eq:sph:minimal:kick_v}\\
+ u_i &\rightarrow u_i + \frac{du_i}{dt} \Delta t \label{eq:sph:minimal:kick_u}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i\right) \label{eq:sph:minimal:kick_P}, \\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i, u_i\right) \label{eq:sph:minimal:kick_c},
+\end{align}
+where we used the pre-defined equation of state to compute the new
+value of the pressure and sound-speed.
+
+\subsubsection{Particle properties prediction}
+
+Inactive particles need to have their quantities predicted forward in
+time in the ``drift'' operation. We update them as follows:
+
+\begin{align}
+ h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
+ \Delta t\right), \label{eq:sph:minimal:drift_h}\\
+ \rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
+ \Delta t\right), \label{eq:sph:minimal:drift_rho} \\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i + \frac{du_i}{dt} \Delta t\right), \label{eq:sph:minimal:drift_P}\\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i, u_i + \frac{du_i}{dt}
+ \Delta t\right) \label{eq:sph:minimal:drift_c},
+\end{align}
+where, as above, the last two updated quantities are obtained using
+the pre-defined equation of state. Note that the thermal energy $u_i$
+itself is \emph{not} updated.
+
+\subsection{Gadget-2 SPH}
+\label{sec:sph:gadget2}
+
+This flavour of SPH is the one implemented in the \gadget-2 code
+\citep{Springel2005}. The basic equations were derived by
+\cite{Springel2002} and also includes a \cite{Balsara1995} switch for
+the suppression of viscosity. The implementation here follows closely the
+presentation of \cite{Springel2005}. Specifically, we use their equations (5), (7),
+(8), (9), (10), (13), (14) and (17). We summarize them here for completeness.
+
+\subsubsection{Density and other fluid properties (\nth{1} neighbour loop)}
+
+For a set of particles $i$ with positions $\vec{x}_i$ with velocities
+$\vec{v}_i$, masses $m_i$, entropic function per unit mass $A_i$ and
+smoothing length $h_i$, we compute the density, derivative of the density with respect
+to $h$ and the ``f-terms'' in a similar way to the minimal-SPH case
+(Eq. \ref{eq:sph:minimal:rho}, \ref{eq:sph:minimal:rho_dh} and
+\ref{eq:sph:minimal:f_i}). From these the pressure and sound-speed can
+be computed using the predefined equation of state:
+
+\begin{align}
+ P_i &= P_{\rm eos}(\rho_i, A_i), \label{eq:sph:gadget2:P}\\
+ c_i &= c_{\rm eos}(\rho_i, A_i). \label{eq:sph:gadget2:c}
+\end{align}
+We additionally compute the divergence and
+curl of the velocity field using standard SPH expressions:
+
+\begin{align}
+ \nabla\cdot\vec{v}_i \equiv\nabla\cdot \vec{v}(\vec{x}_i) &= \frac{1}{\rho_i} \sum_j m_j
+ \vec{v}_{ij}\cdot\nabla_x W(\vec{x}_{ij}, h_i) \label{eq:sph:gadget2:div_v},\\
+ \nabla\times\vec{v}_i \equiv \nabla\times \vec{v}(\vec{x}_i) &= \frac{1}{\rho_i} \sum_j m_j
+ \vec{v}_{ij}\times\nabla_x W(\vec{x}_{ij}, h_i) \label{eq:sph:gadget2:rot_v}.
+\end{align}
+These are used to construct the \cite{Balsara1995} switch $B_i$:
+
+\begin{equation}
+ B_i = \frac{|\nabla\cdot\vec{v}_i|}{|\nabla\cdot\vec{v}_i| +
+ |\nabla\times\vec{v}_i| + 10^{-4}c_i / h_i}, \label{eq:sph:gadget2:balsara}
+\end{equation}
+where the last term in the denominator is added to prevent numerical instabilities.
+
+\subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
+
+The accelerations are computed in a similar fashion to the minimal-SPH
+case with the exception of the viscosity term which is now modified to
+include the switch. Instead of Eq. \ref{eq:sph:minimal:nu_ij}, we get:
+
+\begin{equation}
+\nu_{ij} = -\frac{1}{2}\frac{\alpha \bar B_{ij} \mu_{ij} v_{{\rm sig},ij}}{\bar\rho_{ij}},
+ \label{eq:sph:gadget2:nu_ij}
+\end{equation}
+whilst equations \ref{eq:sph:minimal:v_sig},
+\ref{eq:sph:minimal:dv_dt}, \ref{eq:sph:minimal:dh_dt} and
+\ref{eq:sph:minimal:v_sig_update} remain unchanged. The only other
+change is the equation of motion for the thermodynamic variable which
+now has to be describing the evolution of the entropic function and
+not the evolution of the thermal energy. It reads:
+
+\begin{equation}
+\frac{dA_i}{dt} = \frac{1}{2} \frac{\gamma-1}{\rho^{\gamma-1}} \sum_j
+m_j \nu_{ij}\vec{v}_{ij}\cdot \Wij
+\end{equation}
+
+\subsubsection{Time integration}
+
+The time-step condition is identical to the Minimal-SPH case
+(Eq. \ref{eq:sph:minimal:dt}). The same applies to the integration
+forward in time (Eq. \ref{eq:sph:minimal:kick_v} to
+\ref{eq:sph:minimal:kick_c}) with the exception of the change in
+internal energy (Eq. \ref{eq:sph:minimal:kick_u}) which gets replaced
+by an integration for the the entropy:
+
+
+\begin{align}
+ \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t \label{eq:sph:gadget2:kick_v}\\
+ A_i &\rightarrow A_i + \frac{dA_i}{dt} \Delta t \label{eq:sph:gadget2:kick_A}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:gadget2:kick_P}, \\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:gadget2:kick_c},
+\end{align}
+where, once again, we made use of the equation of state relating
+thermodynaimcal quantities.
+
+\subsubsection{Particle properties prediction}
+
+The prediction step is also identical to the Minimal-SPH case with the
+entropic function replacing the thermal energy.
+
+\begin{align}
+ h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
+ \Delta t\right), \label{eq:sph:gadget2:drift_h}\\
+ \rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
+ \Delta t\right), \label{eq:sph:gadget2:drift_rho} \\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt} \Delta t\right), \label{eq:sph:gadget2:drift_P}\\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt}
+ \Delta t\right) \label{eq:sph:gadget2:drift_c},
+\end{align}
+where, as above, the last two updated quantities are obtained using
+the pre-defined equation of state. Note that the entropic function $A_i$
+itself is \emph{not} updated.
+
+
+
+\subsection{Pressure-Entropy SPH}
+
+\cite{Hopkins2013}
+\label{sec:sph:pe}
+
+\subsubsection{Hydrodynamical accelerations}
+
+\subsubsection{Particle properties prediction}
+
+\subsection{Pressure-Energy SPH}
+\label{sec:sph:pu}
+
+\cite{Hopkins2013}\\
+to be done
+
+\subsection{Anarchy SPH}
+\label{sec:sph:anarchy}
+
+to be done
+
+\bibliographystyle{mnras}
+\bibliography{./bibliography.bib}
+
+\end{document}
diff --git a/theory/SPH/kernel/kernel_definitions.tex b/theory/SPH/kernel/kernel_definitions.tex
index c689f3da01228f4655f24feece961b434374b6a3..b29d23ba58775bc0a1bb2a71ce49fad7b746502c 100644
--- a/theory/SPH/kernel/kernel_definitions.tex
+++ b/theory/SPH/kernel/kernel_definitions.tex
@@ -1,4 +1,4 @@
-\documentclass[usenatbib, useAMS,a4paper]{mnras}
+\documentclass[fleqn, usenatbib, useAMS,a4paper]{mnras}
\usepackage{graphicx}
\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
\usepackage{times}
@@ -27,14 +27,12 @@ The desirable properties of an SPH kernels $W(\vec{x},h)$ are:
\item $W(\vec{x},h)$ should have a finite support and be cheap to
compute.
\end{enumerate}
-
As a consequence, the smoothing kernels used in SPH can
hence be written (in 3D) as
\begin{equation}
W(\vec{x},h) \equiv \frac{1}{H^3}f\left(\frac{|\vec{x}|}{H}\right),
\end{equation}
-
where $H=\gamma h$ is defined below and $f(u)$ is a dimensionless
function, usually a low-order polynomial, such that $f(u \geq 1) = 0$
and normalised such that
@@ -42,7 +40,6 @@ and normalised such that
\begin{equation}
\int f(|\vec{u}|){\rm d}^3u = 1.
\end{equation}
-
$H$ is the kernel's support radius and is used as the ``smoothing
length'' in the Gadget code( {i.e.} $H=h$). This definition is,
however, not very physical and makes comparison of kernels at a
@@ -55,13 +52,11 @@ standard deviation
\sigma^2 \equiv \frac{1}{3}\int \vec{u}^2 W(\vec{u},h) {\rm d}^3u.
\label{eq:sph:sigma}
\end{equation}
-
The smoothing length is then:
\begin{equation}
h\equiv2\sigma.
\label{eq:sph:h}
\end{equation}
-
Each kernel, {\it i.e.} defintion of $f(u)$, will have a different
ratio $\gamma = H/h$. So for a \emph{fixed resolution} $h$, one will
have different kernel support sizes, $H$, and a different number of
@@ -72,18 +67,15 @@ separation:
\begin{equation}
h = \eta \langle x \rangle = \eta \left(\frac{m}{\rho}\right)^{1/3},
\end{equation}
-
where $\rho$ is the local density of the fluid and $m$ the SPH
particle mass.
-
The (weighted) number of neighbours within the kernel support is a
-useful quantity to use in implementations of SPH. It is defined as (in
-3D)
+useful quantity to use in implementations of SPH. It is defined (in
+3D) as:
\begin{equation}
N_{\rm ngb} \equiv \frac{4}{3}\pi \left(\frac{H}{h}\eta\right)^3.
\end{equation}
-
Once the fixed ratio $\gamma= H/h$ is known (via equations
\ref{eq:sph:sigma} and \ref{eq:sph:h}) for a given kernel, the number
of neighbours only depends on the resolution parameter $\eta$. For
@@ -96,23 +88,22 @@ resolution to $\eta=1.2348$ yields the commonly used value $N_{\rm
The \swift kernels are split into two categories, the B-splines
($M_{4,5,6}$) and the Wendland kernels ($C2$, $C4$ and $C6$). In all
cases we impose $f(u>1) = 0$.\\
-
The spline kernels are defined as:
\begin{align}
- f(u) &= C M_n(u), \\
+ f(u) &= C M_n(u), \nonumber\\
M_n(u) &\equiv \frac{1}{2\pi}
\int_{-\infty}^{\infty}
\left(\frac{\sin\left(k/n\right)}{k/n}\right)^n\cos\left(ku\right){\rm
- d}k,
+ d}k \nonumber,
\end{align}
-
-whilst the Wendland kernels read
+whilst the Wendland kernels are constructed from:
\begin{align}
- f(u) &= C \Psi_{i,j}(u), \\
- \Psi_{i,j}(u) &\equiv \mathcal{I}^k\left[\max\left(1-u,0\right)\right],\\
- \mathcal{I}[f](u) &\equiv \int_u^\infty f\left(k\right)k{\rm d}k.
+ f(u) &= C \Psi_{i,j}(u), \nonumber\\
+ \Psi_{i,j}(u) &\equiv
+ \mathcal{I}^k\left[\max\left(1-u,0\right)\right], \nonumber\\
+ \mathcal{I}[f](u) &\equiv \int_u^\infty f\left(k\right)k{\rm d}k. \nonumber
\end{align}
\subsubsection{Cubic spline ($M_4$) kernel}
@@ -230,11 +221,11 @@ All kernels available in \swift are shown on Fig.~\ref{fig:sph:kernels}.
The derivatives of the kernel function have relatively simple
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}|}, \\
- \frac{\partial W(\vec{x},h)}{\partial h} &=&- \frac{1}{h^4}\left[3f\left(\frac{|\vec{x}|}{h}\right) +
+\begin{align}
+ \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].
-\end{eqnarray*}
+\end{align}
Note that for all the kernels listed above, $f'(0) = 0$.
\bibliographystyle{mnras}