diff --git a/.gitignore b/.gitignore
index 3a8146ae84db2d4f2ec5be44853f1d861b5ce5c4..d69de3deae4cbe1dff23fd0f22704f6b4f5e40ed 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,13 +18,17 @@ doc/html/
 doc/latex/
 doc/man/
 doc/Doxyfile
-examples/test
-examples/test_fixdt
-examples/test_fixdt_mpi
-examples/test_mindt
-examples/test_mindt_mpi
-examples/test_mpi
-examples/test_single
+examples/swift
+examples/swift_fixdt
+examples/swift_fixdt_mpi
+examples/swift_mindt
+examples/swift_mindt_mpi
+examples/swift_mpi
+
+tests/testGreetings
+tests/testReading
+tests/input.hdf5
+tests/testSingle
 
 m4/libtool.m4
 m4/ltoptions.m4
diff --git a/Makefile.am b/Makefile.am
index 946534d512007f514eb0c67c2489f22b5bb362c7..9d4f6371b7aa37a1750c5fb8dbed17f9ff48442e 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,6 +1,6 @@
 
 # This file is part of SWIFT.
-# Coypright (c) 2012 pedro.gonnet@durham.ac.uk.
+# Copyright (c) 2012 pedro.gonnet@durham.ac.uk.
 #               2015 matthieu.schaller@durham.ac.uk.
 # 
 # This program is free software: you can redistribute it and/or modify
diff --git a/examples/BigCosmoVolume/makeIC.py b/examples/BigCosmoVolume/makeIC.py
index d5d5392f07671cfc529b76bcb3d00c3f50d3cf52..c4699b5616a634456c245cc33c01987580da4728 100644
--- a/examples/BigCosmoVolume/makeIC.py
+++ b/examples/BigCosmoVolume/makeIC.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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
diff --git a/examples/BigPerturbedBox/makeIC_fcc.py b/examples/BigPerturbedBox/makeIC_fcc.py
index b21cb61f9d6ef9678eeddcc45e07908fb79a4fb0..92bfd9ce24e0018c33c2f68455869ced35b845fa 100644
--- a/examples/BigPerturbedBox/makeIC_fcc.py
+++ b/examples/BigPerturbedBox/makeIC_fcc.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/CosmoVolume/getIC.sh b/examples/CosmoVolume/getIC.sh
index ebee5e63d53d17019d3d8ba2c58e7bed26ff2641..5b951fc247006f0378bcee2b74a9619b7b1acc7d 100755
--- a/examples/CosmoVolume/getIC.sh
+++ b/examples/CosmoVolume/getIC.sh
@@ -1,2 +1,2 @@
 #!/bin/bash
-wget http://www.dur.ac.uk/pedro.gonnet/cosmoVolume.hdf5
+wget http://icc.dur.ac.uk/~jlvc76/Files/SWIFT/cosmoVolume.hdf5
diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex/makeIC.py
index 324dcb65b21527f4e3868571d08fe7a0c9ff98a1..533199869ee325dd75831410f3e9f1181ae30ebd 100644
--- a/examples/GreshoVortex/makeIC.py
+++ b/examples/GreshoVortex/makeIC.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/GreshoVortex/solution.py b/examples/GreshoVortex/solution.py
index a3937d1d144975c8f0023f1b310ace06cbf6f6b2..5b282caf2d7f3311ddb595781fed74c51bb4819f 100644
--- a/examples/GreshoVortex/solution.py
+++ b/examples/GreshoVortex/solution.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 9d4fd5951374d249e76980689d629079749b2c8f..38ba3e7dbf5f904709683ae9ba883c83831f2000 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -1,6 +1,6 @@
 
 # This file is part of SWIFT.
-# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
 #
 # This program is free software: you can redistribute it and/or modify
@@ -29,45 +29,40 @@ MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
 MPI_FLAGS = -DWITH_MPI $(METIS_INCS)
 
 # Set-up the library
-bin_PROGRAMS = test test_fixdt test_mindt test_single
+bin_PROGRAMS = swift swift_fixdt swift_mindt
 
 # Build MPI versions as well?
 if HAVEMPI
-bin_PROGRAMS += test_mpi test_fixdt_mpi test_mindt_mpi
+bin_PROGRAMS += swift_mpi swift_fixdt_mpi swift_mindt_mpi
 endif
 
-# Sources for test
-test_SOURCES = test.c
-test_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep | engine_policy_setaffinity"
-test_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
+# Sources for swift
+swift_SOURCES = main.c
+swift_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep | engine_policy_setaffinity"
+swift_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
 
-# Sources for test_fixdt
-test_fixdt_SOURCES = test.c
-test_fixdt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep | engine_policy_setaffinity"
-test_fixdt_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
+# Sources for swift_fixdt
+swift_fixdt_SOURCES = main.c
+swift_fixdt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep | engine_policy_setaffinity"
+swift_fixdt_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
 
-# Sources for test_mindt
-test_mindt_SOURCES = test.c
-test_mindt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_keep | engine_policy_setaffinity"
-test_mindt_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
+# Sources for swift_mindt
+swift_mindt_SOURCES = main.c
+swift_mindt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_keep | engine_policy_setaffinity"
+swift_mindt_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
 
-# Sources for test_mpi
-test_mpi_SOURCES = test.c
-test_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep"
-test_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
+# Sources for swift_mpi
+swift_mpi_SOURCES = main.c
+swift_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep"
+swift_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
 
-# Sources for test_fixdt_mpi
-test_fixdt_mpi_SOURCES = test.c
-test_fixdt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep"
-test_fixdt_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
+# Sources for swift_fixdt_mpi
+swift_fixdt_mpi_SOURCES = main.c
+swift_fixdt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep"
+swift_fixdt_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
 
-# Sources for test_mindt_mpi
-test_mindt_mpi_SOURCES = test.c
-test_mindt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_keep"
-test_mindt_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
-
-# Sources for test_single
-test_single_SOURCES = test_single.c
-test_single_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
-test_single_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
+# Sources for swift_mindt_mpi
+swift_mindt_mpi_SOURCES = main.c
+swift_mindt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_keep"
+swift_mindt_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
 
diff --git a/examples/PertubedBox/makeIC.py b/examples/PertubedBox/makeIC.py
index fe32c6f10a0f836cfc351c87730fd6a32882bad8..81ae9e5909ec51cb640209ce759d618311ce811f 100644
--- a/examples/PertubedBox/makeIC.py
+++ b/examples/PertubedBox/makeIC.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/SedovBlast/makeIC.py b/examples/SedovBlast/makeIC.py
index abddc803bc6bcb92a0c44eb5ae4e2977a0281aca..9b4b5443f472edf5bb299ed5e7261d115c96293c 100644
--- a/examples/SedovBlast/makeIC.py
+++ b/examples/SedovBlast/makeIC.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/SedovBlast/makeIC_fcc.py b/examples/SedovBlast/makeIC_fcc.py
index 7e684504095e76c39e92f9b2f954b71aa573fb22..88cbaf8042323ea91ed7dd09b1bd63418aff3e3f 100644
--- a/examples/SedovBlast/makeIC_fcc.py
+++ b/examples/SedovBlast/makeIC_fcc.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/SedovBlast/rdf.py b/examples/SedovBlast/rdf.py
index 52b834f3947d1f992093dad9916a6164d5b3b0fa..7f932cc814dc36e14e1bef52e33cf5ed1f527dfd 100644
--- a/examples/SedovBlast/rdf.py
+++ b/examples/SedovBlast/rdf.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
@@ -29,7 +29,7 @@ from numpy import *
 
 # Input values?
 if len(sys.argv) < 3 :
-    print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
+    print "Usage: " , sys.argv[0] , " <filename> <nr. bins>"
     exit()
     
 # Get the input arguments
diff --git a/examples/SedovBlast/solution.py b/examples/SedovBlast/solution.py
index 947c9dde5a02cea1f0edcbc60b67f0e3490625c2..9335e22bdd76585e8d53d3dba9f9a435197f92fc 100644
--- a/examples/SedovBlast/solution.py
+++ b/examples/SedovBlast/solution.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/SodShock/makeIC.py b/examples/SodShock/makeIC.py
index bfd1ede773d0cec42c8b5d0e3a75223fd3abdb5d..8405b192697580e5d9bc586bba1c6fe1265552cf 100644
--- a/examples/SodShock/makeIC.py
+++ b/examples/SodShock/makeIC.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/SodShock/rhox.py b/examples/SodShock/rhox.py
index 36771c5c0d9dcdec99403249ea8a7bb1a4bc58c6..70493be3728cdeb27409a79f616fa3ec5bb9cdfd 100644
--- a/examples/SodShock/rhox.py
+++ b/examples/SodShock/rhox.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
@@ -29,7 +29,7 @@ from numpy import *
 
 # Input values?
 if len(sys.argv) < 3 :
-    print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
+    print "Usage: " , sys.argv[0] , " <filename> <nr. bins>"
     exit()
     
 # Get the input arguments
diff --git a/examples/SodShock/solution.py b/examples/SodShock/solution.py
index a5a1364d90d4af12f3ac051e23748d6228592258..39f25c625232eee9bae0300339955f775f3b46ed 100644
--- a/examples/SodShock/solution.py
+++ b/examples/SodShock/solution.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox/makeIC.py
index 2600d79a4a70de3e69fb76b82e11f6a4b72a73d2..9f740398164a3411453f5a55a011428ececaeb19 100644
--- a/examples/UniformBox/makeIC.py
+++ b/examples/UniformBox/makeIC.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
diff --git a/examples/test.c b/examples/main.c
similarity index 53%
rename from examples/test.c
rename to examples/main.c
index d9986a710a4959dac5e3f7a4d63d066aaf8807f2..2c1087db4a1508828cd076de95e5b79811058b5e 100644
--- a/examples/test.c
+++ b/examples/main.c
@@ -56,493 +56,8 @@
 #define ENGINE_POLICY engine_policy_none
 #endif
 
-/**
- * @brief Mapping function to draw a specific cell (gnuplot).
- */
-
-void map_cells_plot(struct cell *c, void *data) {
-
-  int depth = *(int *)data;
-  double *l = c->loc, *h = c->h;
-
-  if (c->depth <= depth) {
-
-    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2]);
-    printf("%.16e %.16e %.16e\n\n\n", l[0], l[1] + h[1], l[2]);
-
-    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n\n\n", l[0], l[1] + h[1], l[2] + h[2]);
-
-    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n\n\n", l[0], l[1], l[2] + h[2]);
-
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1], l[2] + h[2]);
-
-    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1], l[2]);
-
-    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2]);
-    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
-    printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1] + h[1], l[2]);
-
-    if (!c->split) {
-      for (int k = 0; k < c->count; k++)
-        printf("0 0 0 %.16e %.16e %.16e\n", c->parts[k].x[0], c->parts[k].x[1],
-               c->parts[k].x[2]);
-      printf("\n\n");
-    }
-    /* else
-        for ( int k = 0 ; k < 8 ; k++ )
-            if ( c->progeny[k] != NULL )
-                map_cells_plot( c->progeny[k] , data ); */
-  }
-}
-
-/**
- * @brief Mapping function for checking if each part is in its box.
- */
-
-/* void map_check ( struct part *p , struct cell *c , void *data ) {
-
-    if ( p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
-         p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
-         p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] )
-        printf( "map_check: particle %i is outside of its box.\n" , p->id );
-
-    } */
-
-/**
- * @brief Mapping function for neighbour count.
- */
-
-void map_cellcheck(struct cell *c, void *data) {
-
-  int k, *count = (int *)data;
-  struct part *p;
-
-  __sync_fetch_and_add(count, c->count);
-
-  /* Loop over all parts and check if they are in the cell. */
-  for (k = 0; k < c->count; k++) {
-    p = &c->parts[k];
-    if (p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
-        p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] ||
-        p->x[2] > c->loc[2] + c->h[2]) {
-      printf(
-          "map_cellcheck: particle at [ %.16e %.16e %.16e ] outside of cell [ "
-          "%.16e %.16e %.16e ] - [ %.16e %.16e %.16e ].\n",
-          p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
-          c->loc[0] + c->h[0], c->loc[1] + c->h[1], c->loc[2] + c->h[2]);
-      error("particle out of bounds!");
-    }
-  }
-}
-
-/**
- * @brief Mapping function for maxdepth cell count.
- */
-
-void map_maxdepth(struct cell *c, void *data) {
-
-  int maxdepth = ((int *)data)[0];
-  int *count = &((int *)data)[1];
-
-  // printf( "%e\n" , p->count );
-
-  if (c->depth == maxdepth) *count += 1;
-}
-
-/**
- * @brief Mapping function for neighbour count.
- */
-
-void map_count(struct part *p, struct cell *c, void *data) {
-
-  double *wcount = (double *)data;
-
-  // printf( "%i %e %e\n" , p->id , p->count , p->count_dh );
-
-  *wcount += p->density.wcount;
-}
-
-void map_wcount_min(struct part *p, struct cell *c, void *data) {
-
-  struct part **p2 = (struct part **)data;
-
-  if (p->density.wcount < (*p2)->density.wcount) *p2 = p;
-}
-
-void map_wcount_max(struct part *p, struct cell *c, void *data) {
-
-  struct part **p2 = (struct part **)data;
-
-  if (p->density.wcount > (*p2)->density.wcount) *p2 = p;
-}
-
-void map_h_min(struct part *p, struct cell *c, void *data) {
-
-  struct part **p2 = (struct part **)data;
-
-  if (p->h < (*p2)->h) *p2 = p;
-}
-
-void map_h_max(struct part *p, struct cell *c, void *data) {
-
-  struct part **p2 = (struct part **)data;
-
-  if (p->h > (*p2)->h) *p2 = p;
-}
-
-/**
- * @brief Mapping function for neighbour count.
- */
-
-void map_icount(struct part *p, struct cell *c, void *data) {
-
-  // int *count = (int *)data;
-
-  // printf( "%i\n" , p->icount );
-
-  // *count += p->icount;
-}
-
-/**
- * @brief Mapping function to print the particle position.
- */
-
-void map_dump(struct part *p, struct cell *c, void *data) {
-
-  double *shift = (double *)data;
-
-  printf("%g\t%g\t%g\n", p->x[0] - shift[0], p->x[1] - shift[1],
-         p->x[2] - shift[2]);
-}
-
-/**
- * @brief Compute the average number of pairs per particle using
- *      a brute-force O(N^2) computation.
- *
- * @param dim The space dimensions.
- * @param parts The #part array.
- * @param N The number of parts.
- * @param periodic Periodic boundary conditions flag.
- */
-
-void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
-              int periodic) {
-
-  int i, j, k, count = 0;
-  // int mj, mk;
-  // double maxratio = 1.0;
-  double r2, dx[3], rho = 0.0;
-  double rho_max = 0.0, rho_min = 100;
-
-  /* Loop over all particle pairs. */
-  for (j = 0; j < N; j++) {
-    if (j % 1000 == 0) {
-      printf("pairs_n2: j=%i.\n", j);
-      fflush(stdout);
-    }
-    for (k = j + 1; k < N; k++) {
-      for (i = 0; i < 3; i++) {
-        dx[i] = parts[j].x[i] - parts[k].x[i];
-        if (periodic) {
-          if (dx[i] < -dim[i] / 2)
-            dx[i] += dim[i];
-          else if (dx[i] > dim[i] / 2)
-            dx[i] -= dim[i];
-        }
-      }
-      r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-      if (r2 < parts[j].h * parts[j].h || r2 < parts[k].h * parts[k].h) {
-        runner_iact_density(r2, NULL, parts[j].h, parts[k].h, &parts[j],
-                            &parts[k]);
-        /* if ( parts[j].h / parts[k].h > maxratio )
-            {
-            maxratio = parts[j].h / parts[k].h;
-            mj = j; mk = k;
-            }
-        else if ( parts[k].h / parts[j].h > maxratio )
-            {
-            maxratio = parts[k].h / parts[j].h;
-            mj = j; mk = k;
-            } */
-      }
-    }
-  }
 
-  /* Aggregate the results. */
-  for (k = 0; k < N; k++) {
-    // count += parts[k].icount;
-    rho += parts[k].density.wcount;
-    rho_min = fmin(parts[k].density.wcount, rho_min);
-    rho_min = fmax(parts[k].density.wcount, rho_max);
-  }
 
-  /* Dump the result. */
-  printf("pairs_n2: avg. density per part is %.3f (nr. pairs %.3f).\n",
-         rho / N + 32.0 / 3, ((double)count) / N);
-  printf("pairs_n2: densities are in [ %e , %e ].\n", rho_min / N + 32.0 / 3,
-         rho_max / N + 32.0 / 3);
-  /* printf( "pairs_n2: maximum ratio between parts %i [%e,%e,%e] and %i
-     [%e,%e,%e] is %.3f/%.3f\n" ,
-      mj , parts[mj].x[0] , parts[mj].x[1] , parts[mj].x[2] ,
-      mk , parts[mk].x[0] , parts[mk].x[1] , parts[mk].x[2] ,
-      parts[mj].h , parts[mk].h ); fflush(stdout); */
-  fflush(stdout);
-}
-
-void pairs_single_density(double *dim, long long int pid,
-                          struct part *__restrict__ parts, int N,
-                          int periodic) {
-
-  int i, k;
-  // int mj, mk;
-  // double maxratio = 1.0;
-  double r2, dx[3];
-  float fdx[3];
-  struct part p;
-  // double ih = 12.0/6.25;
-
-  /* Find "our" part. */
-  for (k = 0; k < N && parts[k].id != pid; k++)
-    ;
-  if (k == N) error("Part not found.");
-  p = parts[k];
-  printf("pairs_single: part[%i].id == %lli.\n", k, pid);
-
-  p.rho = 0.0;
-  p.density.wcount = 0.0;
-  // p.icount = 0;
-  p.rho_dh = 0.0;
-
-  /* Loop over all particle pairs. */
-  for (k = 0; k < N; k++) {
-    if (parts[k].id == p.id) continue;
-    for (i = 0; i < 3; i++) {
-      dx[i] = p.x[i] - parts[k].x[i];
-      if (periodic) {
-        if (dx[i] < -dim[i] / 2)
-          dx[i] += dim[i];
-        else if (dx[i] > dim[i] / 2)
-          dx[i] -= dim[i];
-      }
-      fdx[i] = dx[i];
-    }
-    r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
-    if (r2 < p.h * p.h) {
-      runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
-      /* printf( "pairs_simple: interacting particles %lli [%i,%i,%i] and %lli
-         [%i,%i,%i], r=%e.\n" ,
-          pid , (int)(p.x[0]*ih) , (int)(p.x[1]*ih) , (int)(p.x[2]*ih) ,
-          parts[k].id , (int)(parts[k].x[0]*ih) , (int)(parts[k].x[1]*ih) ,
-         (int)(parts[k].x[2]*ih) ,
-          sqrtf(r2) ); */
-    }
-  }
-
-  /* Dump the result. */
-  printf("pairs_single: wcount of part %lli (h=%e) is %f.\n", p.id, p.h,
-         p.density.wcount + 32.0 / 3);
-  fflush(stdout);
-}
-
-void pairs_single_grav(double *dim, long long int pid,
-                       struct gpart *__restrict__ parts, int N, int periodic) {
-
-  int i, k;
-  // int mj, mk;
-  // double maxratio = 1.0;
-  double r2, dx[3];
-  float fdx[3], a[3] = {0.0, 0.0, 0.0}, aabs[3] = {0.0, 0.0, 0.0};
-  struct gpart pi, pj;
-  // double ih = 12.0/6.25;
-
-  /* Find "our" part. */
-  for (k = 0; k < N; k++)
-    if ((parts[k].id > 0 && parts[k].part->id == pid) || parts[k].id == -pid)
-      break;
-  if (k == N) error("Part not found.");
-  pi = parts[k];
-  pi.a[0] = 0.0f;
-  pi.a[1] = 0.0f;
-  pi.a[2] = 0.0f;
-
-  /* Loop over all particle pairs. */
-  for (k = 0; k < N; k++) {
-    if (parts[k].id == pi.id) continue;
-    pj = parts[k];
-    for (i = 0; i < 3; i++) {
-      dx[i] = pi.x[i] - pj.x[i];
-      if (periodic) {
-        if (dx[i] < -dim[i] / 2)
-          dx[i] += dim[i];
-        else if (dx[i] > dim[i] / 2)
-          dx[i] -= dim[i];
-      }
-      fdx[i] = dx[i];
-    }
-    r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
-    runner_iact_grav(r2, fdx, &pi, &pj);
-    a[0] += pi.a[0];
-    a[1] += pi.a[1];
-    a[2] += pi.a[2];
-    aabs[0] += fabsf(pi.a[0]);
-    aabs[1] += fabsf(pi.a[1]);
-    aabs[2] += fabsf(pi.a[2]);
-    pi.a[0] = 0.0f;
-    pi.a[1] = 0.0f;
-    pi.a[2] = 0.0f;
-  }
-
-  /* Dump the result. */
-  message(
-      "acceleration on gpart %lli is a=[ %e %e %e ], |a|=[ %.2e %.2e %.2e ].\n",
-      pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
-}
-
-/**
- * @brief Test the kernel function by dumping it in the interval [0,1].
- *
- * @param N number of intervals in [0,1].
- */
-
-void kernel_dump(int N) {
-
-  int k;
-  float x, w, dw_dx;
-  float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
-  float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
-  // float dw_dx4[4] __attribute__ ((aligned (16)));
-
-  for (k = 0; k <= N; k++) {
-    x = ((float)k) / N;
-    x4[3] = x4[2];
-    x4[2] = x4[1];
-    x4[1] = x4[0];
-    x4[0] = x;
-    kernel_deval(x, &w, &dw_dx);
-    // kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
-    printf(" %e %e %e %e %e %e %e\n", x, w, dw_dx, w4[0], w4[1], w4[2], w4[3]);
-  }
-}
-
-float gadget(float r) {
-  float fac, h_inv, u, r2 = r * r;
-  if (r >= const_epsilon)
-    fac = 1.0f / (r2 * r);
-  else {
-    h_inv = 1. / const_epsilon;
-    u = r * h_inv;
-    if (u < 0.5)
-      fac = const_iepsilon3 * (10.666666666667 + u * u * (32.0 * u - 38.4));
-    else
-      fac = const_iepsilon3 *
-            (21.333333333333 - 48.0 * u + 38.4 * u * u -
-             10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
-  }
-  return const_G * fac;
-}
-
-void gravity_dump(float r_max, int N) {
-
-  int k;
-  float x, w;
-  float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
-  float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
-  // float dw_dx4[4] __attribute__ ((aligned (16)));
-
-  for (k = 1; k <= N; k++) {
-    x = (r_max * k) / N;
-    x4[3] = x4[2];
-    x4[2] = x4[1];
-    x4[1] = x4[0];
-    x4[0] = x;
-    kernel_grav_eval(x, &w);
-    w *= const_G / (x * x * x);
-    // blender_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
-    printf(" %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n", x, w * x, w4[0],
-           w4[1], w4[2], w4[3], gadget(x) * x);
-  }
-}
-
-/**
- * @brief Test the density function by dumping it for two random parts.
- *
- * @param N number of intervals in [0,1].
- */
-
-void density_dump(int N) {
-
-  int k;
-  float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
-  struct part *pi[4], *pj[4], Pi[4], Pj[4];
-
-  /* Init the interaction parameters. */
-  for (k = 0; k < 4; k++) {
-    Pi[k].mass = 1.0f;
-    Pi[k].rho = 0.0f;
-    Pi[k].density.wcount = 0.0f;
-    Pj[k].mass = 1.0f;
-    Pj[k].rho = 0.0f;
-    Pj[k].density.wcount = 0.0f;
-    hi[k] = 1.0;
-    hj[k] = 1.0;
-    pi[k] = &Pi[k];
-    pj[k] = &Pj[k];
-  }
-
-  for (k = 0; k <= N; k++) {
-    r2[3] = r2[2];
-    r2[2] = r2[1];
-    r2[1] = r2[0];
-    r2[0] = ((float)k) / N;
-    Pi[0].density.wcount = 0;
-    Pj[0].density.wcount = 0;
-    runner_iact_density(r2[0], NULL, hi[0], hj[0], &Pi[0], &Pj[0]);
-    printf(" %e %e %e", r2[0], Pi[0].density.wcount, Pj[0].density.wcount);
-    Pi[0].density.wcount = 0;
-    Pj[0].density.wcount = 0;
-    Pi[1].density.wcount = 0;
-    Pj[1].density.wcount = 0;
-    Pi[2].density.wcount = 0;
-    Pj[2].density.wcount = 0;
-    Pi[3].density.wcount = 0;
-    Pj[3].density.wcount = 0;
-    runner_iact_vec_density(r2, NULL, hi, hj, pi, pj);
-    printf(" %e %e %e %e\n", Pi[0].density.wcount, Pi[1].density.wcount,
-           Pi[2].density.wcount, Pi[3].density.wcount);
-  }
-}
-
-/**
- *  Factorize a given integer, attempts to keep larger pair of factors.
- */
-void factor(int value, int *f1, int *f2) {
-  int j;
-  int i;
-
-  j = (int)sqrt(value);
-  for (i = j; i > 0; i--) {
-    if ((value % i) == 0) {
-      *f1 = i;
-      *f2 = value / i;
-      break;
-    }
-  }
-}
 
 /**
  * @brief Main routine that loads a few particles and generates some output.
@@ -696,7 +211,7 @@ int main(int argc, char *argv[]) {
             (long int)sizeof(struct gpart));
   }
 
-  /* Initilaize unit system */
+  /* Initialize unit system */
   initUnitSystem(&us);
   if (myrank == 0) {
     message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
@@ -781,7 +296,7 @@ int main(int argc, char *argv[]) {
     // message( "cutoffs in [ %g %g ]." , s.h_min , s.h_max ); fflush(stdout);
   }
 
-  /* Verify that each particle is in it's propper cell. */
+  /* Verify that each particle is in it's proper cell. */
   if (myrank == 0) {
     icount = 0;
     space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
@@ -890,7 +405,7 @@ int main(int argc, char *argv[]) {
 #endif
     }
 
-    /* Dump a line of agregate output. */
+    /* Dump a line of aggregate output. */
     /*     if (myrank == 0) { */
     /*       printf("%i %e %.16e %.16e %.16e %.3e %.3e %i %.3e %.3e", j, e.time,
      */
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index 3e01be7da31ff87ba137cae394a203d5a6de7b25..52f1a90592120e61bd3e789ee0cb8fe7edb63a8d 100644
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  #                    Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
diff --git a/examples/runs.sh b/examples/runs.sh
index ff8e3339d69bd414c9c55b4ed60d532db121d55c..e31dc5a090eb512f9e5ae3b301610cd51839d9bb 100755
--- a/examples/runs.sh
+++ b/examples/runs.sh
@@ -34,31 +34,31 @@ do
     # Sod-Shock runs
     if [ ! -e SodShock_mindt_${cpu}.dump ]
     then
-        ./test_mindt -c 1.0 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1.0 > SodShock_${cpu}.dump
+        ./swift_mindt -c 1.0 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1.0 > SodShock_${cpu}.dump
     fi
     if [ ! -e SodShock_fixed_${cpu}.dump ]
     then
-        ./test_fixdt -r 1000 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1e-4 > SodShock_fixed_${cpu}.dump
+        ./swift_fixdt -r 1000 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1e-4 > SodShock_fixed_${cpu}.dump
     fi
     
     # Sedov blast
     if [ ! -e SedovBlast_mindt_${cpu}.dump ]
     then
-        ./test_mindt -c 0.2 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 1e-10 > SedovBlast_${cpu}.dump
+        ./swift_mindt -c 0.2 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 1e-10 > SedovBlast_${cpu}.dump
     fi
     if [ ! -e SedovBlast_fixed_${cpu}.dump ]
     then
-        ./test_fixdt -r 4096 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 5e-5 > SedovBlast_fixed_${cpu}.dump
+        ./swift_fixdt -r 4096 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 5e-5 > SedovBlast_fixed_${cpu}.dump
     fi
     
     # Cosmological volume
     if [ ! -e CosmoVolume_mindt_${cpu}.dump ]
     then
-        ./test_mindt -c 0.01 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 0.01 > CosmoVolume_${cpu}.dump
+        ./swift_mindt -c 0.01 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 0.01 > CosmoVolume_${cpu}.dump
     fi
     if [ ! -e CosmoVolume_fixed_${cpu}.dump ]
     then
-        ./test_fixdt -r 256 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 1e-8 > CosmoVolume_fixed_${cpu}.dump
+        ./swift_fixdt -r 256 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 1e-8 > CosmoVolume_fixed_${cpu}.dump
     fi
 
 done
diff --git a/examples/runs_mpi_cv.sh b/examples/runs_mpi_cv.sh
index 4a009bbdbeb18d003b7572cc1dcb7b7e3fa19a2c..bffa2d56a6a8f3624c0919db8a046dd2db0b80b4 100755
--- a/examples/runs_mpi_cv.sh
+++ b/examples/runs_mpi_cv.sh
@@ -23,43 +23,43 @@ for cpu in $(seq 1 $CPN)
 do
     if [ ! -e ${PREFIX}_${QUEUE}_1x${cpu}.dump ]
     then
-        bsub -oo ${PREFIX}_${QUEUE}_1x${cpu}.dump -q ${QUEUE} -P ${PROJECT} -x -n 1 -R "span[ptile=1]" ./test_fixdt -r $RUNS -t $cpu -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+        bsub -oo ${PREFIX}_${QUEUE}_1x${cpu}.dump -q ${QUEUE} -P ${PROJECT} -x -n 1 -R "span[ptile=1]" ./swift_fixdt -r $RUNS -t $cpu -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
     fi
 done
 
 # Multi-node runs
 if [ ! -e ${PREFIX}_${QUEUE}_2x${cpu}.dump ]
 then
-    bsub -oo ${PREFIX}_${QUEUE}_2x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 2 -R "span[ptile=1]" mpirun -np 2 ./test_fixdt_mpi -r $RUNS -t $CPN -g "2 1 1" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+    bsub -oo ${PREFIX}_${QUEUE}_2x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 2 -R "span[ptile=1]" mpirun -np 2 ./swift_fixdt_mpi -r $RUNS -t $CPN -g "2 1 1" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
 fi
 
 if [ ! -e ${PREFIX}_${QUEUE}_4x${cpu}.dump ]
 then
-    bsub -oo ${PREFIX}_${QUEUE}_4x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 4 -R "span[ptile=1]" mpirun -np 4 ./test_fixdt_mpi -r $RUNS -t $CPN -g "2 2 1" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+    bsub -oo ${PREFIX}_${QUEUE}_4x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 4 -R "span[ptile=1]" mpirun -np 4 ./swift_fixdt_mpi -r $RUNS -t $CPN -g "2 2 1" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
 fi
 
 if [ ! -e ${PREFIX}_${QUEUE}_8x${cpu}.dump ]
 then
-    bsub -oo ${PREFIX}_${QUEUE}_8x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 8 -R "span[ptile=1]" mpirun -np 8 ./test_fixdt_mpi -r $RUNS -t $CPN -g "2 2 2" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+    bsub -oo ${PREFIX}_${QUEUE}_8x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 8 -R "span[ptile=1]" mpirun -np 8 ./swift_fixdt_mpi -r $RUNS -t $CPN -g "2 2 2" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
 fi
 
 if [ ! -e ${PREFIX}_${QUEUE}_16x${cpu}.dump ]
 then
-    bsub -oo ${PREFIX}_${QUEUE}_16x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 16 -R "span[ptile=1]" mpirun -np 16 ./test_fixdt_mpi -r $RUNS -t $CPN -g "4 2 2" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+    bsub -oo ${PREFIX}_${QUEUE}_16x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 16 -R "span[ptile=1]" mpirun -np 16 ./swift_fixdt_mpi -r $RUNS -t $CPN -g "4 2 2" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
 fi
 
 if [ ! -e ${PREFIX}_${QUEUE}_32x${cpu}.dump ]
 then
-    bsub -oo ${PREFIX}_${QUEUE}_32x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 32 -R "span[ptile=1]" mpirun -np 32 ./test_fixdt_mpi -r $RUNS -t $CPN -g "4 4 2" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+    bsub -oo ${PREFIX}_${QUEUE}_32x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 32 -R "span[ptile=1]" mpirun -np 32 ./swift_fixdt_mpi -r $RUNS -t $CPN -g "4 4 2" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
 fi
 
 if [ ! -e ${PREFIX}_${QUEUE}_64x${cpu}.dump ]
 then
-    bsub -oo ${PREFIX}_${QUEUE}_64x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 64 -R "span[ptile=1]" mpirun -np 64 ./test_fixdt_mpi -r $RUNS -t $CPN -g "4 4 4" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+    bsub -oo ${PREFIX}_${QUEUE}_64x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 64 -R "span[ptile=1]" mpirun -np 64 ./swift_fixdt_mpi -r $RUNS -t $CPN -g "4 4 4" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
 fi
 
 if [ ! -e ${PREFIX}_${QUEUE}_128x${cpu}.dump ]
 then
-    bsub -oo ${PREFIX}_${QUEUE}_128x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 128 -R "span[ptile=1]" mpirun -np 128 ./test_fixdt_mpi -r $RUNS -t $CPN -g "8 4 4" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
+    bsub -oo ${PREFIX}_${QUEUE}_128x${CPN}.dump -q ${QUEUE} -P ${PROJECT} -x -W 02:00 -n 128 -R "span[ptile=1]" mpirun -np 128 ./swift_fixdt_mpi -r $RUNS -t $CPN -g "8 4 4" -f ${INPUT} -m 0.705 -w 6000 -z 300 -d 1e-8
 fi
 
diff --git a/src/Makefile.am b/src/Makefile.am
index 46841c6394f5c6c024f1428a06d1025fa6604962..0fdc51571424e2a9ad91d1d7c714d91446702161 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,6 +1,6 @@
 
 # This file is part of SWIFT.
-# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
 #
 # This program is free software: you can redistribute it and/or modify
@@ -35,12 +35,13 @@ endif
 # List required headers
 include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
-    common_io.h single_io.h multipole.h
+    common_io.h single_io.h multipole.h map.h tools.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
-    units.c common_io.c single_io.c multipole.c version.c
+    units.c common_io.c single_io.c multipole.c version.c map.c \
+    kernel.c tools.c
 
 # Include files for distribution, not installation.
 noinst_HEADERS = atomic.h cycle.h error.h inline.h kernel.h vector.h \
diff --git a/src/atomic.h b/src/atomic.h
index 16b268c4c799cd1ca8c38a3382df912a9d618614..818d210e60a7aacdf61d12b60623ce87e62c9ed2 100644
--- a/src/atomic.h
+++ b/src/atomic.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/cell.c b/src/cell.c
index 0bcf25255bb9952ccb79a44ea337c5555f3c22ea..579028eca02be9f1f8ffa47c4806e330fdcc09e8 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
@@ -299,7 +299,7 @@ int cell_glocktree(struct cell *c) {
 }
 
 /**
- * @brief Unock a cell's parents.
+ * @brief Unlock a cell's parents.
  *
  * @param c The #cell.
  */
diff --git a/src/cell.h b/src/cell.h
index f147f927413faccb31111b0736daa2c978bdbdf0..4cc09bdecd8838fc579c4fb22ea28fd14a0e2416 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/common_io.c b/src/common_io.c
index a09f9c9ac37d280a4ae2838b3b7e36556e81047a..6ed22a582016f2bd6506a193a75c21988e1c5553 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
@@ -218,7 +218,7 @@ void writeStringAttribute(hid_t grp, char* name, const char* str, int length) {
 
 /**
  * @brief Writes a double value as an attribute
- * @param grp The groupm in which to write
+ * @param grp The group in which to write
  * @param name The name of the attribute
  * @param data The value to write
  */
@@ -228,7 +228,7 @@ void writeAttribute_d(hid_t grp, char* name, double data) {
 
 /**
  * @brief Writes a float value as an attribute
- * @param grp The groupm in which to write
+ * @param grp The group in which to write
  * @param name The name of the attribute
  * @param data The value to write
  */
@@ -238,7 +238,7 @@ void writeAttribute_f(hid_t grp, char* name, float data) {
 
 /**
  * @brief Writes an int value as an attribute
- * @param grp The groupm in which to write
+ * @param grp The group in which to write
  * @param name The name of the attribute
  * @param data The value to write
  */
@@ -249,7 +249,7 @@ void writeAttribute_i(hid_t grp, char* name, int data) {
 
 /**
  * @brief Writes a long value as an attribute
- * @param grp The groupm in which to write
+ * @param grp The group in which to write
  * @param name The name of the attribute
  * @param data The value to write
  */
@@ -259,7 +259,7 @@ void writeAttribute_l(hid_t grp, char* name, long data) {
 
 /**
  * @brief Writes a string value as an attribute
- * @param grp The groupm in which to write
+ * @param grp The group in which to write
  * @param name The name of the attribute
  * @param str The string to write
  */
diff --git a/src/common_io.h b/src/common_io.h
index c58ffae4d76dafe71faf44395f7f280120f4f0fe..c221ad3aaf84deb83c0067ee4ece729b98003430 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
diff --git a/src/const.h b/src/const.h
index ccccf6fa89884328efb33fcb018e0b17228fceff..d62673332c2d19872cad1d9c62ca6eff7f0e013a 100644
--- a/src/const.h
+++ b/src/const.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (ptcedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (ptcedro.gonnet@durham.ac.uk)
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
diff --git a/src/debug.c b/src/debug.c
index 06504455331c98f76a02ea724df1885236ee058d..0d8da8cae01ea00d10b094715e7c2d3387b3b2d7 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Matthieu Schaller (matthieu.schaller@durham.ac.uk),
+ * Copyright (c) 2013 Matthieu Schaller (matthieu.schaller@durham.ac.uk),
  *                    Pedro Gonnet (pedro.gonnet@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
diff --git a/src/debug.h b/src/debug.h
index 5c9eb7b32db8ef5d07d97b867c68ec09b848e337..83461df45e3c0fb137557fba5fdf68cac9d4915a 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2012 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
diff --git a/src/engine.c b/src/engine.c
index bc08795fab7db6dd327a8499bee936131761e9b9..f9cf36af3d5dd61df58a309537dd28dc67ff3bf7 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
@@ -127,7 +127,7 @@ void engine_mkghosts(struct engine *e, struct cell *c, struct cell *super) {
 }
 
 /**
- * @brief Redistribute the particles amongst the nodes accorind
+ * @brief Redistribute the particles amongst the nodes according
  *      to their cell's node IDs.
  *
  * @param e The #engine.
@@ -455,7 +455,7 @@ void engine_repartition(struct engine *e) {
     for (k = 0; k < nr_cells; k++) weights_v[k] *= scale;
   }
 
-/* Merge the weights arrays accross all nodes. */
+/* Merge the weights arrays across all nodes. */
 #if IDXTYPEWIDTH == 32
   if ((res = MPI_Reduce((nodeID == 0) ? MPI_IN_PLACE : weights_v, weights_v,
                         nr_cells, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD)) !=
@@ -578,7 +578,7 @@ void engine_repartition(struct engine *e) {
   /* Now comes the tricky part: Exchange particles between all nodes.
      This is done in two steps, first allreducing a matrix of
      how many particles go from where to where, then re-allocating
-     the parts array, and emiting the sends and receives.
+     the parts array, and emitting the sends and receives.
      Finally, the space, tasks, and proxies need to be rebuilt. */
 
   /* Redistribute the particles between the nodes. */
@@ -767,7 +767,7 @@ void engine_exchange_cells(struct engine *e) {
     proxy_cells_exch2(&e->proxies[pid]);
   }
 
-  /* Wait for all the sends to have finnished too. */
+  /* Wait for all the sends to have finished too. */
   if (MPI_Waitall(nr_proxies, reqs_out, MPI_STATUSES_IGNORE) != MPI_SUCCESS)
     error("MPI_Waitall on sends failed.");
 
@@ -788,7 +788,7 @@ void engine_exchange_cells(struct engine *e) {
                            e->proxies[pid].cells_in[j], e->s);
   }
 
-  /* Wait for all the sends to have finnished too. */
+  /* Wait for all the sends to have finished too. */
   if (MPI_Waitall(nr_proxies, reqs_out, MPI_STATUSES_IGNORE) != MPI_SUCCESS)
     error("MPI_Waitall on sends failed.");
 
@@ -885,11 +885,11 @@ int engine_exchange_strays(struct engine *e, int offset, int *ind, int N) {
     proxy_parts_exch2(&e->proxies[pid]);
   }
 
-  /* Wait for all the sends to have finnished too. */
+  /* Wait for all the sends to have finished too. */
   if (MPI_Waitall(e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) != MPI_SUCCESS)
     error("MPI_Waitall on sends failed.");
 
-  /* Count the total number of incomming particles and make sure we have
+  /* Count the total number of incoming particles and make sure we have
      enough space to accommodate them. */
   int count_in = 0;
   for (k = 0; k < e->nr_proxies; k++) count_in += e->proxies[k].nr_parts_in;
@@ -956,7 +956,7 @@ int engine_exchange_strays(struct engine *e, int offset, int *ind, int N) {
     }
   }
 
-  /* Wait for all the sends to have finnished too. */
+  /* Wait for all the sends to have finished too. */
   if (nr_out > 0)
     if (MPI_Waitall(2 * e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) !=
         MPI_SUCCESS)
@@ -1214,7 +1214,7 @@ void engine_maketasks(struct engine *e) {
     /* Get a handle on the proxy. */
     struct proxy *p = &e->proxies[pid];
 
-    /* Loop through the proxy's incomming cells and add the
+    /* Loop through the proxy's incoming cells and add the
        recv tasks. */
     for (k = 0; k < p->nr_cells_in; k++)
       engine_addtasks_recv(e, p->cells_in[k], NULL, NULL);
@@ -1252,7 +1252,7 @@ int engine_marktasks(struct engine *e) {
   struct cell *ci, *cj;
   // ticks tic = getticks();
 
-  /* Muc less to do here if we're on a fixed time-step. */
+  /* Much less to do here if we're on a fixed time-step. */
   if (!(e->policy & engine_policy_multistep)) {
 
     /* Run through the tasks and mark as skip or not. */
@@ -1455,7 +1455,7 @@ void engine_prepare(struct engine *e) {
   int buff;
   if (MPI_Allreduce(&rebuild, &buff, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD) !=
       MPI_SUCCESS)
-    error("Failed to aggreggate the rebuild flag accross nodes.");
+    error("Failed to aggregate the rebuild flag across nodes.");
   rebuild = buff;
 // message( "rebuild allreduce took %.3f ms." , (double)(getticks() -
 // tic)/CPU_TPS*1000 );
@@ -1506,7 +1506,7 @@ void engine_barrier(struct engine *e, int tid) {
   /* Wait for the barrier to open. */
   while (e->barrier_launch == 0 || tid >= e->barrier_launchcount)
     if (pthread_cond_wait(&e->barrier_cond, &e->barrier_mutex) != 0)
-      error("Eror waiting for barrier to close.");
+      error("Error waiting for barrier to close.");
 
   /* This thread has been launched. */
   e->barrier_running += 1;
@@ -1748,6 +1748,39 @@ void engine_step(struct engine *e) {
 
   TIMER_TIC2
 
+  /* Get the maximum dt. */
+  if (e->policy & engine_policy_multistep) {
+    dt_step = 2.0f * dt;
+    for (k = 0; k < 32 && (e->step & (1 << k)) == 0; k++) dt_step *= 2;
+  } else
+    dt_step = FLT_MAX;
+
+  /* Set the maximum dt. */
+  e->dt_step = dt_step;
+  e->s->dt_step = dt_step;
+  // message( "dt_step set to %.3e (dt=%.3e)." , dt_step , e->dt );
+  // fflush(stdout);
+
+  // printParticle( parts , 432626 );
+
+  /* First kick. */
+  if (e->step == 0 || !(e->policy & engine_policy_fixdt)) {
+    TIMER_TIC
+    engine_launch(e, (e->nr_threads > 8) ? 8 : e->nr_threads,
+                  (1 << task_type_kick1) | (1 << task_type_link));
+    TIMER_TOC(timer_kick1);
+  }
+
+  /* Check if all the kick1 threads have executed. */
+  /* for ( k = 0 ; k < e->sched.nr_tasks ; k++ )
+      if ( e->sched.tasks[k].type == task_type_kick1 &&
+           e->sched.tasks[k].toc == 0 )
+          error( "Not all kick1 tasks completed." ); */
+
+  // for(k=0; k<10; ++k)
+  //   printParticle(parts, k);
+  // printParticle( e->s->parts , 3392063069037 , e->s->nr_parts );
+
   /* Re-distribute the particles amongst the nodes? */
   if (e->forcerepart) engine_repartition(e);
 
@@ -1828,6 +1861,8 @@ if ( e->nodeID == 0 )
     message( "nr_parts=%i." , nr_parts ); */
 #endif
 
+  e->dt_min = dt_min;
+  e->dt_max = dt_max;
   e->count_step = count;
   e->ekin = ekin;
   e->epot = epot;
@@ -2055,7 +2090,7 @@ void engine_split(struct engine *e, int *grid) {
  * @param nr_queues The number of task queues to create.
  * @param nr_nodes The number of MPI ranks.
  * @param nodeID The MPI rank of this node.
- * @param policy The queueing policy to use.
+ * @param policy The queuing policy to use.
  * @param timeBegin Time at the begininning of the simulation.
  * @param timeEnd Time at the end of the simulation.
  */
@@ -2065,6 +2100,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
                  float timeBegin, float timeEnd) {
 
   int k;
+  float dt_min = dt;
 #if defined(HAVE_SETAFFINITY)
   int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
   int i, j, cpuid[nr_cores];
diff --git a/src/engine.h b/src/engine.h
index f3e28b5e4d38dc73ecef9eaba0c7f1163b69fcca..6ef59767a358ce950a2222201e1087352e995ec8 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/error.h b/src/error.h
index 0893b81757a0955975f6e01265c02315df94d81c..7e5c290203786bce33d244f0aa2dace621ec374d 100644
--- a/src/error.h
+++ b/src/error.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
diff --git a/src/inline.h b/src/inline.h
index 06728cb87f5e342b22d4a4a861cbd83ea6af31d9..c4dd9d59becabf54895b60cf3ac7ba809ac150c5 100644
--- a/src/inline.h
+++ b/src/inline.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
diff --git a/src/intrinsics.h b/src/intrinsics.h
index 4bb97e351f0a1e6dafe1c686e90f5038d98e7217..21b8c8e68bc45d8799db496ff30ac0cfb289acea 100644
--- a/src/intrinsics.h
+++ b/src/intrinsics.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2015 Matthieu Schaller matthieu.schaller@durham.ac.uk)
+ * 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
@@ -24,7 +24,7 @@
  * significant bit position. If x is 0, the result is undefined.
  *
  * This is a wrapper for the GCC intrinsics with an implementation (from
- * Hacker's Delight) if the compiler intrisincs are not available.
+ * Hacker's Delight) if the compiler intrinsics are not available.
  */
 __attribute__((always_inline))
     INLINE static int intrinsics_clz(unsigned int x) {
@@ -64,7 +64,7 @@ __attribute__((always_inline))
  * @brief Returns the number of 1-bits in x.
  *
  * This is a wrapper for the GCC intrinsics with an implementation (from
- * Hacker's Delight) if the compiler intrisincs are not available.
+ * Hacker's Delight) if the compiler intrinsics are not available.
  */
 __attribute__((always_inline))
     INLINE static int intrinsics_popcount(unsigned int x) {
diff --git a/src/kernel.c b/src/kernel.c
new file mode 100644
index 0000000000000000000000000000000000000000..58f5b0c9fdaa62663c65d5af18afe0a15a813834
--- /dev/null
+++ b/src/kernel.c
@@ -0,0 +1,99 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ *                    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 <math.h>
+#include <stdio.h>
+
+#include "kernel.h"
+
+/**
+ * @brief Test the SPH kernel function by dumping it in the interval [0,1].
+ *
+ * @param N number of intervals in [0,1].
+ */
+void SPH_kernel_dump(int N) {
+
+  int k;
+  float x, w, dw_dx;
+  float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+  float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+  // float dw_dx4[4] __attribute__ ((aligned (16)));
+
+  for (k = 0; k <= N; k++) {
+    x = ((float)k) / N;
+    x4[3] = x4[2];
+    x4[2] = x4[1];
+    x4[1] = x4[0];
+    x4[0] = x;
+    kernel_deval(x, &w, &dw_dx);
+    // kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
+    printf(" %e %e %e %e %e %e %e\n", x, w, dw_dx, w4[0], w4[1], w4[2], w4[3]);
+  }
+}
+
+/**
+ * @brief The Gadget-2 gravity kernel function
+ *
+ * @param r The distance between particles
+ */
+float gadget(float r) {
+  float fac, h_inv, u, r2 = r * r;
+  if (r >= const_epsilon)
+    fac = 1.0f / (r2 * r);
+  else {
+    h_inv = 1. / const_epsilon;
+    u = r * h_inv;
+    if (u < 0.5)
+      fac = const_iepsilon3 * (10.666666666667 + u * u * (32.0 * u - 38.4));
+    else
+      fac = const_iepsilon3 *
+            (21.333333333333 - 48.0 * u + 38.4 * u * u -
+             10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
+  }
+  return const_G * fac;
+}
+
+/**
+ * @brief Test the gravity kernel function by dumping it in the interval [0,1].
+ *
+ * @param r_max The radius up to which the kernel is dumped.
+ * @param N number of intervals in [0,1].
+ */
+void gravity_kernel_dump(float r_max, int N) {
+
+  int k;
+  float x, w;
+  float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+  float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+  // float dw_dx4[4] __attribute__ ((aligned (16)));
+
+  for (k = 1; k <= N; k++) {
+    x = (r_max * k) / N;
+    x4[3] = x4[2];
+    x4[2] = x4[1];
+    x4[1] = x4[0];
+    x4[0] = x;
+    kernel_grav_eval(x, &w);
+    w *= const_G / (x * x * x);
+    // blender_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
+    printf(" %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n", x, w * x, w4[0],
+           w4[1], w4[2], w4[3], gadget(x) * x);
+  }
+}
diff --git a/src/kernel.h b/src/kernel.h
index 0fc232597e1e9917d17f068407acc85b37659d42..d55e073479d1585f3677389939fa8e7c53dbe1c9 100644
--- a/src/kernel.h
+++ b/src/kernel.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
@@ -609,4 +609,8 @@ __attribute__((always_inline)) INLINE static void kernel_eval(float x,
 
 #endif  // Kernel choice
 
+/* Some cross-check functions */
+void SPH_kernel_dump(int N);
+void gravity_kernel_dump(float r_max, int N);
+
 #endif  // SWIFT_KERNEL_H
diff --git a/src/lock.h b/src/lock.h
index 19a4e74bf82d3b6bb8e305388ca42929cc9d719e..90e9f90602c120ddd10f4cdefb9b08cedbf45e0f 100644
--- a/src/lock.h
+++ b/src/lock.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/map.c b/src/map.c
new file mode 100644
index 0000000000000000000000000000000000000000..dcaa53465767d414cc54fe05940069ae5ff06d77
--- /dev/null
+++ b/src/map.c
@@ -0,0 +1,198 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ *                    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 <stdio.h>
+#include <stdlib.h>
+#include "error.h"
+#include "map.h"
+
+/**
+ * @brief Mapping function to draw a specific cell (gnuplot).
+ */
+
+void map_cells_plot(struct cell *c, void *data) {
+
+  int depth = *(int *)data;
+  double *l = c->loc, *h = c->h;
+
+  if (c->depth <= depth) {
+
+    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2]);
+    printf("%.16e %.16e %.16e\n\n\n", l[0], l[1] + h[1], l[2]);
+
+    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n\n\n", l[0], l[1] + h[1], l[2] + h[2]);
+
+    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n\n\n", l[0], l[1], l[2] + h[2]);
+
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1], l[2] + h[2]);
+
+    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0], l[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1], l[2]);
+
+    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2]);
+    printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
+    printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1] + h[1], l[2]);
+
+    if (!c->split) {
+      for (int k = 0; k < c->count; k++)
+        printf("0 0 0 %.16e %.16e %.16e\n", c->parts[k].x[0], c->parts[k].x[1],
+               c->parts[k].x[2]);
+      printf("\n\n");
+    }
+    /* else
+        for ( int k = 0 ; k < 8 ; k++ )
+            if ( c->progeny[k] != NULL )
+                map_cells_plot( c->progeny[k] , data ); */
+  }
+}
+
+/**
+ * @brief Mapping function for checking if each part is in its box.
+ */
+
+/* void map_check ( struct part *p , struct cell *c , void *data ) {
+
+    if ( p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
+         p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
+         p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] )
+        printf( "map_check: particle %i is outside of its box.\n" , p->id );
+
+    } */
+
+/**
+ * @brief Mapping function for neighbour count.
+ */
+
+void map_cellcheck(struct cell *c, void *data) {
+
+  int k, *count = (int *)data;
+  struct part *p;
+
+  __sync_fetch_and_add(count, c->count);
+
+  /* Loop over all parts and check if they are in the cell. */
+  for (k = 0; k < c->count; k++) {
+    p = &c->parts[k];
+    if (p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
+        p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] ||
+        p->x[2] > c->loc[2] + c->h[2]) {
+      printf(
+          "map_cellcheck: particle at [ %.16e %.16e %.16e ] outside of cell [ "
+          "%.16e %.16e %.16e ] - [ %.16e %.16e %.16e ].\n",
+          p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
+          c->loc[0] + c->h[0], c->loc[1] + c->h[1], c->loc[2] + c->h[2]);
+      error("particle out of bounds!");
+    }
+  }
+}
+
+/**
+ * @brief Mapping function for maxdepth cell count.
+ */
+
+void map_maxdepth(struct cell *c, void *data) {
+
+  int maxdepth = ((int *)data)[0];
+  int *count = &((int *)data)[1];
+
+  // printf( "%e\n" , p->count );
+
+  if (c->depth == maxdepth) *count += 1;
+}
+
+/**
+ * @brief Mapping function for neighbour count.
+ */
+
+void map_count(struct part *p, struct cell *c, void *data) {
+
+  double *wcount = (double *)data;
+
+  // printf( "%i %e %e\n" , p->id , p->count , p->count_dh );
+
+  *wcount += p->density.wcount;
+}
+
+void map_wcount_min(struct part *p, struct cell *c, void *data) {
+
+  struct part **p2 = (struct part **)data;
+
+  if (p->density.wcount < (*p2)->density.wcount) *p2 = p;
+}
+
+void map_wcount_max(struct part *p, struct cell *c, void *data) {
+
+  struct part **p2 = (struct part **)data;
+
+  if (p->density.wcount > (*p2)->density.wcount) *p2 = p;
+}
+
+void map_h_min(struct part *p, struct cell *c, void *data) {
+
+  struct part **p2 = (struct part **)data;
+
+  if (p->h < (*p2)->h) *p2 = p;
+}
+
+void map_h_max(struct part *p, struct cell *c, void *data) {
+
+  struct part **p2 = (struct part **)data;
+
+  if (p->h > (*p2)->h) *p2 = p;
+}
+
+/**
+ * @brief Mapping function for neighbour count.
+ */
+
+void map_icount(struct part *p, struct cell *c, void *data) {
+
+  // int *count = (int *)data;
+
+  // printf( "%i\n" , p->icount );
+
+  // *count += p->icount;
+}
+
+/**
+ * @brief Mapping function to print the particle position.
+ */
+
+void map_dump(struct part *p, struct cell *c, void *data) {
+
+  double *shift = (double *)data;
+
+  printf("%g\t%g\t%g\n", p->x[0] - shift[0], p->x[1] - shift[1],
+         p->x[2] - shift[2]);
+}
diff --git a/src/map.h b/src/map.h
new file mode 100644
index 0000000000000000000000000000000000000000..0fe04ad0158499626a16a0b18637ff190cf59c97
--- /dev/null
+++ b/src/map.h
@@ -0,0 +1,41 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ *                    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/>.
+ *
+ ******************************************************************************/
+
+
+/* Includes. */
+#ifndef SWIFT_MAP_H
+#define SWIFT_MAP_H
+
+#include "part.h"
+#include "cell.h"
+
+void map_cells_plot(struct cell *c, void *data);
+void map_check(struct part *p, struct cell *c, void *data);
+void map_cellcheck(struct cell *c, void *data);
+void map_maxdepth(struct cell *c, void *data);
+void map_count(struct part *p, struct cell *c, void *data);
+void map_wcount_min(struct part *p, struct cell *c, void *data);
+void map_wcount_max(struct part *p, struct cell *c, void *data);
+void map_h_min(struct part *p, struct cell *c, void *data);
+void map_h_max(struct part *p, struct cell *c, void *data);
+void map_icount(struct part *p, struct cell *c, void *data);
+void map_dump(struct part *p, struct cell *c, void *data);
+
+#endif /* SWIFT_MAP_H */
diff --git a/src/multipole.c b/src/multipole.c
index 439e9cd5f0218bddf28d228de6eb3bb14a2d6735..fc4701fead68e64e8742730325931282f80e8934 100644
--- a/src/multipole.c
+++ b/src/multipole.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/multipole.h b/src/multipole.h
index ffa5d713f507b85f6ae9216cfe81d4fc49316345..91ba6df965ce9d3b088d538411b7f0a8555ba0e4 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 63c9c26b08e497bc809304f8c46942babb4336a9..3fff596718ef41caba9c81a71e912b9d19ac192a 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
@@ -47,14 +47,14 @@
  * @param type The #DATA_TYPE of the attribute.
  * @param N The number of particles.
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part_c A (char*) pointer on the first occurence of the field of
+ * @param part_c A (char*) pointer on the first occurrence of the field of
  *interest in the parts array
  * @param importance If COMPULSORY, the data must be present in the IC file. If
  *OPTIONAL, the array will be zeroed when the data is not present.
  *
- * @todo A better version using HDF5 hyperslabs to read the file directly into
+ * @todo A better version using HDF5 hyper-slabs to read the file directly into
  *the part array
- * will be written once the strucutres have been stabilized.
+ * will be written once the structures have been stabilized.
  *
  * Calls #error() if an error occurs.
  */
@@ -102,7 +102,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
   temp = malloc(N * dim * typeSize);
   if (temp == NULL) error("Unable to allocate memory for temporary buffer");
 
-  /* Prepare information for hyperslab */
+  /* Prepare information for hyper-slab */
   if (dim > 1) {
     rank = 2;
     shape[0] = N;
@@ -120,7 +120,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
   /* Create data space in memory */
   h_memspace = H5Screate_simple(rank, shape, NULL);
 
-  /* Select hyperslab in file */
+  /* Select hyper-slab in file */
   h_filespace = H5Dget_space(h_data);
   H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);
 
@@ -185,7 +185,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
  * in the file.
  *
  * @warning Can not read snapshot distributed over more than 1 file !!!
- * @todo Read snaphsots distributed in more than one file.
+ * @todo Read snapshots distributed in more than one file.
  *
  * Calls #error() if an error occurs.
  *
@@ -309,14 +309,14 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
  * @param N_total Total number of particles across all cores
  * @param offset Offset in the array where this mpi task starts writing
- * @param part_c A (char*) pointer on the first occurence of the field of
+ * @param part_c A (char*) pointer on the first occurrence of the field of
  *interest in the parts array
  * @param us The UnitSystem currently in use
  * @param convFactor The UnitConversionFactor for this array
  *
- * @todo A better version using HDF5 hyperslabs to write the file directly from
+ * @todo A better version using HDF5 hyper-slabs to write the file directly from
  *the part array
- * will be written once the strucutres have been stabilized.
+ * will be written once the structures have been stabilized.
  *
  * Calls #error() if an error occurs.
  */
@@ -416,8 +416,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
   conversionString(buffer, us, convFactor);
   writeAttribute_d(h_data, "CGS conversion factor",
                    conversionFactor(us, convFactor));
-  writeAttribute_f(h_data, "h-scale exponant", hFactor(us, convFactor));
-  writeAttribute_f(h_data, "a-scale exponant", aFactor(us, convFactor));
+  writeAttribute_f(h_data, "h-scale exponent", hFactor(us, convFactor));
+  writeAttribute_f(h_data, "a-scale exponent", aFactor(us, convFactor));
   writeAttribute_s(h_data, "Conversion factor", buffer);
 
   /* Free and close everything */
@@ -441,7 +441,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
  * @param N_total Total number of particles across all cores
  * @param mpi_rank The MPI task rank calling the function
  * @param offset Offset in the array where this mpi task starts writing
- * @param part A (char*) pointer on the first occurence of the field of interest
+ * @param part A (char*) pointer on the first occurrence of the field of interest
  *in the parts array
  * @param field The name (code name) of the field to read from.
  * @param us The UnitSystem currently in use
@@ -511,7 +511,7 @@ void write_output_parallel(struct engine* e, struct UnitSystem* us,
   MPI_Bcast(&N_total_d, 1, MPI_DOUBLE, mpi_size - 1, comm);
   if (N_total_d > 1.e15)
     error(
-        "Error while computing the offest for parallel output: Simulation has "
+        "Error while computing the offset for parallel output: Simulation has "
         "more than 10^15 particles.\n");
   N_total = (long long)N_total_d;
   offset = (long long)offset_d;
diff --git a/src/parallel_io.h b/src/parallel_io.h
index fa46a230ab73e52a3b471dc0b157f5cf0f99ef73..5ee2714a18dbabdd5fc8c07c0b444e04a87ad4f8 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2012 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
diff --git a/src/part.h b/src/part.h
index d0e29d3c45218bb7d90076dd6d87e85472167098..97977dadc988f361c7c6aae5b512ca857ad384da 100644
--- a/src/part.h
+++ b/src/part.h
@@ -1,7 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- *                    Matthieu Schaller (matthieu.schaller@durhm.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/proxy.c b/src/proxy.c
index 257a11593d674a3f8a5674670ac3a12f203e1eb2..7d2e546bf945ca18c2195ea2801d1b2058cb2f58 100644
--- a/src/proxy.c
+++ b/src/proxy.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@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
@@ -142,7 +142,7 @@ void proxy_addcell_in(struct proxy *p, struct cell *c) {
   if (p->nr_cells_in == p->size_cells_in) {
     p->size_cells_in *= proxy_buffgrow;
     if ((temp = malloc(sizeof(struct cell *) * p->size_cells_in)) == NULL)
-      error("Failed to allocate ingoing cell list.");
+      error("Failed to allocate incoming cell list.");
     memcpy(temp, p->cells_in, sizeof(struct cell *) * p->nr_cells_in);
     free(p->cells_in);
     p->cells_in = temp;
diff --git a/src/proxy.h b/src/proxy.h
index 8cb08d0a66095597227b52b317f3808190cdc45f..3cd33e0f0819ee1ecac53213630445b39c809dea 100644
--- a/src/proxy.h
+++ b/src/proxy.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@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
@@ -40,7 +40,7 @@ struct proxy {
   /* ID of the node this proxy represents. */
   int mynodeID, nodeID;
 
-  /* Incomming cells. */
+  /* Incoming cells. */
   struct cell **cells_in;
   struct pcell *pcells_in;
   int nr_cells_in, size_cells_in, size_pcells_in;
diff --git a/src/queue.c b/src/queue.c
index 3fa0096bf0fab8ecc6ec2508d5a7c2529451e54d..4e5330518c3d067a3e8a8ce7806e3560a41cc8ef 100644
--- a/src/queue.c
+++ b/src/queue.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/queue.h b/src/queue.h
index 533007684fa41a4a25a10a14c504358926d0fe06..d81627b4c4a709498c31c6c6152d3650388020e2 100644
--- a/src/queue.h
+++ b/src/queue.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/runner.c b/src/runner.c
index 999f65e2e6a74655ba1bddf4ed0533224e6c6d77..c4835a311f4406daf8777e919c6a921ab2b978c6 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
@@ -599,9 +599,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
             p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
 
         /* If no derivative, double the smoothing length. */
-        if (wcount_dh == 0.0f) {
-          h_corr = p->h;
-        }
+        if (wcount_dh == 0.0f) h_corr = p->h;
+
         /* Otherwise, compute the smoothing length update (Newton step). */
         else {
           h_corr = (kernel_nwneigh - wcount) / wcount_dh;
@@ -693,6 +692,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
     count = redo;
     if (count > 0) {
 
+      // error( "Bad smoothing length, fixing this isn't implemented yet." );
+
       /* Climb up the cell hierarchy. */
       for (finger = c; finger != NULL; finger = finger->parent) {
 
@@ -758,6 +759,131 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
 
   TIMER_TIC
 
+  /* Init idt to avoid compiler stupidity. */
+  idt = (dt > 0) ? 1.0f / dt : 0.0f;
+  hdt = dt / 2;
+
+  /* Loop over the particles and kick them. */
+  for (k = 0; k < nr_parts; k++) {
+
+    /* Get a handle on the part. */
+    p = &parts[k];
+    xp = &xparts[k];
+
+    /* Get local copies of particle data. */
+    pdt = p->dt;
+    m = p->mass;
+    x[0] = p->x[0];
+    x[1] = p->x[1];
+    x[2] = p->x[2];
+    v_hdt[0] = xp->v_hdt[0];
+    v_hdt[1] = xp->v_hdt[1];
+    v_hdt[2] = xp->v_hdt[2];
+    u_hdt = xp->u_hdt;
+
+    /* Update the particle's data (if active). */
+    if (pdt <= dt_step) {
+
+      /* Increase the number of particles updated. */
+      count += 1;
+
+      /* Scale the derivatives as they're freshly computed. */
+      h = p->h;
+      h_dt = p->force.h_dt *= h * 0.333333333f;
+      xp->omega = 1.0f + h * p->rho_dh / p->rho * 0.3333333333f;
+
+      /* Compute the new time step. */
+      u_dt = p->force.u_dt;
+      dt_cfl = const_cfl * h / p->force.v_sig;
+      dt_h_change =
+          (h_dt != 0.0f) ? fabsf(const_ln_max_h_change * h / h_dt) : FLT_MAX;
+      dt_u_change =
+          (u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / u_dt) : FLT_MAX;
+      dt_new = fminf(dt_cfl, fminf(dt_h_change, dt_u_change));
+      if (pdt == 0.0f)
+        p->dt = pdt = dt_new;
+      else
+        p->dt = pdt = fminf(dt_new, 2.0f * pdt);
+
+      /* Update positions and energies at the full step. */
+      p->v[0] = v_hdt[0] + hdt * p->a[0];
+      p->v[1] = v_hdt[1] + hdt * p->a[1];
+      p->v[2] = v_hdt[2] + hdt * p->a[2];
+      p->u = u_hdt + hdt * u_dt;
+
+      /* Set the new particle-specific time step. */
+      if (dt > 0.0f) {
+        float dt_curr = dt;
+        j = (int)(pdt * idt);
+        while (j > 1) {
+          dt_curr *= 2.0f;
+          j >>= 1;
+        }
+        xp->dt_curr = dt_curr;
+      }
+    }
+
+    /* Get the smallest/largest dt. */
+    dt_min = fminf(dt_min, pdt);
+    dt_max = fmaxf(dt_max, pdt);
+
+    /* Collect total energy. */
+    ekin += 0.5 * m *
+            (v_hdt[0] * v_hdt[0] + v_hdt[1] * v_hdt[1] + v_hdt[2] * v_hdt[2]);
+    epot += m * u_hdt;
+
+    /* Collect momentum */
+    mom[0] += m * v_hdt[0];
+    mom[1] += m * v_hdt[1];
+    mom[2] += m * v_hdt[2];
+
+    /* Collect angular momentum */
+    ang[0] += m * (x[1] * v_hdt[2] - x[2] * v_hdt[1]);
+    ang[1] += m * (x[2] * v_hdt[0] - x[0] * v_hdt[2]);
+    ang[2] += m * (x[0] * v_hdt[1] - x[1] * v_hdt[0]);
+
+    /* Collect entropic function */
+    // lent += u * pow( p->rho, 1.f-const_gamma );
+  }
+
+#ifdef TIMER_VERBOSE
+  message("runner %02i: %i parts at depth %i took %.3f ms.", r->id, c->count,
+          c->depth, ((double)TIMER_TOC(timer_kick2)) / CPU_TPS * 1000);
+  fflush(stdout);
+#else
+  TIMER_TOC(timer_kick2);
+#endif
+
+  /* Store the computed values in the cell. */
+  c->dt_min = dt_min;
+  c->dt_max = dt_max;
+  c->updated = count;
+  c->ekin = ekin;
+  c->epot = epot;
+  c->mom[0] = mom[0];
+  c->mom[1] = mom[1];
+  c->mom[2] = mom[2];
+  c->ang[0] = ang[0];
+  c->ang[1] = ang[1];
+  c->ang[2] = ang[2];
+}
+
+/**
+ * @brief Mapping function to set dt_min and dt_max, do the first
+ * kick.
+ */
+
+void runner_dokick1(struct runner *r, struct cell *c) {
+
+  int j, k;
+  struct engine *e = r->e;
+  float pdt, dt_step = e->dt_step, dt = e->dt, hdt = dt / 2;
+  float dt_min, dt_max, h_max, dx, dx_max;
+  float a[3], v[3], u, u_dt, h, h_dt, w, rho;
+  double x[3], x_old[3];
+  struct part *restrict p, *restrict parts = c->parts;
+  struct xpart *restrict xp, *restrict xparts = c->xparts;
+
   /* No children? */
   if (!c->split) {
 
@@ -817,6 +943,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
       /* Predict gradient term */
       p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (rho * xp->omega);
     }
+
   }
 
   if (timer) {
@@ -828,6 +955,12 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
     TIMER_TOC(timer_drift);
 #endif
   }
+
+  /* Store the values. */
+  c->dt_min = dt_min;
+  c->dt_max = dt_max;
+  c->h_max = h_max;
+  c->dx_max = dx_max;
 }
 
 /**
@@ -860,6 +993,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
   if (!c->split) {
 
     /* Init the min/max counters. */
+    dt_min = FLT_MAX;
+    dt_max = 0.0f;
     h_max = 0.0f;
     dx_max = 0.0f;
 
@@ -869,6 +1004,11 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
       /* Get a handle on the part. */
       p = &parts[k];
       xp = &xparts[k];
+
+      /* Get local copies of particle data. */
+      pdt = p->dt;
+      u_dt = p->force.u_dt;
+      h = p->h;
       m = p->mass;
       x[0] = p->x[0], x[1] = p->x[1], x[2] = p->x[2];
 
@@ -964,6 +1104,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
   else {
 
     /* Init with the first non-null child. */
+    dt_min = FLT_MAX;
+    dt_max = 0.0f;
     h_max = 0.0f;
     dx_max = 0.0f;
     updated = 0;
@@ -981,6 +1123,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
       if (c->progeny[k] != NULL) {
         struct cell *cp = c->progeny[k];
         runner_dokick(r, cp, 0);
+        dt_min = fminf(dt_min, cp->dt_min);
+        dt_max = fmaxf(dt_max, cp->dt_max);
         h_max = fmaxf(h_max, cp->h_max);
         dx_max = fmaxf(dx_max, cp->dx_max);
         updated += cp->count;
@@ -998,6 +1142,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
   }
 
   /* Store the values. */
+  c->dt_min = dt_min;
+  c->dt_max = dt_max;
   c->h_max = h_max;
   c->dx_max = dx_max;
   c->updated = count;
diff --git a/src/runner.h b/src/runner.h
index b9c2ed9379ccb8f04058367a9366e6a5ff87b3b2..2de7a5931c5af4c7c52de1b024cfa2dd41d1f603 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index bc3f5fe0f8fdfb723a2264a3a34c1163c6338e1b..0c319243131243be07d8c842a97d0dae9b640be3 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
@@ -1299,7 +1299,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
 
           } else {
 
-            /* Add this interaction to the non-summetric queue. */
+            /* Add this interaction to the non-symmetric queue. */
             r2q1[icount1] = r2;
             dxq1[3 * icount1 + 0] = dx[0];
             dxq1[3 * icount1 + 1] = dx[1];
@@ -1818,7 +1818,7 @@ void DOSUB1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
 
     }
 
-    /* Otherwsie, compute self-interaction. */
+    /* Otherwise, compute self-interaction. */
     else
       DOSELF1(r, ci);
 
@@ -2093,7 +2093,7 @@ void DOSUB2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
 
     }
 
-    /* Otherwsie, compute self-interaction. */
+    /* Otherwise, compute self-interaction. */
     else
       DOSELF2(r, ci);
 
@@ -2382,7 +2382,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
 
     }
 
-    /* Otherwsie, compute self-interaction. */
+    /* Otherwise, compute self-interaction. */
     else
       DOSELF_SUBSET(r, ci, parts, ind, count);
 
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index ac7267f6e036757310622040d56c2e1ee0213c4f..85329c95a4723cf01b92a52b5500cac05020b876 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@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
@@ -241,7 +241,7 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
 
 /**
  * @brief Compute the recursive downward sweep, i.e. apply the multipole
- *        accelleration on all the particles.
+ *        acceleration on all the particles.
  *
  * @param r The #runner.
  * @param c The top-level #cell.
@@ -254,7 +254,7 @@ void runner_dograv_down(struct runner *r, struct cell *c) {
   /* Split? */
   if (c->split) {
 
-    /* Apply this cell's accelleration on the multipoles below. */
+    /* Apply this cell's acceleration on the multipoles below. */
     for (int k = 0; k < 8; k++)
       if (c->progeny[k] != NULL) {
         struct multipole *mp = &c->progeny[k]->multipole;
@@ -272,7 +272,7 @@ void runner_dograv_down(struct runner *r, struct cell *c) {
   /* No, leaf node. */
   else {
 
-    /* Apply the multipole accelleration to all gparts. */
+    /* Apply the multipole acceleration to all gparts. */
     for (int k = 0; k < c->gcount; k++) {
       struct gpart *p = &c->gparts[k];
       p->a[0] += m->a[0];
@@ -589,7 +589,7 @@ void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj,
     theta = (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) /
             (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]);
 
-    /* Split the interacton? */
+    /* Split the interaction? */
     if (theta < const_theta_max * const_theta_max) {
 
       /* Are both ci and cj split? */
diff --git a/src/runner_iact.h b/src/runner_iact.h
index e1561132af3fff847989af34c268dfb1069ed40d..f1828eaf106b0ee94435fcaa90ea1feb3485b5d8 100644
--- a/src/runner_iact.h
+++ b/src/runner_iact.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
@@ -421,7 +421,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Apply balsara switch */
   Pi_ij *= (pi->force.balsara + pj->force.balsara);
 
-  /* Termal conductivity */
+  /* Thermal conductivity */
   v_sig_u = sqrtf(2.f * (const_hydro_gamma - 1.f) *
                   fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
   tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
@@ -622,7 +622,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
   Pi_ij.v = balsara.v * alpha_ij.v * v_sig.v * omega_ij.v / (pirho.v + pjrho.v);
   Pi_ij.v *= (wi_dr.v + wj_dr.v);
 
-  /* Termal conductivity */
+  /* Thermal conductivity */
   v_sig_u.v = vec_sqrt(vec_set1(2.f * (const_hydro_gamma - 1.f)) *
                        vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) /
                        (pirho.v + pjrho.v));
@@ -737,7 +737,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   /* Apply balsara switch */
   Pi_ij *= (pi->force.balsara + pj->force.balsara);
 
-  /* Termal conductivity */
+  /* Thermal conductivity */
   v_sig_u = sqrtf(2.f * (const_hydro_gamma - 1.f) *
                   fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
   tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
@@ -929,7 +929,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
   Pi_ij.v = balsara.v * alpha_ij.v * v_sig.v * omega_ij.v / (pirho.v + pjrho.v);
   Pi_ij.v *= (wi_dr.v + wj_dr.v);
 
-  /* Termal conductivity */
+  /* Thermal conductivity */
   v_sig_u.v = vec_sqrt(vec_set1(2.f * (const_hydro_gamma - 1.f)) *
                        vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) /
                        (pirho.v + pjrho.v));
diff --git a/src/runner_iact_grav.h b/src/runner_iact_grav.h
index 2fd30c1c3854db56564300f0a3e1a13a6dc31251..e62be446e8263bf02e3fd73f902b28cb1c3b16cf 100644
--- a/src/runner_iact_grav.h
+++ b/src/runner_iact_grav.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
@@ -53,7 +53,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav(
   /* Scale the acceleration. */
   acc *= const_G * ir * ir * ir;
 
-  /* Aggregate the accellerations. */
+  /* Aggregate the accelerations. */
   for (k = 0; k < 3; k++) {
     w = acc * dx[k];
     pi->a[k] -= w * mj;
@@ -101,7 +101,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_grav(
   /* Scale the acceleration. */
   acc.v *= vec_set1(const_G) * ir.v * ir.v * ir.v;
 
-  /* Aggregate the accellerations. */
+  /* Aggregate the accelerations. */
   for (k = 0; k < 3; k++) {
     w.v = acc.v * dx[k].v;
     ai.v = w.v * mj.v;
diff --git a/src/runner_iact_legacy.h b/src/runner_iact_legacy.h
index 42044caafc8686f4ce5f639bacc16f297f266868..5716cf4a7a84a9d9fcd4cef38a307e483389d4ab 100644
--- a/src/runner_iact_legacy.h
+++ b/src/runner_iact_legacy.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
diff --git a/src/scheduler.c b/src/scheduler.c
index a70319c64e66b449ceed6eb0ce7805ef80527563..bbf1d11447b5d77e4fac24e6d59869e450e030a8 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
@@ -186,7 +186,7 @@ void scheduler_splittasks(struct scheduler *s) {
           /* Take a step back (we're going to recycle the current task)... */
           redo = 1;
 
-          /* Add the self taks. */
+          /* Add the self task. */
           for (k = 0; ci->progeny[k] == NULL; k++)
             ;
           t->ci = ci->progeny[k];
@@ -608,7 +608,7 @@ void scheduler_splittasks(struct scheduler *s) {
           /* Ignore this task if the cell has no gparts. */
           if (ci->gcount == 0 || cj->gcount == 0) t->type = task_type_none;
 
-          /* Split the interacton? */
+          /* Split the interaction? */
           else if (theta < const_theta_max * const_theta_max) {
 
             /* Are both ci and cj split? */
@@ -716,7 +716,7 @@ void scheduler_ranktasks(struct scheduler *s) {
   struct task *t, *tasks = s->tasks;
   int *tid = s->tasks_ind, nr_tasks = s->nr_tasks;
 
-  /* Run throught the tasks and get all the waits right. */
+  /* Run through the tasks and get all the waits right. */
   for (i = 0, k = 0; k < nr_tasks; k++) {
     tid[k] = k;
     for (j = 0; j < tasks[k].nr_unlock_tasks; j++)
@@ -774,7 +774,7 @@ void scheduler_reset(struct scheduler *s, int size) {
   /* Do we need to re-allocate? */
   if (size > s->size) {
 
-    /* Free exising task lists if necessary. */
+    /* Free existing task lists if necessary. */
     if (s->tasks != NULL) free(s->tasks);
     if (s->tasks_ind != NULL) free(s->tasks_ind);
 
@@ -813,7 +813,7 @@ void scheduler_reweight(struct scheduler *s) {
   float wscale = 0.001;
   // ticks tic;
 
-  /* Run throught the tasks backwards and set their waits and
+  /* Run through the tasks backwards and set their waits and
      weights. */
   // tic = getticks();
   for (k = nr_tasks - 1; k >= 0; k--) {
@@ -902,7 +902,7 @@ void scheduler_start(struct scheduler *s, unsigned int mask) {
   struct task *t, *tasks = s->tasks;
   // ticks tic;
 
-  /* Run throught the tasks and set their waits. */
+  /* Run through the tasks and set their waits. */
   // tic = getticks();
   for (k = nr_tasks - 1; k >= 0; k--) {
     t = &tasks[tid[k]];
@@ -991,7 +991,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
           MPI_Error_string(err, buff, &len);
           error("Failed to emit irecv for particle data (%s).", buff);
         }
-        // message( "recieving %i parts with tag=%i from %i to %i." ,
+        // message( "receiving %i parts with tag=%i from %i to %i." ,
         //     t->ci->count , t->flags , t->ci->nodeID , s->nodeID );
         // fflush(stdout);
         qid = 1 % s->nr_queues;
@@ -1130,7 +1130,7 @@ struct task *scheduler_unlock(struct scheduler *s, struct task *t) {
  * @brief Get a task, preferably from the given queue.
  *
  * @param s The #scheduler.
- * @param qid The ID of the prefered #queue.
+ * @param qid The ID of the preferred #queue.
  * @param super the super-cell
  *
  * @return A pointer to a #task or @c NULL if there are no available tasks.
@@ -1161,7 +1161,7 @@ struct task *scheduler_gettask(struct scheduler *s, int qid,
         if (res != NULL) break;
       }
 
-      /* If unsucessful, try stealing from the other queues. */
+      /* If unsuccessful, try stealing from the other queues. */
       if (s->flags & scheduler_flag_steal) {
         int count = 0, qids[nr_queues];
         for (k = 0; k < nr_queues; k++)
diff --git a/src/scheduler.h b/src/scheduler.h
index 620b712885a1653397b3e9fd0e632cc0e562cf19..fd30c76f56ec7ffceba8af10a96127ca2ffc0135 100644
--- a/src/scheduler.h
+++ b/src/scheduler.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/serial_io.c b/src/serial_io.c
index c98135d24055c1a69bf836172cd3b2518703bef7..c61774fed2e5991e5a1ed06bac572224b09ade9c 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
@@ -55,14 +55,14 @@
  * @param type The #DATA_TYPE of the attribute.
  * @param N The number of particles.
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part_c A (char*) pointer on the first occurence of the field of
+ * @param part_c A (char*) pointer on the first occurrence of the field of
  *interest in the parts array
  * @param importance If COMPULSORY, the data must be present in the IC file. If
  *OPTIONAL, the array will be zeroed when the data is not present.
  *
- * @todo A better version using HDF5 hyperslabs to read the file directly into
+ * @todo A better version using HDF5 hyper-slabs to read the file directly into
  *the part array
- * will be written once the strucutres have been stabilized.
+ * will be written once the structures have been stabilized.
  *
  * Calls #error() if an error occurs.
  */
@@ -110,7 +110,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
   temp = malloc(N * dim * typeSize);
   if (temp == NULL) error("Unable to allocate memory for temporary buffer");
 
-  /* Prepare information for hyperslab */
+  /* Prepare information for hyper-slab */
   if (dim > 1) {
     rank = 2;
     shape[0] = N;
@@ -128,7 +128,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
   /* Create data space in memory */
   h_memspace = H5Screate_simple(rank, shape, NULL);
 
-  /* Select hyperslab in file */
+  /* Select hyper-slab in file */
   h_filespace = H5Dget_space(h_data);
   H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);
 
@@ -202,7 +202,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
  * in the file.
  *
  * @warning Can not read snapshot distributed over more than 1 file !!!
- * @todo Read snaphsots distributed in more than one file.
+ * @todo Read snapshots distributed in more than one file.
  *
  * Calls #error() if an error occurs.
  *
@@ -227,7 +227,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
     /* message("Opening file '%s' as IC.", fileName); */
     h_file = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT);
     if (h_file < 0)
-      error("Error while opening file '%s' for inital read.", fileName);
+      error("Error while opening file '%s' for initial read.", fileName);
 
     /* Open header to read simulation properties */
     /* message("Reading runtime parameters..."); */
@@ -382,8 +382,8 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, char* name,
   conversionString(buffer, us, convFactor);
   writeAttribute_d(h_data, "CGS conversion factor",
                    conversionFactor(us, convFactor));
-  writeAttribute_f(h_data, "h-scale exponant", hFactor(us, convFactor));
-  writeAttribute_f(h_data, "a-scale exponant", aFactor(us, convFactor));
+  writeAttribute_f(h_data, "h-scale exponent", hFactor(us, convFactor));
+  writeAttribute_f(h_data, "a-scale exponent", aFactor(us, convFactor));
   writeAttribute_s(h_data, "Conversion factor", buffer);
 
   H5Dclose(h_data);
@@ -400,7 +400,7 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, char* name,
  * @param type The #DATA_TYPE of the array.
  * @param N The number of particles to write.
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part_c A (char*) pointer on the first occurence of the field of
+ * @param part_c A (char*) pointer on the first occurrence of the field of
  *interest in the parts array
  * @param us The UnitSystem currently in use
  * @param convFactor The UnitConversionFactor for this array
@@ -431,7 +431,7 @@ void writeArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
   for (i = 0; i < N; ++i)
     memcpy(&temp_c[i * copySize], part_c + i * partSize, copySize);
 
-  /* Construct information for the hyperslab */
+  /* Construct information for the hyper-slab */
   if (dim > 1) {
     rank = 2;
     shape[0] = N;
@@ -487,7 +487,7 @@ void writeArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
  * @param type The #DATA_TYPE of the array.
  * @param N The number of particles to write.
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part A (char*) pointer on the first occurence of the field of interest
+ * @param part A (char*) pointer on the first occurrence of the field of interest
  *in the parts array
  * @param field The name (code name) of the field to read from.
  * @param us The UnitSystem currently in use
@@ -540,7 +540,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
   MPI_Bcast(&N_total_d, 1, MPI_DOUBLE, mpi_size - 1, comm);
   if (N_total_d > 1.e15)
     error(
-        "Error while computing the offest for parallel output: Simulation has "
+        "Error while computing the offset for parallel output: Simulation has "
         "more than 10^15 particles.\n");
   N_total = (long long)N_total_d;
   offset = (long long)offset_d;
diff --git a/src/serial_io.h b/src/serial_io.h
index bb05fc61bdca1b0db36386e6773a034cc17ea6b9..e5ecca9c8cbafbbaf2e555c5c216b494e25cc922 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2012 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
diff --git a/src/single_io.c b/src/single_io.c
index f26563d3e548af279a98f46a1b0e858114d3c25c..7e92949e24f6618861050597e2e9335c26c8e1b8 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
@@ -50,14 +50,14 @@
  * @param type The #DATA_TYPE of the attribute.
  * @param N The number of particles.
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part_c A (char*) pointer on the first occurence of the field of
+ * @param part_c A (char*) pointer on the first occurrence of the field of
  *interest in the parts array
  * @param importance If COMPULSORY, the data must be present in the IC file. If
  *OPTIONAL, the array will be zeroed when the data is not present.
  *
- * @todo A better version using HDF5 hyperslabs to read the file directly into
+ * @todo A better version using HDF5 hyper-slabs to read the file directly into
  *the part array
- * will be written once the strucutres have been stabilized.
+ * will be written once the structures have been stabilized.
  *
  * Calls #error() if an error occurs.
  */
@@ -158,7 +158,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
  * in the file.
  *
  * @warning Can not read snapshot distributed over more than 1 file !!!
- * @todo Read snaphsots distributed in more than one file.
+ * @todo Read snapshots distributed in more than one file.
  *
  * Calls #error() if an error occurs.
  *
@@ -256,14 +256,14 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
  * @param type The #DATA_TYPE of the array.
  * @param N The number of particles to write.
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part_c A (char*) pointer on the first occurence of the field of
+ * @param part_c A (char*) pointer on the first occurrence of the field of
  *interest in the parts array
  * @param us The UnitSystem currently in use
  * @param convFactor The UnitConversionFactor for this array
  *
- * @todo A better version using HDF5 hyperslabs to write the file directly from
+ * @todo A better version using HDF5 hyper-slabs to write the file directly from
  *the part array
- * will be written once the strucutres have been stabilized.
+ * will be written once the structures have been stabilized.
  *
  * Calls #error() if an error occurs.
  */
@@ -333,8 +333,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
   conversionString(buffer, us, convFactor);
   writeAttribute_d(h_data, "CGS conversion factor",
                    conversionFactor(us, convFactor));
-  writeAttribute_f(h_data, "h-scale exponant", hFactor(us, convFactor));
-  writeAttribute_f(h_data, "a-scale exponant", aFactor(us, convFactor));
+  writeAttribute_f(h_data, "h-scale exponent", hFactor(us, convFactor));
+  writeAttribute_f(h_data, "a-scale exponent", aFactor(us, convFactor));
   writeAttribute_s(h_data, "Conversion factor", buffer);
 
   /* Free and close everything */
@@ -353,7 +353,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
  * @param type The #DATA_TYPE of the array.
  * @param N The number of particles to write.
  * @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part A (char*) pointer on the first occurence of the field of interest
+ * @param part A (char*) pointer on the first occurrence of the field of interest
  *in the parts array
  * @param field The name (code name) of the field to read from.
  * @param us The UnitSystem currently in use
diff --git a/src/single_io.h b/src/single_io.h
index 91d229178bbd45df3ba358172d1f52c70008adb7..f6689901106a2cd5d85a873d4047e1d21edd3547 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2012 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
diff --git a/src/space.c b/src/space.c
index 603f7aec42bfc7048f61146191cc1ea5af568207..2c0e6a7ec36116b28f29c1f728e50b5149eaf651 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
@@ -179,7 +179,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
     float buff;
     if (MPI_Allreduce(&h_max, &buff, 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD) !=
         MPI_SUCCESS)
-      error("Failed to aggreggate the rebuild flag accross nodes.");
+      error("Failed to aggregate the rebuild flag across nodes.");
     h_max = buff;
   }
 #endif
@@ -993,7 +993,7 @@ void space_split(struct space *s, struct cell *c) {
     c->t_end_max = t_end_max;
   }
 
-  /* Set ownership accorind to the start of the parts array. */
+  /* Set ownership according to the start of the parts array. */
   c->owner = ((c->parts - s->parts) % s->nr_parts) * s->nr_queues / s->nr_parts;
 }
 
@@ -1089,7 +1089,7 @@ struct cell *space_getcell(struct space *s) {
 void space_init(struct space *s, double dim[3], struct part *parts, int N,
                 int periodic, double h_max, int verbose) {
 
-  /* Store eveything in the space. */
+  /* Store everything in the space. */
   s->dim[0] = dim[0];
   s->dim[1] = dim[1];
   s->dim[2] = dim[2];
diff --git a/src/space.h b/src/space.h
index c12ec46be968d713618c41db5ab1385ed147d33e..cace0a6d6d69e35990309934586c63d0f402f814 100644
--- a/src/space.h
+++ b/src/space.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/swift.h b/src/swift.h
index b302bca9b007cec47c96e1ab07770a3a3dc84966..1b7d86c769bb33f2677d81a25c286d88cbc2782c 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
@@ -32,6 +32,7 @@
 #include "engine.h"
 #include "error.h"
 #include "lock.h"
+#include "map.h"
 #include "multipole.h"
 #include "parallel_io.h"
 #include "part.h"
@@ -44,6 +45,7 @@
 #include "task.h"
 #include "timers.h"
 #include "units.h"
+#include "tools.h"
 #include "version.h"
 
 #ifdef LEGACY_GADGET2_SPH
diff --git a/src/task.c b/src/task.c
index 949caab56c4c4d8a0e3c73d05014ebc5ad68657a..642e6833c48236395f20c20b8e83f93d74757ea0 100644
--- a/src/task.c
+++ b/src/task.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/task.h b/src/task.h
index 852b376fb42fc1c08fd4de1e888ec4dedf3ee751..6d9a3f6b2c90c38b3a55b4f15e833892fbeccf78 100644
--- a/src/task.h
+++ b/src/task.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/timers.c b/src/timers.c
index 01a77d7804241f108b092f7d6857c90be3861cd0..2501d347c8cea608650ece4c2883dab85ceee058 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/timers.h b/src/timers.h
index e6fc55d98876c4704554171e45807ada818bebe0..4ae56e5e672fa4a5da947312f424c7cef8e61fd2 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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
diff --git a/src/tools.c b/src/tools.c
new file mode 100644
index 0000000000000000000000000000000000000000..5f8abcbeb6d339bedc3df71e5ccf651535b8bc0c
--- /dev/null
+++ b/src/tools.c
@@ -0,0 +1,285 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Copyright (c) 2015 Peter W. Draper (p.w.draper@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 <math.h>
+#include <stdlib.h>
+#include <stddef.h>
+#include <stdio.h>
+
+#include "error.h"
+#include "part.h"
+#include "cell.h"
+#include "tools.h"
+#include "swift.h"
+
+
+/**
+ *  Factorize a given integer, attempts to keep larger pair of factors.
+ */
+void factor(int value, int *f1, int *f2) {
+  int j;
+  int i;
+
+  j = (int)sqrt(value);
+  for (i = j; i > 0; i--) {
+    if ((value % i) == 0) {
+      *f1 = i;
+      *f2 = value / i;
+      break;
+    }
+  }
+}
+
+
+
+/**
+ * @brief Compute the average number of pairs per particle using
+ *      a brute-force O(N^2) computation.
+ *
+ * @param dim The space dimensions.
+ * @param parts The #part array.
+ * @param N The number of parts.
+ * @param periodic Periodic boundary conditions flag.
+ */
+
+void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
+              int periodic) {
+
+  int i, j, k, count = 0;
+  // int mj, mk;
+  // double maxratio = 1.0;
+  double r2, dx[3], rho = 0.0;
+  double rho_max = 0.0, rho_min = 100;
+
+  /* Loop over all particle pairs. */
+  for (j = 0; j < N; j++) {
+    if (j % 1000 == 0) {
+      printf("pairs_n2: j=%i.\n", j);
+      fflush(stdout);
+    }
+    for (k = j + 1; k < N; k++) {
+      for (i = 0; i < 3; i++) {
+        dx[i] = parts[j].x[i] - parts[k].x[i];
+        if (periodic) {
+          if (dx[i] < -dim[i] / 2)
+            dx[i] += dim[i];
+          else if (dx[i] > dim[i] / 2)
+            dx[i] -= dim[i];
+        }
+      }
+      r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+      if (r2 < parts[j].h * parts[j].h || r2 < parts[k].h * parts[k].h) {
+        runner_iact_density(r2, NULL, parts[j].h, parts[k].h, &parts[j],
+                            &parts[k]);
+        /* if ( parts[j].h / parts[k].h > maxratio )
+            {
+            maxratio = parts[j].h / parts[k].h;
+            mj = j; mk = k;
+            }
+        else if ( parts[k].h / parts[j].h > maxratio )
+            {
+            maxratio = parts[k].h / parts[j].h;
+            mj = j; mk = k;
+            } */
+      }
+    }
+  }
+
+  /* Aggregate the results. */
+  for (k = 0; k < N; k++) {
+    // count += parts[k].icount;
+    rho += parts[k].density.wcount;
+    rho_min = fmin(parts[k].density.wcount, rho_min);
+    rho_min = fmax(parts[k].density.wcount, rho_max);
+  }
+
+  /* Dump the result. */
+  printf("pairs_n2: avg. density per part is %.3f (nr. pairs %.3f).\n",
+         rho / N + 32.0 / 3, ((double)count) / N);
+  printf("pairs_n2: densities are in [ %e , %e ].\n", rho_min / N + 32.0 / 3,
+         rho_max / N + 32.0 / 3);
+  /* printf( "pairs_n2: maximum ratio between parts %i [%e,%e,%e] and %i
+     [%e,%e,%e] is %.3f/%.3f\n" ,
+      mj , parts[mj].x[0] , parts[mj].x[1] , parts[mj].x[2] ,
+      mk , parts[mk].x[0] , parts[mk].x[1] , parts[mk].x[2] ,
+      parts[mj].h , parts[mk].h ); fflush(stdout); */
+  fflush(stdout);
+}
+
+void pairs_single_density(double *dim, long long int pid,
+                          struct part *__restrict__ parts, int N,
+                          int periodic) {
+
+  int i, k;
+  // int mj, mk;
+  // double maxratio = 1.0;
+  double r2, dx[3];
+  float fdx[3];
+  struct part p;
+  // double ih = 12.0/6.25;
+
+  /* Find "our" part. */
+  for (k = 0; k < N && parts[k].id != pid; k++)
+    ;
+  if (k == N) error("Part not found.");
+  p = parts[k];
+  printf("pairs_single: part[%i].id == %lli.\n", k, pid);
+
+  p.rho = 0.0;
+  p.density.wcount = 0.0;
+  // p.icount = 0;
+  p.rho_dh = 0.0;
+
+  /* Loop over all particle pairs. */
+  for (k = 0; k < N; k++) {
+    if (parts[k].id == p.id) continue;
+    for (i = 0; i < 3; i++) {
+      dx[i] = p.x[i] - parts[k].x[i];
+      if (periodic) {
+        if (dx[i] < -dim[i] / 2)
+          dx[i] += dim[i];
+        else if (dx[i] > dim[i] / 2)
+          dx[i] -= dim[i];
+      }
+      fdx[i] = dx[i];
+    }
+    r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
+    if (r2 < p.h * p.h) {
+      runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
+      /* printf( "pairs_simple: interacting particles %lli [%i,%i,%i] and %lli
+         [%i,%i,%i], r=%e.\n" ,
+          pid , (int)(p.x[0]*ih) , (int)(p.x[1]*ih) , (int)(p.x[2]*ih) ,
+          parts[k].id , (int)(parts[k].x[0]*ih) , (int)(parts[k].x[1]*ih) ,
+         (int)(parts[k].x[2]*ih) ,
+          sqrtf(r2) ); */
+    }
+  }
+
+  /* Dump the result. */
+  printf("pairs_single: wcount of part %lli (h=%e) is %f.\n", p.id, p.h,
+         p.density.wcount + 32.0 / 3);
+  fflush(stdout);
+}
+
+void pairs_single_grav(double *dim, long long int pid,
+                       struct gpart *__restrict__ parts, int N, int periodic) {
+
+  int i, k;
+  // int mj, mk;
+  // double maxratio = 1.0;
+  double r2, dx[3];
+  float fdx[3], a[3] = {0.0, 0.0, 0.0}, aabs[3] = {0.0, 0.0, 0.0};
+  struct gpart pi, pj;
+  // double ih = 12.0/6.25;
+
+  /* Find "our" part. */
+  for (k = 0; k < N; k++)
+    if ((parts[k].id > 0 && parts[k].part->id == pid) || parts[k].id == -pid)
+      break;
+  if (k == N) error("Part not found.");
+  pi = parts[k];
+  pi.a[0] = 0.0f;
+  pi.a[1] = 0.0f;
+  pi.a[2] = 0.0f;
+
+  /* Loop over all particle pairs. */
+  for (k = 0; k < N; k++) {
+    if (parts[k].id == pi.id) continue;
+    pj = parts[k];
+    for (i = 0; i < 3; i++) {
+      dx[i] = pi.x[i] - pj.x[i];
+      if (periodic) {
+        if (dx[i] < -dim[i] / 2)
+          dx[i] += dim[i];
+        else if (dx[i] > dim[i] / 2)
+          dx[i] -= dim[i];
+      }
+      fdx[i] = dx[i];
+    }
+    r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
+    runner_iact_grav(r2, fdx, &pi, &pj);
+    a[0] += pi.a[0];
+    a[1] += pi.a[1];
+    a[2] += pi.a[2];
+    aabs[0] += fabsf(pi.a[0]);
+    aabs[1] += fabsf(pi.a[1]);
+    aabs[2] += fabsf(pi.a[2]);
+    pi.a[0] = 0.0f;
+    pi.a[1] = 0.0f;
+    pi.a[2] = 0.0f;
+  }
+
+  /* Dump the result. */
+  message(
+      "acceleration on gpart %lli is a=[ %e %e %e ], |a|=[ %.2e %.2e %.2e ].\n",
+      pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
+}
+
+
+/**
+ * @brief Test the density function by dumping it for two random parts.
+ *
+ * @param N number of intervals in [0,1].
+ */
+
+void density_dump(int N) {
+
+  int k;
+  float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
+  struct part *pi[4], *pj[4], Pi[4], Pj[4];
+
+  /* Init the interaction parameters. */
+  for (k = 0; k < 4; k++) {
+    Pi[k].mass = 1.0f;
+    Pi[k].rho = 0.0f;
+    Pi[k].density.wcount = 0.0f;
+    Pj[k].mass = 1.0f;
+    Pj[k].rho = 0.0f;
+    Pj[k].density.wcount = 0.0f;
+    hi[k] = 1.0;
+    hj[k] = 1.0;
+    pi[k] = &Pi[k];
+    pj[k] = &Pj[k];
+  }
+
+  for (k = 0; k <= N; k++) {
+    r2[3] = r2[2];
+    r2[2] = r2[1];
+    r2[1] = r2[0];
+    r2[0] = ((float)k) / N;
+    Pi[0].density.wcount = 0;
+    Pj[0].density.wcount = 0;
+    runner_iact_density(r2[0], NULL, hi[0], hj[0], &Pi[0], &Pj[0]);
+    printf(" %e %e %e", r2[0], Pi[0].density.wcount, Pj[0].density.wcount);
+    Pi[0].density.wcount = 0;
+    Pj[0].density.wcount = 0;
+    Pi[1].density.wcount = 0;
+    Pj[1].density.wcount = 0;
+    Pi[2].density.wcount = 0;
+    Pj[2].density.wcount = 0;
+    Pi[3].density.wcount = 0;
+    Pj[3].density.wcount = 0;
+    runner_iact_vec_density(r2, NULL, hi, hj, pi, pj);
+    printf(" %e %e %e %e\n", Pi[0].density.wcount, Pi[1].density.wcount,
+           Pi[2].density.wcount, Pi[3].density.wcount);
+  }
+}
diff --git a/src/tools.h b/src/tools.h
new file mode 100644
index 0000000000000000000000000000000000000000..7354a51030fef3a1fc8445afce770f56d5866fd6
--- /dev/null
+++ b/src/tools.h
@@ -0,0 +1,32 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Copyright (c) 2015 Peter W. Draper (p.w.draper@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/>.
+ *
+ ******************************************************************************/
+
+void factor(int value, int *f1, int *f2);
+void density_dump(int N);
+void pairs_single_grav(double *dim, long long int pid,
+                       struct gpart *__restrict__ parts, int N, int periodic);
+void pairs_single_density(double *dim, long long int pid,
+                          struct part *__restrict__ parts, int N,
+                          int periodic);
+
+void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
+              int periodic);
+
diff --git a/src/units.c b/src/units.c
index e802d73355c9f387ec5f40ba71e46d564a15ff25..8c9fd14452e9e1fdfe029ac89d22d7cd43aa0ef7 100644
--- a/src/units.c
+++ b/src/units.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
  *
  * This program is free software: you can redistribute it and/or modify
@@ -300,7 +300,7 @@ float aFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
 }
 
 /**
- * @brief Returns a string containg the exponants of the base units making up
+ * @brief Returns a string containing the exponents of the base units making up
  * the conversion factors
  */
 void conversionString(char* buffer, struct UnitSystem* us,
@@ -316,7 +316,7 @@ void conversionString(char* buffer, struct UnitSystem* us,
  * @brief Returns the conversion factor for a given unit (expressed in terms of
  * the 5 fundamental units) in the chosen unit system
  * @param us The unit system used
- * @param baseUnitsExponants The exponant of each base units required to form
+ * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
 double generalConversionFactor(struct UnitSystem* us,
@@ -334,7 +334,7 @@ double generalConversionFactor(struct UnitSystem* us,
  * @brief Returns the h factor exponentiation for a given unit (expressed in
  * terms of the 5 fundamental units)
  * @param us The unit system used
- * @param baseUnitsExponants The exponant of each base units required to form
+ * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
 float generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
@@ -351,7 +351,7 @@ float generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
  * @brief Returns the scaling factor exponentiation for a given unit (expressed
  * in terms of the 5 fundamental units)
  * @param us The unit system used
- * @param baseUnitsExponants The exponant of each base units required to form
+ * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
 float generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
@@ -363,12 +363,12 @@ float generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
 }
 
 /**
- * @brief Returns a string containg the exponants of the base units making up
+ * @brief Returns a string containing the exponents of the base units making up
  * the conversion factors (expressed in terms of the 5 fundamental units)
  * @param buffer The buffer in which to write (The buffer must be long enough,
  * 140 chars at most)
- * @param us The UnistSystem in use.
- * @param baseUnitsExponants The exponant of each base units required to form
+ * @param us The UnitsSystem in use.
+ * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
 void generalConversionString(char* buffer, struct UnitSystem* us,
diff --git a/src/units.h b/src/units.h
index 0523247b04f1d09c9974bae5afab8070590d77b1..1b977529784c1ef3069e1e932b16fd0b87073786 100644
--- a/src/units.h
+++ b/src/units.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2012 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
@@ -143,14 +143,14 @@ float generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
 float aFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
 
 /**
- * @brief Returns a string containg the exponants of the base units making up
+ * @brief Returns a string containing the exponents of the base units making up
  * the conversion factors (expressed in terms of the 5 fundamental units)
  */
 void generalConversionString(char* buffer, struct UnitSystem* us,
                              float baseUnitsExponants[5]);
 
 /**
- * @brief Returns a string containg the exponants of the base units making up
+ * @brief Returns a string containing the exponents of the base units making up
  * the conversion factors
  */
 void conversionString(char* buffer, struct UnitSystem* us,
diff --git a/src/vector.h b/src/vector.h
index cddc933ab1c64e537ce68fcb0a77e74bbe6818d3..8629b39725ebb7d800047116f45fb7ec6c6e7a1d 100644
--- a/src/vector.h
+++ b/src/vector.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
  *               2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
diff --git a/tests/Makefile.am b/tests/Makefile.am
index dbdd20d7a4c26daf05c6bd20a4b47b6491a72dc4..9ab592e6ccde6194ffcc85c510554d195e7da9c4 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,5 +1,5 @@
 # This file is part of SWIFT.
-# Coypright (c) 2015 matthieu.schaller@durham.ac.uk.
+# Copyright (c) 2015 matthieu.schaller@durham.ac.uk.
 # 
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -20,11 +20,10 @@ AM_CFLAGS = -I../src -DCPU_TPS=2.67e9 $(HDF5_CPPFLAGS)
 AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
 
 # List of programs and scripts to run in the test suite
-TESTS = testGreetings testReading.sh testTimeIntegration
-
+TESTS = testGreetings testReading.sh testSingle testTimeIntegration
 
 # List of test programs to compile
-check_PROGRAMS = testGreetings testReading testTimeIntegration
+check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration
 
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
@@ -32,3 +31,8 @@ testGreetings_SOURCES = testGreetings.c
 testReading_SOURCES = testReading.c
 
 testTimeIntegration_SOURCES = testTimeIntegration.c
+
+# Sources for test_single
+testSingle_SOURCES = testSingle.c
+testSingle_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
+testSingle_LDADD =  ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
diff --git a/tests/makeInput.py b/tests/makeInput.py
index 3cb802e381855ee27fd4ec520276b51a900d634a..8d8c07faf1e528869616bf93fbb7ffe8bf196af0 100644
--- a/tests/makeInput.py
+++ b/tests/makeInput.py
@@ -1,6 +1,6 @@
 ###############################################################################
  # This file is part of SWIFT.
- # Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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
diff --git a/examples/test_single.c b/tests/testSingle.c
similarity index 98%
rename from examples/test_single.c
rename to tests/testSingle.c
index 36c219eebf214feaf775283834578c5a487d5e90..02c52160d0d442496629f1bb3947f89524964fb8 100644
--- a/examples/test_single.c
+++ b/tests/testSingle.c
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  * 
  * This program is free software: you can redistribute it and/or modify