diff --git a/.gitignore b/.gitignore
index f554fa05e392948ce823b0be52d355e5c2e8e333..29e6a7672653a7a8d752de4e8035da95cd30b526 100644
--- a/.gitignore
+++ b/.gitignore
@@ -105,7 +105,7 @@ tests/testDump
tests/testLogger
tests/benchmarkInteractions
tests/testGravityDerivatives
-tests/testPotential
+tests/testPotentialSelf
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b80fc3a5608f50eb4b6fc68f95f74112be152f68..3a4f451a4a3f69907772ce63d0155fa18ba60444 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -26,7 +26,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
- testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotential
+ testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -36,7 +36,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
- testGravityDerivatives testPotential
+ testGravityDerivatives testPotentialSelf
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -100,7 +100,7 @@ testLogger_SOURCES = testLogger.c
testGravityDerivatives_SOURCES = testGravityDerivatives.c
-testPotential_SOURCES = testPotential.c
+testPotentialSelf_SOURCES = testPotentialSelf.c
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c
new file mode 100644
index 0000000000000000000000000000000000000000..285cf11ecf8d8a44bda87b2e1196aa3baac605cf
--- /dev/null
+++ b/tests/testPotentialSelf.c
@@ -0,0 +1,186 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 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/>.
+ *
+ ******************************************************************************/
+#include "../config.h"
+
+/* Some standard headers. */
+#include <fenv.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Local headers. */
+#include "runner_doiact_grav.h"
+#include "swift.h"
+
+const int num_tests = 100;
+const double eps = 0.02;
+
+/**
+ * @brief Check that a and b are consistent (up to some relative error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param s String used to identify this check in messages
+ */
+void check_value(double a, double b, const char *s) {
+ if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6)
+ error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
+}
+
+/* Definitions of the potential and force that match
+ exactly the theory document */
+double S(double x) { return exp(x) / (1. + exp(x)); }
+
+double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); }
+
+double potential(double mass, double r, double H, double rlr) {
+
+ const double u = r / H;
+ const double x = r / rlr;
+ double pot;
+ if (u > 1.)
+ pot = -mass / r;
+ else
+ pot = -mass *
+ (-3. * u * u * u * u * u * u * u + 15. * u * u * u * u * u * u -
+ 28. * u * u * u * u * u + 21. * u * u * u * u - 7. * u * u + 3.) /
+ H;
+
+ return pot * (2. - 2. * S(2. * x));
+}
+
+double acceleration(double mass, double r, double H, double rlr) {
+
+ const double u = r / H;
+ const double x = r / rlr;
+ double acc;
+ if (u > 1.)
+ acc = -mass / (r * r * r);
+ else
+ acc = -mass * (21. * u * u * u * u * u - 90. * u * u * u * u +
+ 140. * u * u * u - 84. * u * u + 14.) /
+ (H * H * H);
+
+ return r * acc * (4. * x * S_prime(2 * x) - 2. * S(2. * x) + 2.);
+}
+
+int main() {
+
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+ /* Initialise a few things to get us going */
+ struct engine e;
+ e.max_active_bin = num_time_bins;
+ e.time = 0.1f;
+ e.ti_current = 8;
+ e.time_base = 1e-10;
+
+ struct space s;
+ s.width[0] = 0.2;
+ s.width[1] = 0.2;
+ s.width[2] = 0.2;
+ e.s = &s;
+
+ struct gravity_props props;
+ props.a_smooth = 1.25;
+ e.gravity_properties = &props;
+
+ struct runner r;
+ bzero(&r, sizeof(struct runner));
+ r.e = &e;
+
+ const double rlr = props.a_smooth * s.width[0];
+
+ /* Init the cache for gravity interaction */
+ gravity_cache_init(&r.ci_gravity_cache, num_tests * 2);
+
+ /* Let's create one cell with a massive particle and a bunch of test particles
+ */
+ struct cell c;
+ bzero(&c, sizeof(struct cell));
+ c.width[0] = 1.;
+ c.width[1] = 1.;
+ c.width[2] = 1.;
+ c.loc[0] = 0.;
+ c.loc[1] = 0.;
+ c.loc[2] = 0.;
+ c.gcount = 1 + num_tests;
+ c.ti_old_gpart = 8;
+ c.ti_gravity_end_min = 8;
+ c.ti_gravity_end_max = 8;
+
+ posix_memalign((void **)&c.gparts, gpart_align,
+ c.gcount * sizeof(struct gpart));
+ bzero(c.gparts, c.gcount * sizeof(struct gpart));
+
+ /* Create the massive particle */
+ c.gparts[0].x[0] = 0.;
+ c.gparts[0].x[1] = 0.5;
+ c.gparts[0].x[2] = 0.5;
+ c.gparts[0].mass = 1.;
+ c.gparts[0].epsilon = eps;
+ c.gparts[0].time_bin = 1;
+ c.gparts[0].type = swift_type_dark_matter;
+ c.gparts[0].id_or_neg_offset = 1;
+#ifdef SWIFT_DEBUG_CHECKS
+ c.gparts[0].ti_drift = 8;
+#endif
+
+ /* Create the mass-less particles */
+ for (int n = 1; n < num_tests + 1; ++n) {
+
+ struct gpart *gp = &c.gparts[n];
+
+ gp->x[0] = n / ((double)num_tests);
+ gp->x[1] = 0.5;
+ gp->x[2] = 0.5;
+ gp->mass = 0.;
+ gp->epsilon = eps;
+ gp->time_bin = 1;
+ gp->type = swift_type_dark_matter;
+ gp->id_or_neg_offset = n + 1;
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->ti_drift = 8;
+#endif
+ }
+
+ /* Now compute the forces */
+ runner_doself_grav_pp_truncated(&r, &c);
+
+ /* Verify everything */
+ for (int n = 1; n < num_tests + 1; ++n) {
+ const struct gpart *gp = &c.gparts[n];
+
+ double pot_true = potential(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+ double acc_true =
+ acceleration(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+
+ // message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true);
+
+ check_value(gp->potential, pot_true, "potential");
+ check_value(gp->a_grav[0], acc_true, "acceleration");
+ }
+
+ free(c.gparts);
+ return 0;
+}