diff --git a/.gitignore b/.gitignore
index a8778a3b82b2cf89f2169f2580d16b0975bb41ba..0c4d8bddadb17afcc29aa3dc1d4dba6cb6fd2f5b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -102,6 +102,9 @@ theory/paper_pasc/pasc_paper.pdf
 theory/Multipoles/fmm.pdf
 theory/Multipoles/fmm_standalone.pdf
 theory/Multipoles/potential.pdf
+theory/Multipoles/potential_long.pdf
+theory/Multipoles/potential_short.pdf
+theory/Multipoles/force_short.pdf
 
 m4/libtool.m4
 m4/ltoptions.m4
diff --git a/src/engine.c b/src/engine.c
index 6bd2fb5ade526e4ff7fd363fa655fd096afcd827..89bb6c693e63b2705e1b8a9371713edc0945acb1 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -909,6 +909,9 @@ void engine_repartition(struct engine *e) {
 #else
   if (e->reparttype->type != REPART_NONE)
     error("SWIFT was not compiled with MPI and METIS support.");
+
+  /* Clear the repartition flag. */
+  e->forcerepart = 0;
 #endif
 }
 
@@ -923,8 +926,9 @@ void engine_repartition_trigger(struct engine *e) {
 
   /* Do nothing if there have not been enough steps since the last
    * repartition, don't want to repeat this too often or immediately after
-   * a repartition step. */
-  if (e->step - e->last_repartition >= 2) {
+   * a repartition step. Also nothing to do when requested. */
+  if (e->step - e->last_repartition >= 2 &&
+      e->reparttype->type != REPART_NONE) {
 
     /* Old style if trigger is >1 or this is the second step (want an early
      * repartition following the initial repartition). */
@@ -985,8 +989,9 @@ void engine_repartition_trigger(struct engine *e) {
     if (e->forcerepart) e->last_repartition = e->step;
   }
 
-  /* We always reset CPU time for next check. */
-  e->cputime_last_step = clocks_get_cputime_used();
+  /* We always reset CPU time for next check, unless it will not be used. */
+  if (e->reparttype->type != REPART_NONE)
+    e->cputime_last_step = clocks_get_cputime_used();
 #endif
 }
 
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 276d6d837bcc164bbc006906b2e632041e033cae..1e7554f7d84220b8c962d60cc4538c685b5bad52 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -94,7 +94,7 @@ void hydro_props_print(const struct hydro_props *p) {
   message("Hydrodynamic kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
           p->eta_neighbours, p->target_neighbours);
 
-  message("Hydrodynamic tolerance in h: %.5f (+/- %.4f neighbours).",
+  message("Hydrodynamic relative tolerance in h: %.5f (+/- %.4f neighbours).",
           p->h_tolerance, p->delta_neighbours);
 
   message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
diff --git a/theory/Multipoles/cells.pdf b/theory/Multipoles/cells.pdf
index b1c144ce83e44b42f3fb7dab48b675cca19a288e..d621f6f1023d71503f698b69694d980ef27814e6 100644
Binary files a/theory/Multipoles/cells.pdf and b/theory/Multipoles/cells.pdf differ
diff --git a/theory/Multipoles/fmm_summary.tex b/theory/Multipoles/fmm_summary.tex
index dd191d6c99da9f73c8fcac40facf4425adcbff8c..1ff9ab88ada6836d6118c7cfd74e39f4d1c504b3 100644
--- a/theory/Multipoles/fmm_summary.tex
+++ b/theory/Multipoles/fmm_summary.tex
@@ -6,7 +6,7 @@ evaluate for each particle in a system the potential and associated
 forces generated by all the other particles. Mathematically, this means
 evaluate
 \begin{equation}
-  \phi(\mathbf{x}_a) = \sum_{b \neq a} G m_b\phi(\mathbf{x}_a -
+  \phi(\mathbf{x}_a) = \sum_{b \neq a} G m_b\varphi(\mathbf{x}_a -
   \mathbf{x}_b)\qquad \forall~a\in N
   \label{eq:fmm:n_body}
 \end{equation}
@@ -23,10 +23,11 @@ arising from nearby particles. \\
 
 In what follows, we use the compact multi-index notation of
 \cite{Dehnen2014} (repeated in appendix \ref{sec:multi_index_notation}
-for completeness) to simplify expressions. $\mathbf{k}$, $\mathbf{m}$
-and $\mathbf{n}$ are multi-indices and $\mathbf{r}$, $\mathbf{R}$,
-$\mathbf{x}$, $\mathbf{y}$ and $\mathbf{z}$ are vectors, whilst $a$
-and $b$ are particle indices.\\
+for completeness) to simplify expressions and ease
+comparisons. $\mathbf{k}$, $\mathbf{m}$ and $\mathbf{n}$ are
+multi-indices and $\mathbf{r}$, $\mathbf{R}$, $\mathbf{x}$,
+$\mathbf{y}$ and $\mathbf{z}$ are vectors, whilst $a$ and $b$ are
+particle indices.\\
 
 \begin{figure}
 \includegraphics[width=\columnwidth]{cells.pdf}
@@ -46,34 +47,34 @@ with centres of mass $\mathbf{z}_A$ and  $\mathbf{z}_B$
 respectively, as shown on Fig.~\ref{fig:fmm:cells}, the potential
 generated by $b$ at the location of $a$ can be rewritten as
 \begin{align}
-  \phi(\mathbf{x}_a - \mathbf{x}_b)
-  &= \phi\left(\mathbf{x}_a - \mathbf{z}_A - \mathbf{x}_b +
+  \varphi(\mathbf{x}_a - \mathbf{x}_b)
+  &= \varphi\left(\mathbf{x}_a - \mathbf{z}_A - \mathbf{x}_b +
   \mathbf{z}_B + \mathbf{z}_A - \mathbf{z}_B\right)  \nonumber \\
-  &= \phi\left(\mathbf{r}_a - \mathbf{r}_b + \mathbf{R}\right)
+  &= \varphi\left(\mathbf{r}_a - \mathbf{r}_b + \mathbf{R}\right)
   \nonumber \\
   &= \sum_\mathbf{k} \frac{1}{\mathbf{k}!} \left(\mathbf{r}_a -
-  \mathbf{r}_b\right)^{\mathbf{k}} \nabla^{\mathbf{k}}\phi(\mathbf{R})
+  \mathbf{r}_b\right)^{\mathbf{k}} \nabla^{\mathbf{k}}\varphi(\mathbf{R})
   \nonumber \\
   &= \sum_\mathbf{k} \frac{1}{\mathbf{k}!} \sum_{\mathbf{n} <
     \mathbf{k}} \binom{\mathbf{k}}{\mathbf{n}} \mathbf{r}_a^{\mathbf{n}}
   \left(-\mathbf{r}_b\right)^{\mathbf{k} - \mathbf{n}}
-  \nabla^{\mathbf{k}}\phi(\mathbf{R})\nonumber \\
+  \nabla^{\mathbf{k}}\varphi(\mathbf{R})\nonumber \\
   &= \sum_\mathbf{n} \frac{1}{\mathbf{n}!} \mathbf{r}_a^{\mathbf{n}}
   \sum_\mathbf{m} \frac{1}{\mathbf{m}!}
-  \left(-\mathbf{r}_b\right)^\mathbf{m} \nabla^{\mathbf{n}+\mathbf{m}} \phi(\mathbf{R}),
+  \left(-\mathbf{r}_b\right)^\mathbf{m} \nabla^{\mathbf{n}+\mathbf{m}} \varphi(\mathbf{R}),
 \end{align}
-where we used the Taylor expansion of $\phi$ around $\mathbf{R} \equiv
+where we used the Taylor expansion of $\varphi$ around $\mathbf{R} \equiv
 \mathbf{z}_A - \mathbf{z}_B$ on the third line, used $\mathbf{r}_a
 \equiv \mathbf{x}_a - \mathbf{z}_A$, $\mathbf{r}_b \equiv \mathbf{x}_b
 - \mathbf{z}_B$ throughout and defined $\mathbf{m} \equiv
 \mathbf{k}-\mathbf{n}$ on the last line. Expanding the series only up
 to order $p$, we get
 \begin{equation}
-  \phi(\mathbf{x}_a - \mathbf{x}_b) \approx \sum_{\mathbf{n}}^{p}
+  \varphi(\mathbf{x}_a - \mathbf{x}_b) \approx \sum_{\mathbf{n}}^{p}
   \frac{1}{\mathbf{n}!} \mathbf{r}_a^{\mathbf{n}} \sum_{\mathbf{m}}^{p
     -|\mathbf{n}|} 
   \frac{1}{\mathbf{m}!} \left(-\mathbf{r}_b\right)^\mathbf{m}
-  \nabla^{\mathbf{n}+\mathbf{m}} \phi(\mathbf{R}),
+  \nabla^{\mathbf{n}+\mathbf{m}} \varphi(\mathbf{R}),
   \label{eq:fmm:fmm_one_part}
 \end{equation}
 with the approximation converging as $p\rightarrow\infty$ towards the
@@ -82,13 +83,13 @@ correct value provided $|\mathbf{R}|<|\mathbf{r}_a +
 combine their contributions to the potential at location
 $\mathbf{x}_a$ in cell $A$, we get
 \begin{align}
-  \phi_{BA}(\mathbf{x}_a) &= \sum_{b\in B} m_b\phi(\mathbf{x}_a -
+  \phi_{BA}(\mathbf{x}_a) &= \sum_{b\in B}G m_b\varphi(\mathbf{x}_a -
   \mathbf{x}_b)  \label{eq:fmm:fmm_one_cell}  \\
-  &\approx \sum_{\mathbf{n}}^{p}
+  &\approx G\sum_{\mathbf{n}}^{p}
   \frac{1}{\mathbf{n}!} \mathbf{r}_a^{\mathbf{n}} \sum_{\mathbf{m}}
     ^{p -|\mathbf{n}|}
   \frac{1}{\mathbf{m}!} \sum_{b\in B} m_b\left(-\mathbf{r}_b\right)^\mathbf{m}
-  \nabla^{\mathbf{n}+\mathbf{m}} \phi(\mathbf{R}) \nonumber. 
+  \nabla^{\mathbf{n}+\mathbf{m}} \varphi(\mathbf{R}) \nonumber. 
 \end{align}
 This last equation forms the basis of the FMM. The algorithm
 decomposes the equation into three separated sums evaluated at
@@ -106,12 +107,12 @@ compute the second kernel, \textsc{M2L} (multipole to local
 expansion), which corresponds to the interaction of a cell with
 another one:
 \begin{equation}
-  \mathsf{F}_{\mathbf{n}}(\mathbf{z}_A) = \sum_{\mathbf{m}}^{p -|\mathbf{n}|}
+  \mathsf{F}_{\mathbf{n}}(\mathbf{z}_A) = G\sum_{\mathbf{m}}^{p -|\mathbf{n}|}
   \mathsf{M}_{\mathbf{m}}(\mathbf{z}_B)
   \mathsf{D}_{\mathbf{n}+\mathbf{m}}(\mathbf{R}), \label{eq:fmm:M2L} 
 \end{equation}
 where $\mathsf{D}_{\mathbf{n}+\mathbf{m}}(\mathbf{R}) \equiv
-\nabla^{\mathbf{n}+\mathbf{m}} \phi(\mathbf{R})$ is an order $n+m$
+\nabla^{\mathbf{n}+\mathbf{m}} \varphi(\mathbf{R})$ is an order $n+m$
 derivative of the potential. This is the computationally expensive
 step of the FMM algorithm as the number of operations in a naive
 implementation using cartesian coordinates scales as
@@ -165,8 +166,8 @@ multiplications (provided $\mathsf{D}$ can be evaluated quickly, see
 Sec.~\ref{ssec:grav_derivatives}), which are extremely efficient
 instructions on modern architectures. However, the fully expanded sums
 can lead to rather large and prone to typo expressions. To avoid any
-misshap, we use a \texttt{python} to generate C code in which all the
-sums are unrolled and correct by construction. In \swift, we
+mishaps, we use a \texttt{python} script to generate C code in which
+all the sums are unrolled and correct by construction. In \swift, we
 implemented the kernels up to order $p=5$, as it proved to be accurate
 enough for our purpose, but this could be extended to higher order
 easily. This implies storing $56$ numbers per cell for each
diff --git a/theory/Multipoles/gravity_derivatives.tex b/theory/Multipoles/gravity_derivatives.tex
index 4d9cf6c918c2214ddce432a0067cf4e2f9aaf682..e4569ef960fae5e92343f1d99902a5c14fd6ee5c 100644
--- a/theory/Multipoles/gravity_derivatives.tex
+++ b/theory/Multipoles/gravity_derivatives.tex
@@ -2,36 +2,36 @@
 \label{ssec:grav_derivatives}
 
 The calculation of all the
-$\mathsf{D}_\mathbf{n}(x,y,z) \equiv \nabla^{\mathbf{n}}\phi(x,y,z)$ terms up
+$\mathsf{D}_\mathbf{n}(x,y,z) \equiv \nabla^{\mathbf{n}}\varphi(x,y,z)$ terms up
 to the relevent order can be quite tedious and it is beneficial to
 automatize the whole setup. Ideally, one would like to have an
 expression for each of these terms that is only made of multiplications
 and additions of each of the coordinates and the inverse distance. We
-achieve this by writing $\phi$ as a composition of functions
-$\phi(u(x,y,z))$ and apply the \textit{Fa\`a di Bruno}
+achieve this by writing $\varphi$ as a composition of functions
+$\varphi(u(x,y,z))$ and apply the \textit{Fa\`a di Bruno}
 formula \citep[i.e. the ``chain rule'' for higher order derivatives,
  see e.g.][]{Hardy2006} to construct our terms:
 \begin{equation}
 \label{eq:faa_di_bruno}
-\frac{\partial^n}{\partial x_1 \cdots \partial x_n} \phi(u)
-= \sum_{A} \phi^{(|A|)}(u) \prod_{B \in
+\frac{\partial^n}{\partial x_1 \cdots \partial x_n} \varphi(u)
+= \sum_{A} \varphi^{(|A|)}(u) \prod_{B \in
 A} \frac{\partial^{|B|}}{\prod_{c\in B}\partial x_c} u(x,y,z),
 \end{equation}
 where $A$ is the set of all partitions of $\lbrace1,\cdots, n\rbrace$,
 $B$ is a block of a partition in the set $A$ and $|\cdot|$ denotes the
-cardinality of a set. For generic functions $\phi$ and $u$ this
+cardinality of a set. For generic functions $\varphi$ and $u$ this
 formula yields an untracktable number of terms; an 8th-order
 derivative will have $4140$ (!)  terms in the sum\footnote{The number
   of terms in the sum is given by the Bell number of the same
   order.}. \\ For the un-softened gravitational potential, we choose to write
 \begin{align}
-   \phi(x,y,z) &= 1 / \sqrt{u(x,y,z)}, \\
+   \varphi(x,y,z) &= 1 / \sqrt{u(x,y,z)}, \\
    u(x,y,z) &= x^2 + y^2 + z^2.
 \end{align}
-This choice allows to have derivatives of any order of $\phi(u)$ that
+This choice allows to have derivatives of any order of $\varphi(u)$ that
 can be easily expressed and only depend on powers of $u$:
 \begin{equation}
-\phi^{(n)}(u) = (-1)^n\cdot\frac{(2n-1)!!}{2^n}\cdot\frac{1}{u^{n+\frac{1}{2}}},
+\varphi^{(n)}(u) = (-1)^n\cdot\frac{(2n-1)!!}{2^n}\cdot\frac{1}{u^{n+\frac{1}{2}}},
 \end{equation}
 where $!!$ denotes the semi-factorial. More importantly, this
 choice of decomposition allows us to have non-zero derivatives of $u$
@@ -39,7 +39,7 @@ only up to second order in $x$, $y$ or $z$. The number of non-zero
 terms in eq. \ref{eq:faa_di_bruno} is hence drastically reduced. For
 instance, when computing $\mathsf{D}_{(4,1,3)}(\mathbf{r}) \equiv
 \frac{\partial^8}{\partial x^4 \partial y \partial z^3}
-\phi(u(x,y,z))$, $4100$ of the $4140$ terms will involve at least one
+\varphi(u(x,y,z))$, $4100$ of the $4140$ terms will involve at least one
 zero-valued derivative (e.g. $\partial^3/\partial x^3$ or
 $\partial^2/\partial x\partial y$) of $u$. Furthermore, among the 40
 remaining terms, many will involve the same combination of derivatives
diff --git a/theory/Multipoles/mesh_summary.tex b/theory/Multipoles/mesh_summary.tex
index c9e8bc4577e1040ff7b6daadea25736ad6508ff9..3069257c8845804d9a307cc54fffec5e36e4ae8c 100644
--- a/theory/Multipoles/mesh_summary.tex
+++ b/theory/Multipoles/mesh_summary.tex
@@ -1,2 +1,39 @@
 \subsection{Coupling the FMM to a mesh for periodic long-range forces}
 \label{ssec:mesh_summary}
+
+\begin{equation}
+  S(x) = \frac{e^x}{1 + e^x}
+\end{equation}
+
+\begin{align}
+  \varphi_s(r) &= \frac{1}{r}\left[2 - 2S\left(\frac{2r}{r_s}\right)\right] \nonumber\\
+  &= \frac{1}{r}\left[2 - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}}\right] 
+\end{align}
+\begin{align}
+  |\mathbf{f}_s(r)| &= \frac{1}{r^2}\left[\frac{4r}{r_s}S'\left(\frac{2r}{r_s}\right) - 2S\left(\frac{2r}{r_s}\right) + 2\right] \nonumber \\
+  &= \frac{1}{r^2}\left[\frac{4r}{r_s}\frac{e^{\frac{2r}{r_s}}}{(1+e^{\frac{2r}{r_s}})^2} - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}} + 2\right]
+\end{align}
+
+\begin{equation}
+  \tilde\varphi_l(k) = \frac{1}{k^2}\left[\frac{\upi}{2}kr_s\textrm{csch}\left(\frac{\upi}{2}kr_s\right) \right]
+\end{equation}
+
+\begin{figure}
+\includegraphics[width=\columnwidth]{potential_short.pdf}
+\caption{aa}
+\label{fig:fmm:potential_short}
+\end{figure}
+
+
+\begin{figure}
+\includegraphics[width=\columnwidth]{force_short.pdf}
+\caption{bb}
+\label{fig:fmm:force_short}
+\end{figure}
+
+
+\begin{figure}
+\includegraphics[width=\columnwidth]{potential_long.pdf}
+\caption{cc}
+\label{fig:fmm:potential_long}
+\end{figure}
diff --git a/theory/Multipoles/plot_mesh.py b/theory/Multipoles/plot_mesh.py
new file mode 100644
index 0000000000000000000000000000000000000000..6706016f73b4b6251c6d517ec89eacbb7a469417
--- /dev/null
+++ b/theory/Multipoles/plot_mesh.py
@@ -0,0 +1,267 @@
+###############################################################################
+ # 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 *
+from scipy import integrate
+from scipy import special
+from scipy.optimize import curve_fit
+from scipy.optimize import fsolve
+from matplotlib.font_manager import FontProperties
+import numpy
+import math
+
+params = {'axes.labelsize': 9,
+'axes.titlesize': 10,
+'font.size': 10,
+'legend.fontsize': 10,
+'xtick.labelsize': 8,
+'ytick.labelsize': 8,
+'text.usetex': True,
+'figure.figsize' : (3.15,3.15),
+'figure.subplot.left'    : 0.12,
+'figure.subplot.right'   : 0.99  ,
+'figure.subplot.bottom'  : 0.09  ,
+'figure.subplot.top'     : 0.99  ,
+'figure.subplot.wspace'  : 0.  ,
+'figure.subplot.hspace'  : 0.  ,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+colors=['#4477AA', '#CC6677', '#DDCC77', '#117733']
+
+
+# Parameters
+r_s = 2.
+r_min = 1e-2
+r_max = 1.5e2
+
+# Radius
+r = logspace(log10(r_min), log10(r_max), 401)
+r_rs = r / r_s
+
+k = logspace(log10(r_min/r_s**2), log10(r_max/r_s**2), 401)
+k_rs = k * r_s
+
+# Newtonian solution
+phi_newton = 1. / r
+phit_newton = 1. / k**2
+force_newton = 1. / r**2
+
+def my_exp(x):
+    return 1. + x + (x**2 / 2.) + (x**3 / 6.) + (x**4 / 24.) + (x**5 / 120.)# + (x**6 / 720.)
+    #return exp(x)
+
+def csch(x): # hyperbolic cosecant
+    return 1. / sinh(x)
+
+def sigmoid(x):
+    return my_exp(x) / (my_exp(x) + 1.)
+
+def d_sigmoid(x):
+    return my_exp(x) / ((my_exp(x) + 1)**2)
+
+def swift_corr(x):
+    return 2 * sigmoid( 4 * x ) - 1
+
+#figure()
+#x = linspace(-4, 4, 100)
+#plot(x, special.erf(x), '-', color=colors[0])
+#plot(x, swift_corr(x), '-', color=colors[1])
+#plot(x, x, '-', color=colors[2])
+#ylim(-1.1, 1.1)
+#xlim(-4.1, 4.1)
+#savefig("temp.pdf")
+
+# Correction in real space
+corr_short_gadget2 = special.erf(r / (2.*r_s))
+corr_short_swift = swift_corr(r / (2.*r_s)) 
+eta_short_gadget2 = special.erfc(r / 2.*r_s) + (r / (r_s * math.sqrt(math.pi))) * exp(-r**2 / (4.*r_s**2))
+eta_short_swift = 4. * (r / r_s) * d_sigmoid(2. * r / r_s) - 2. * sigmoid(2 * r / r_s) + 2.
+
+# Corection in Fourier space
+corr_long_gadget2 = exp(-k**2*r_s**2)
+corr_long_swift = math.pi * k * r_s * csch(0.5 * math.pi * r_s * k) / 2.
+
+# Shortrange term
+phi_short_gadget2 = (1.  / r ) * (1. - corr_short_gadget2)
+phi_short_swift = (1.  / r ) * (1. - corr_short_swift)
+force_short_gadget2 = (1. / r**2) * eta_short_gadget2
+force_short_swift = (1. / r**2) * eta_short_swift
+
+# Long-range term
+phi_long_gadget2 = (1.  / r ) * corr_short_gadget2
+phi_long_swift = (1.  / r ) * corr_short_swift
+phit_long_gadget2 = corr_long_gadget2 / k**2
+phit_long_swift = corr_long_swift / k**2
+
+
+
+
+figure()
+
+# Potential
+subplot(311, xscale="log", yscale="log")
+
+plot(r_rs, phi_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
+plot(r_rs, phi_short_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(r_rs, phi_short_swift, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+
+xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
+ylim(1.1/r_max, 0.9/r_min)
+ylabel("$\\varphi_s(r)$", labelpad=-3)
+
+legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+
+# Correction
+subplot(312, xscale="log", yscale="log")
+plot(r_rs, np.ones(np.size(r)), '--', lw=1.4, color=colors[0])
+plot(r_rs, 1. - corr_short_gadget2, '-', lw=1.4, color=colors[2])
+plot(r_rs, 1. - corr_short_swift, '-', lw=1.4, color=colors[3])
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
+plot([1., 1.], [-1e5, 1e5], 'k-', alpha=0.5, lw=0.5)
+
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
+ylim(3e-3, 1.5)
+#ylabel("$\\chi_s(r)$", labelpad=-3)
+ylabel("$\\varphi_s(r) \\times r$", labelpad=-2)
+
+# 1 - Correction
+subplot(313, xscale="log", yscale="log")
+plot(r_rs, corr_short_gadget2, '-', lw=1.4, color=colors[2])
+plot(r_rs, corr_short_swift, '-', lw=1.4, color=colors[3])
+
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r)), 'k:', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
+
+xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
+ylim(3e-3, 1.5)
+#ylabel("$1 - \\chi_s(r)$", labelpad=-2)
+ylabel("$1 - \\varphi_s(r) \\times r$", labelpad=-2)
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+xlabel("$r / r_s$", labelpad=-3)
+
+savefig("potential_short.pdf")
+
+##################################################################################################
+
+
+# Force
+figure()
+subplot(311, xscale="log", yscale="log")
+
+plot(r_rs, force_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
+plot(r_rs, force_short_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(r_rs, force_short_swift, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+
+xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
+ylim(1.1/r_max**2, 0.9/r_min**2)
+ylabel("$|\\mathbf{f}_s(r)|$", labelpad=-3)
+yticks([1e-4, 1e-2, 1e0, 1e2], ["$10^{-4}$", "$10^{-2}$", "$10^{0}$", "$10^{2}$"])
+
+legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+
+# Correction
+subplot(312, xscale="log", yscale="log")
+plot(r_rs, np.ones(np.size(r)), '--', lw=1.4, color=colors[0])
+plot(r_rs, eta_short_gadget2, '-', lw=1.4, color=colors[2])
+plot(r_rs, eta_short_swift, '-', lw=1.4, color=colors[3])
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
+plot([1., 1.], [-1e5, 1e5], 'k-', alpha=0.5, lw=0.5)
+
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
+ylim(3e-3, 1.5)
+#ylabel("$\\eta_s(r)$", labelpad=-3)
+ylabel("$|\\mathbf{f}_s(r)|\\times r^2$", labelpad=-2)
+
+# 1 - Correction
+subplot(313, xscale="log", yscale="log")
+plot(r_rs, 1. - eta_short_gadget2, '-', lw=1.4, color=colors[2])
+plot(r_rs, 1. - eta_short_swift, '-', lw=1.4, color=colors[3])
+
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r)), 'k:', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
+
+xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
+ylim(3e-3, 1.5)
+#ylabel("$1 - \\eta_s(r)$", labelpad=-2)
+ylabel("$1 - |\\mathbf{f}_s(r)|\\times r^2$", labelpad=-3)
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+xlabel("$r / r_s$", labelpad=-3)
+
+savefig("force_short.pdf")
+
+##################################################################################################
+
+figure()
+subplot(311, xscale="log", yscale="log")
+
+# Potential
+plot(k_rs, phit_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
+plot(k_rs, phit_long_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(k_rs, phit_long_swift, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot(k_rs, -phit_long_swift, ':', lw=1.4, color=colors[3])
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+
+legend(loc="lower left", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+
+xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
+ylim(1.1/r_max**2, 0.9/r_min**2)
+ylabel("$\\tilde{\\varphi_l}(k)$", labelpad=-3)
+yticks([1e-4, 1e-2, 1e0, 1e2], ["$10^{-4}$", "$10^{-2}$", "$10^{0}$", "$10^{2}$"])
+
+subplot(312, xscale="log", yscale="log")
+
+# Potential normalized
+plot(k_rs, phit_newton * k**2, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
+plot(k_rs, phit_long_gadget2 * k**2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(k_rs, phit_long_swift * k**2, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
+
+xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
+ylim(3e-3, 1.5)
+ylabel("$k^2 \\times \\tilde{\\varphi_l}(k)$", labelpad=-3)
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+
+subplot(313, xscale="log", yscale="log")
+
+plot(k_rs, 1. - phit_long_gadget2 * k**2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(k_rs, 1. - phit_long_swift * k**2, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r)), 'k:', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
+
+xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
+ylim(3e-3, 1.5)
+ylabel("$1 - k^2 \\times \\tilde{\\varphi_l}(k)$", labelpad=-3)
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+
+xlabel("$k \\times r_s$", labelpad=0)
+
+savefig("potential_long.pdf")
diff --git a/theory/Multipoles/plot_potential.py b/theory/Multipoles/plot_potential.py
index 620cc5712c546fe1fc7a22dd9b3fb5483f70226e..8761314572cdbda1304cdf882f920651b58be08e 100644
--- a/theory/Multipoles/plot_potential.py
+++ b/theory/Multipoles/plot_potential.py
@@ -141,7 +141,7 @@ plot([epsilon, epsilon], [-10, 10], 'k-', alpha=0.5, lw=0.5)
 plot([epsilon/plummer_equivalent_factor, epsilon/plummer_equivalent_factor], [0, 10], 'k-', alpha=0.5, lw=0.5)
 
 ylim(0, 2.3)
-ylabel("$\\phi(r)$", labelpad=1)
+ylabel("$\\varphi(r)$", labelpad=1)
 #yticks([0., 0.5, 1., 1.5, 2., 2.5], ["$%.1f$"%(0.*epsilon), "$%.1f$"%(0.5*epsilon), "$%.1f$"%(1.*epsilon), "$%.1f$"%(1.5*epsilon), "$%.1f$"%(2.*epsilon)])
 
 xlim(0,r_max_plot)
@@ -163,6 +163,6 @@ xticks([0., 0.5, 1., 1.5, 2., 2.5], ["$%.1f$"%(0./epsilon), "", "$%.1f$"%(1./eps
 xlabel("$r/H$", labelpad=-7)
 
 ylim(0, 0.95)
-ylabel("$|\\overrightarrow{\\nabla}\\phi(r)|$", labelpad=0)
+ylabel("$|\\overrightarrow{\\nabla}\\varphi(r)|$", labelpad=0)
 
 savefig("potential.pdf")
diff --git a/theory/Multipoles/potential_derivatives.tex b/theory/Multipoles/potential_derivatives.tex
index 9010f06e91e8e0df0d36b539dee09d7b5a67bede..5c7b1e6566d7d51b5d27ea3c24d785571e1ad692 100644
--- a/theory/Multipoles/potential_derivatives.tex
+++ b/theory/Multipoles/potential_derivatives.tex
@@ -7,7 +7,7 @@ the notation $\mathbf{r}=(r_x, r_y, r_z)$, $r = |\mathbf{r}|$ and
 $u=r/H$. Starting from the potential (Eq. \ref{eq:fmm:potential},
 reproduced here for clarity), 
 \begin{align}
-\mathsf{D}_{000}(\mathbf{r}) = \phi (\mathbf{r},H) = 
+\mathsf{D}_{000}(\mathbf{r}) = \varphi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
 \frac{1}{H} \left(-3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3\right) & \mbox{if} & u < 1,\\
 \frac{1}{r} & \mbox{if} & u \geq 1, 
@@ -19,7 +19,7 @@ we can construct the higher order terms by successively applying the
 relevant ones here split by order.
 
 \begin{align}
-\mathsf{D}_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \phi (\mathbf{r},H) = 
+\mathsf{D}_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \varphi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
 -\frac{r_x}{H^3} \left(21u^5 - 90u^4 + 140u^3 - 84u^2 + 14\right) & \mbox{if} & u < 1,\\
 -\frac{r_x}{r^3} & \mbox{if} & u \geq 1, 
@@ -30,7 +30,7 @@ relevant ones here split by order.
 \noindent\rule{6cm}{0.4pt}
 
 \begin{align}
-\mathsf{D}_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \phi (\mathbf{r},H) = 
+\mathsf{D}_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \varphi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
 \frac{r_x^2}{H^5}\left(-105u^3+360u^2-420u+168\right) -
 \frac{1}{H^3} \left(21u^5 - 90u^4 + 140u^3 - 84u^2 + 14\right) & \mbox{if} & u < 1,\\
@@ -40,7 +40,7 @@ relevant ones here split by order.
 \end{align}
 
 \begin{align}
-\mathsf{D}_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \phi (\mathbf{r},H) = 
+\mathsf{D}_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \varphi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
 \frac{r_xr_y}{H^5}\left(-105u^3+360u^2-420u+168\right) & \mbox{if} & u < 1,\\
 3\frac{r_xr_y}{r^5} & \mbox{if} & u \geq 1, 
@@ -51,7 +51,7 @@ relevant ones here split by order.
 \noindent\rule{6cm}{0.4pt}
 
 \begin{align}
-\mathsf{D}_{300}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \phi (\mathbf{r},H) = 
+\mathsf{D}_{300}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \varphi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
 -\frac{r_x^3}{H^7} \left(315u - 720 + 420u^{-1}\right) +
 \frac{3r_x}{H^5}\left(-105u^3+360u^2-420u+168\right) & \mbox{if} & u < 1,\\
@@ -61,7 +61,7 @@ relevant ones here split by order.
 \end{align}
 
 \begin{align}
-\mathsf{D}_{210}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \phi (\mathbf{r},H) = 
+\mathsf{D}_{210}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \varphi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
 -\frac{r_x^2r_y}{H^7} \left(315u - 720 + 420u^{-1}\right) +
 \frac{r_y}{H^5}\left(-105u^3+360u^2-420u+168\right) & \mbox{if} & u < 1,\\
@@ -72,7 +72,7 @@ relevant ones here split by order.
 
 
 \begin{align}
-\mathsf{D}_{111}(\mathbf{r}) = \frac{\partial^3}{\partial r_x\partial r_y\partial r_z} \phi (\mathbf{r},H) = 
+\mathsf{D}_{111}(\mathbf{r}) = \frac{\partial^3}{\partial r_x\partial r_y\partial r_z} \varphi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
 -\frac{r_xr_yr_z}{H^7} \left(315u - 720 + 420u^{-1}\right) & \mbox{if} & u < 1,\\
 -15\frac{r_xr_yr_z}{r^7} & \mbox{if} & u \geq 1, 
diff --git a/theory/Multipoles/potential_softening.tex b/theory/Multipoles/potential_softening.tex
index bdfc27da290d23a8c5fbd1423c442f3948616001..aa9ee12340a3492a19dcf9048548952ef7e141e1 100644
--- a/theory/Multipoles/potential_softening.tex
+++ b/theory/Multipoles/potential_softening.tex
@@ -10,9 +10,10 @@ $\epsilon$. Instead of the commonly used spline kernel of
 is cheaper to compute and has a very similar overall shape. The C2
 kernel has the advantage of being branch-free leading to an expression
 which is faster to evaluate using vector units available on modern
-architectures; it also does not require any division. We set
-$\tilde\delta(\mathbf{x}) = \rho(|\mathbf{x}|) = W(|\mathbf{x}|,
-3\epsilon_{\rm Plummer})$, with $W(r, H)$ given by
+architectures; it also does not require any divisions to evaluate the
+softened forces. We set $\tilde\delta(\mathbf{x}) =
+\rho(|\mathbf{x}|) = W(|\mathbf{x}|, 3\epsilon_{\rm Plummer})$, with
+$W(r, H)$ given by
 
 \begin{align}
 W(r,H) &= \frac{21}{2\pi H^3} \times \nonumber \\
@@ -22,9 +23,9 @@ W(r,H) &= \frac{21}{2\pi H^3} \times \nonumber \\
 \end{array}
 \right.
 \end{align}
-and $u = r/H$. The potential $\phi(r,H)$ corresponding to this density distribution reads
+and $u = r/H$. The potential $\varphi(r,H)$ corresponding to this density distribution reads
 \begin{align}
-\phi = 
+\varphi = 
 \left\lbrace\begin{array}{rcl}
 \frac{1}{H} (-3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3) & \mbox{if} & u < 1,\\
 \frac{1}{r} & \mbox{if} & u \geq 1.
diff --git a/theory/Multipoles/run.sh b/theory/Multipoles/run.sh
index 3b2156a40c4309a1fb21980fb14022f963ce4f51..fc376188ad2e69d2879ce963ddc7069c736fc8b7 100755
--- a/theory/Multipoles/run.sh
+++ b/theory/Multipoles/run.sh
@@ -1,9 +1,14 @@
 #!/bin/bash
-if [! -e potential.pdf ]
+if [ ! -e potential.pdf ]
 then
-    echo "Generating figures..."
+    echo "Generating 1st figure..."
     python plot_potential.py
 fi
+if [ ! -e potential_short.pdf ]
+then
+    echo "Generating 2nd figures..."
+    python plot_mesh.py
+fi
 echo "Generating PDF..."
 pdflatex -jobname=fmm fmm_standalone.tex
 bibtex fmm.aux