diff --git a/examples/SedovBlast/rdf.py b/examples/SedovBlast/rdf.py
deleted file mode 100644
index a8474939e2732f1b5291fc088f4a19f51d52ce5c..0000000000000000000000000000000000000000
--- a/examples/SedovBlast/rdf.py
+++ /dev/null
@@ -1,120 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Coypright (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/>.
- # 
- ##############################################################################
-
-import h5py
-import random
-import sys
-import math
-from numpy import *
-
-# Reads the HDF5 output of SWIFT and generates a radial density profile
-# of the different physical values.
-
-# Input values?
-if len(sys.argv) < 3 :
-    print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
-    exit()
-    
-# Get the input arguments
-fileName = sys.argv[1];
-nr_bins = int( sys.argv[2] );
-
-
-# Open the file
-fileName = sys.argv[1];
-file = h5py.File( fileName , 'r' )
-
-# Get the space dimensions.
-grp = file[ "/Header" ]
-boxsize = grp.attrs[ 'BoxSize' ]
-boxsize = boxsize[0]
-
-# Get the particle data
-grp = file.get( '/PartType0' )
-ds = grp.get( 'Coordinates' )
-coords = ds[...]
-ds = grp.get( 'Velocities' )
-v = ds[...]
-ds = grp.get( 'Masses' )
-m = ds[...]
-ds = grp.get( 'SmoothingLength' )
-h = ds[...]
-ds = grp.get( 'InternalEnergy' )
-u = ds[...]
-ds = grp.get( 'ParticleIDs' )
-ids = ds[...]
-ds = grp.get( 'Density' )
-rho = ds[...]
-
-# Set the origin to be the middle of the box
-origin = [ boxsize / 2 , boxsize / 2 , boxsize / 2 ]
-
-# Get the maximum radius
-r_max = boxsize / 2
-
-# Init the bins
-nr_parts = coords.shape[0]
-bins_v = zeros( nr_bins )
-bins_m = zeros( nr_bins )
-bins_h = zeros( nr_bins )
-bins_u = zeros( nr_bins )
-bins_rho = zeros( nr_bins )
-bins_count = zeros( nr_bins )
-bins_P = zeros( nr_bins )
-
-# Loop over the particles and fill the bins.
-for i in range( nr_parts ):
-
-    # Get the box index.
-    r = coords[i] - origin
-    r = sqrt( r[0]*r[0] + r[1]*r[1] + r[2]*r[2] )
-    ind = floor( r / r_max * nr_bins )
-    
-    # Update the bins
-    if ( ind < nr_bins ):
-        bins_count[ind] += 1
-        bins_v[ind] += sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
-        bins_m[ind] += m[i]
-        bins_h[ind] += h[i]
-        bins_u[ind] += u[i]
-        bins_rho[ind] += rho[i]
-        bins_P[ind] += (2.0/3)*u[i]*rho[i]
-    
-# Loop over the bins and dump them
-print "# bucket left right count v m h u rho"
-for i in range( nr_bins ):
-
-    # Normalize by the bin volume.
-    r_left = r_max * i / nr_bins
-    r_right = r_max * (i+1) / nr_bins
-    vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
-    ivol = 1.0 / vol
-
-    print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
-        ( i , r_left , r_right , \
-          bins_count[i] * ivol , \
-          bins_v[i] / bins_count[i] , \
-          bins_m[i] * ivol , \
-          bins_h[i] / bins_count[i] , \
-          bins_u[i] / bins_count[i] , \
-          bins_rho[i] / bins_count[i] ,
-          bins_P[i] / bins_count[i] )
-    
-    
diff --git a/examples/SedovBlast/solution.py b/examples/SedovBlast/solution.py
index 1b2dcf5957bc6684ebaa88d4ea08bee68e56cd6f..cd26aa4d408de58ac71c570cd131a017f13ec792 100644
--- a/examples/SedovBlast/solution.py
+++ b/examples/SedovBlast/solution.py
@@ -33,7 +33,7 @@ P_0 = 1.e-5         # Background Pressure
 E_0 = 1.e2          # Energy of the explosion
 gamma = 5./3.     # Gas polytropic index
 
-t = 0.2           # Time of the solution
+t = 0.15           # Time of the solution
 
 N = 1000          # Number of radial points
 R_max = 3.        # Maximal radius
diff --git a/examples/SodShock/makeIC.py b/examples/SodShock/makeIC.py
index 8278862d4927e49f0e972bedc974b7fd76237d87..bd7bef3d4d18baa271d68e258d5a88161f56eb4f 100644
--- a/examples/SodShock/makeIC.py
+++ b/examples/SodShock/makeIC.py
@@ -37,36 +37,21 @@ fileName = "sodShock.hdf5"
 #---------------------------------------------------
 
 #Read in high density glass
-glass1 = h5py.File("../Glass/glass_200000.hdf5")
-pos1 = glass1["/PartType1/Coordinates"][:,:]
-pos1 = pos1 / 20. # Particles are in [0:0.5, 0:0.5, 0:0.5]
-
-toRemove = array(0)
-for i in range(size(pos1)/3-1):
-    if pos1[i,0] == pos1[i+1,0]:
-        if pos1[i,1] == pos1[i+1,1]:
-            if pos1[i,2] == pos1[i+1,2]:
-                toRemove = append(toRemove, i)
-
-pos1 = delete(pos1, toRemove, 0)
+# glass1 = h5py.File("../Glass/glass_200000.hdf5")
+glass1 = h5py.File("glass_001.hdf5")
+pos1 = glass1["/PartType0/Coordinates"][:,:]
+pos1 = pos1 / 2. # Particles are in [0:0.5, 0:0.5, 0:0.5]
 
 
 #Read in high density glass
-glass2 = h5py.File("../Glass/glass_50000.hdf5")
-pos2 = glass2["/PartType1/Coordinates"][:,:]
-pos2 = pos2 / 20. # Particles are in [0:0.5, 0:0.5, 0:0.5]
-
-toRemove = array(0)
-for i in range(size(pos2)/3-1):
-    if pos2[i,0] == pos2[i+1,0]:
-        if pos2[i,1] == pos2[i+1,1]:
-            if pos2[i,2] == pos2[i+1,2]:
-                toRemove = append(toRemove, i)
-pos2 = delete(pos2, toRemove, 0)
+# glass2 = h5py.File("../Glass/glass_50000.hdf5")
+glass2 = h5py.File("glass_002.hdf5")
+pos2 = glass2["/PartType0/Coordinates"][:,:]
+pos2 = pos2 / 2. # Particles are in [0:0.5, 0:0.5, 0:0.5]
 
 
 #Generate high density region
-rho1 = size(pos1) / size(pos2) 
+rho1 = 1.
 coord1 = append(pos1, pos1 + [0, 0.5, 0], 0)
 coord1 = append(coord1, pos1 + [0, 0, 0.5], 0)
 coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
@@ -77,7 +62,7 @@ u1 = ones((N1, 1)) * P1 / ((gamma - 1.) * rho1)
 m1 = ones((N1, 1)) * 0.5 * rho1 / N1
 
 #Generate low density region
-rho2 = 1.
+rho2 = 0.25
 coord2 = append(pos2, pos2 + [0, 0.5, 0], 0)
 coord2 = append(coord2, pos2 + [0, 0, 0.5], 0)
 coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
@@ -133,4 +118,7 @@ ds[()] = ids[:]
 
 file.close()
 
+print "particle 31075: " , coords[31075,:]
+print "particle 31076: " , coords[31076,:]
+
 
diff --git a/examples/SodShock/rhox.py b/examples/SodShock/rhox.py
deleted file mode 100644
index 36771c5c0d9dcdec99403249ea8a7bb1a4bc58c6..0000000000000000000000000000000000000000
--- a/examples/SodShock/rhox.py
+++ /dev/null
@@ -1,115 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Coypright (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/>.
- # 
- ##############################################################################
-
-import h5py
-import random
-import sys
-import math
-from numpy import *
-
-# Reads the HDF5 output of SWIFT and generates a radial density profile
-# of the different physical values.
-
-# Input values?
-if len(sys.argv) < 3 :
-    print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
-    exit()
-    
-# Get the input arguments
-fileName = sys.argv[1];
-nr_bins = int( sys.argv[2] );
-
-
-# Open the file
-fileName = sys.argv[1];
-file = h5py.File( fileName , 'r' )
-
-# Get the space dimensions.
-grp = file[ "/Header" ]
-boxsize = grp.attrs[ 'BoxSize' ]
-boxsize = boxsize[0]
-
-# Get the particle data
-grp = file.get( '/PartType0' )
-ds = grp.get( 'Coordinates' )
-coords = ds[...]
-ds = grp.get( 'Velocities' )
-v = ds[...]
-# ds = grp.get( 'Mass' )
-# m = ds[...]
-ds = grp.get( 'SmoothingLength' )
-h = ds[...]
-ds = grp.get( 'InternalEnergy' )
-u = ds[...]
-ds = grp.get( 'ParticleIDs' )
-ids = ds[...]
-ds = grp.get( 'Density' )
-rho = ds[...]
-
-# Get the maximum radius
-r_max = boxsize
-
-# Init the bins
-nr_parts = coords.shape[0]
-bins_v = zeros( nr_bins )
-bins_m = zeros( nr_bins )
-bins_h = zeros( nr_bins )
-bins_u = zeros( nr_bins )
-bins_rho = zeros( nr_bins )
-bins_count = zeros( nr_bins )
-bins_P = zeros( nr_bins )
-
-# Loop over the particles and fill the bins.
-for i in range( nr_parts ):
-
-    # Get the box index.
-    r = coords[i,0]
-    ind = floor( r / r_max * nr_bins )
-    
-    # Update the bins
-    bins_count[ind] += 1
-    bins_v[ind] += v[i,0] # sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
-    # bins_m[ind] += m[i]
-    bins_h[ind] += h[i]
-    bins_u[ind] += u[i]
-    bins_rho[ind] += rho[i]
-    bins_P[ind] += (2.0/3)*u[i]*rho[i]
-    
-# Loop over the bins and dump them
-print "# bucket left right count v m h u rho"
-for i in range( nr_bins ):
-
-    # Normalize by the bin volume.
-    r_left = r_max * i / nr_bins
-    r_right = r_max * (i+1) / nr_bins
-    vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
-    ivol = 1.0 / vol
-
-    print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
-        ( i , r_left , r_right , \
-          bins_count[i] * ivol , \
-          bins_v[i] / bins_count[i] , \
-          bins_m[i] * ivol , \
-          bins_h[i] / bins_count[i] , \
-          bins_u[i] / bins_count[i] , \
-          bins_rho[i] / bins_count[i] ,
-          bins_P[i] / bins_count[i] )
-    
-    
diff --git a/examples/SodShock/solution.py b/examples/SodShock/solution.py
index 56f280e9fb4c032fb6d888668af42d52fdba2d7f..a5a1364d90d4af12f3ac051e23748d6228592258 100644
--- a/examples/SodShock/solution.py
+++ b/examples/SodShock/solution.py
@@ -29,11 +29,11 @@ from numpy import *
 
 
 # Parameters
-rho_L = 4.
-P_L = 1.
+rho_L = 1
+P_L = 1
 v_L = 0.
 
-rho_R = 1.
+rho_R = 0.25
 P_R = 0.1795
 v_R = 0.
 
@@ -42,8 +42,8 @@ gamma = 5./3. # Polytropic index
 t = 0.12  # Time of the evolution
 
 N = 1000  # Number of points
-x_min = -0.15
-x_max = 0.15
+x_min = -0.25
+x_max = 0.25
 
 
 # ---------------------------------------------------------------
diff --git a/examples/test.c b/examples/test.c
index 7155f233aaa4bde4ea1e5624bfb731e51d478ad3..d96b901964caaececbe5a7fc194969f55a2a9496 100644
--- a/examples/test.c
+++ b/examples/test.c
@@ -885,7 +885,7 @@ int main ( int argc , char *argv[] ) {
         engine_step( &e , 0 );
 
         if(j % 100 == 0)
-	  write_output(&e);
+            write_output(&e);
 	
         
         /* Dump the first few particles. */
@@ -913,7 +913,7 @@ int main ( int argc , char *argv[] ) {
         for ( k = 0 ; k < s.nr_parts ; k++ )
             if ( s.parts[k].dt < p->dt )
                 p = &s.parts[k];
-        printf( "main: particle %lli/%i at [ %e %e %e ] has smallest dt=%.3e (h=%.3e,u=%.3e).\n" ,
+        printf( "main: particle %lli/%i at [ %e %e %e ] has smallest dt=%e (h=%.3e,u=%.3e).\n" ,
 	        p->id , (int)(p - s.parts) , p->x[0] , p->x[1] , p->x[2] , p->dt , p->h , p->u ); */
         
         /* Output. */
@@ -982,6 +982,9 @@ int main ( int argc , char *argv[] ) {
     printf( "main: particle %lli/%i at [ %e %e %e ] (h=%e) has maximum wcount %.3f.\n" ,
 	    p->id , (int)(p - s.parts) , p->x[0] , p->x[1] , p->x[2] , p->h , p->density.wcount );
     
+    /* Write final output. */
+    write_output( &e );
+    
     /* Get the average interactions per particle. */
     // icount = 0;
     // space_map_parts( &s , &map_icount , &icount );
diff --git a/examples/test_single.c b/examples/test_single.c
index 3d1657e1fe481992d5cf9a6b24bfdda383bf4d85..f2aae8095d40d4ec5cb8e2617eb688b2039efa67 100644
--- a/examples/test_single.c
+++ b/examples/test_single.c
@@ -64,23 +64,28 @@ int main ( int argc , char *argv[] ) {
 
     int k, N = 100;
     struct part p1, p2;
-    float x, w, dwdx, r2, dx[3] = { 0.0f , 0.0f , 0.0f };
+    float x, w, dwdx, r2, dx[3] = { 0.0f , 0.0f , 0.0f }, gradw[3];
     
     /* Init the particles. */
     for ( k = 0 ; k < 3 ; k++ ) {
         p1.a[k] = 0.0f; p1.v[k] = 0.0f; p1.x[k] = 0.0;
         p2.a[k] = 0.0f; p2.v[k] = 0.0f; p2.x[k] = 0.0;
         }
+    p1.v[0] = 100.0f;
+    p1.id = 0; p2.id = 1;
+    p1.density.wcount = 48.0f; p2.density.wcount = 48.0f;
     p1.rho = 1.0f; p1.mass = 9.7059e-4; p1.h = 0.222871287 / 2;
     p2.rho = 1.0f; p2.mass = 9.7059e-4; p2.h = 0.222871287 / 2;
-    p1.force.c = 0.0f; p1.force.balsara = 0.0f;
-    p2.force.c = 0.0f; p2.force.balsara = 0.0f;
+    p1.force.c = 0.0040824829f; p1.force.balsara = 0.0f;
+    p2.force.c = 58.8972740361f; p2.force.balsara = 0.0f;
     p1.u = 1.e-5 / ((const_gamma - 1.)*p1.rho);
     p2.u = 1.e-5 / ((const_gamma - 1.)*p2.rho) + 100.0f / ( 33 * p2.mass );
     p1.force.POrho2 = p1.u * ( const_gamma - 1.0f ) / p1.rho;
     p2.force.POrho2 = p2.u * ( const_gamma - 1.0f ) / p2.rho;
     
     /* Dump a header. */
+    printParticle_single( &p1 );
+    printParticle_single( &p2 );
     printf( "# r a_1 udt_1 a_2 udt_2\n" );
     
     /* Loop over the different radii. */
@@ -91,28 +96,31 @@ int main ( int argc , char *argv[] ) {
         r2 = dx[0]*dx[0];
         
         /* Clear the particle fields. */
-        /* p1.a[0] = 0.0f; p1.force.u_dt = 0.0f;
-        p2.a[0] = 0.0f; p2.force.u_dt = 0.0f; */
+        p1.a[0] = 0.0f; p1.force.u_dt = 0.0f;
+        p2.a[0] = 0.0f; p2.force.u_dt = 0.0f;
         
         /* Interact the particles. */
-        // runner_iact_force( r2 , dx , p1.h , p2.h , &p1 , &p2 );
+        runner_iact_force( r2 , dx , p1.h , p2.h , &p1 , &p2 );
         
         /* Clear the particle fields. */
-        p1.rho = 0.0f; p1.density.wcount = 0.0f;
-        p2.rho = 0.0f; p2.density.wcount = 0.0f;
+        /* p1.rho = 0.0f; p1.density.wcount = 0.0f;
+        p2.rho = 0.0f; p2.density.wcount = 0.0f; */
         
         /* Interact the particles. */
-        runner_iact_density( r2 , dx , p1.h , p2.h , &p1 , &p2 );
+        // runner_iact_density( r2 , dx , p1.h , p2.h , &p1 , &p2 );
         
         /* Evaluate just the kernel. */
         x = fabsf( dx[0] ) / p1.h;
         kernel_deval( x , &w , &dwdx );
+        gradw[0] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[0] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
+        gradw[1] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[1] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
+        gradw[2] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[2] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
         
         /* Output the results. */
-        printf( "%.3e %.3e %.3e %.3e %.3e %.3e %.3e\n" ,
-            // -dx[0] , p1.a[0] , p1.force.u_dt , p2.a[0] , p2.force.u_dt ,
-            -dx[0] , p1.rho , p1.density.wcount , p2.rho , p2.density.wcount ,
-            w , dwdx );
+        printf( "%.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e\n" ,
+            -dx[0] , p1.a[0] , p1.a[1] , p1.a[2] , p1.force.u_dt , 
+            /// -dx[0] , p1.rho , p1.density.wcount , p2.rho , p2.density.wcount ,
+            w , dwdx , gradw[0] , gradw[1] , gradw[2] );
     
         } /* loop over radii. */