diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index 566ec8e94ba2ddf56607a771d760192e721d17db..c5a0fc1d7ba29d91a33b6ba9616ef13be1677061 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -36,113 +36,117 @@ import os
 
 print "Plotting..."
 
-# Read Quickshed accelerations
-data=loadtxt("interaction_dump.dat")
-id = data[:,0]
-posx = data[:,2]
-
-accx_u=data[:,5]
-accy_u=data[:,6]
-accz_u=data[:,7]
-
-accx_s=data[:,8]
-accy_s=data[:,9]
-accz_s=data[:,10]
-
-
-
-
-# Build error ------------------------------------------------
-
-errx_s = (accx_s - accx_u )/abs(accx_u) 
-erry_s = (accy_s - accy_u )/abs(accy_u) 
-errz_s = (accz_s - accz_u )/abs(accz_u) 
-
-e_errx_s = errx_s#[abs(errx_s) > 0.001]
-e_erry_s = erry_s#[abs(erry_s) > 0.001]
-e_errz_s = errz_s#[abs(errz_s) > 0.001]
-
-# Statistics
-meanx_s = mean(errx_s[abs(errx_s) < 0.1])
-stdx_s = std(errx_s[abs(errx_s) < 0.1])
-meany_s = mean(erry_s[abs(erry_s) < 0.1])
-stdy_s = std(erry_s[abs(erry_s) < 0.1])
-meanz_s = mean(errz_s[abs(errz_s) < 0.1])
-stdz_s = std(errz_s[abs(errz_s) < 0.1])
-
-
-
-# Plot -------------------------------------------------------
-figure(frameon=True)
-
-subplot(311, title="Acceleration along X")
-#plot(id[abs(errx_s) > 0.001], e_errx_s , 'rx')
-plot(posx, e_errx_s , 'rx')
-#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" )
-ylim(-0.2, 0.2)
-#xlim(0, size(id)-1)
-grid()
-
-subplot(312, title="Acceleration along Y")
-#plot(id[abs(erry_s) > 0.001], e_erry_s , 'rx')
-plot(posx, e_erry_s , 'rx')
-#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" )
-ylim(-0.2, 0.2)
-#xlim(0, size(id)-1)
-
-grid()
-
-subplot(313, title="Acceleration along Z")
-#plot(id[abs(errz_s) > 0.001], e_errz_s , 'rx', label="Sorted")
-plot(posx, e_errz_s , 'rx', label="Sorted")
-#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanz_bh, stdz_bh, meanz_new, stdz_new), backgroundcolor="w", va="top", ha="right" )
-legend(loc="upper right")
-
-ylim(-0.2, 0.2)
-#xlim(0, size(id)-1)
-grid()
-
-savefig("accelerations.png")
-
-
-
-
-# Plot -------------------------------------------------------
-bins = linspace(-3, 3, 10000)
-
-e_errx_s = errx_s[(abs(errx_s) > 0.001) & (abs(errx_s) < 1.)]
-e_erry_s = erry_s[(abs(erry_s) > 0.001) & (abs(erry_s) < 1.)]
-e_errz_s = errz_s[(abs(errz_s) > 0.001) & (abs(errz_s) < 1.)]
-
-
-
-figure(frameon=True)
-subplot(311, title="Acceleration along X")#, yscale='log')
-hist(e_errx_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Sorted")
-legend(loc="upper right")
-xlim(-0.12, 0.12)
-
-subplot(312, title="Acceleration along Y")
-hist(e_erry_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
-xlim(-0.12, 0.12)
-
-subplot(313, title="Acceleration along Z")
-hist(e_errz_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
-xlim(-0.12, 0.12)
-
-savefig("histogram.png")
-
-
-
-
-
-print "E_error for sorted: ( x= %4.3e"%std(e_errx_s), "y= %4.3e"%std(e_erry_s), "z= %4.3e"%std(e_errz_s), ") Combined YZ = %4.3e"%(( std(e_erry_s) + std(e_errz_s) )/2.)
-
-#print std(e_errx_s), std(e_erry_s), std(e_errz_s)
-
-
-#print "Min   for sorted: ( x= %4.3e"%min(errx_s), "y= %4.3e"%min(erry_s), "z= %4.3e"%min(errz_s), ") Combined YZ = %4.3e"%(min( min(erry_s), min(errz_s) ))
-#print "Max   for sorted: ( x= %4.3e"%max(errx_s), "y= %4.3e"%max(erry_s), "z= %4.3e"%max(errz_s), ") Combined YZ = %4.3e"%(max( max(erry_s), max(errz_s) ))
+names = ["side", "edge", "corner"]
+
+for orientation in range( 3 ):
+
+    # Read Quickshed accelerations
+    data=loadtxt( "interaction_dump_%d.dat"%orientation )
+    id = data[:,0]
+    pos = data[:,2]
+    
+    accx_u=data[:,6]
+    accy_u=data[:,7]
+    accz_u=data[:,8]
+    
+    accx_s=data[:,9]
+    accy_s=data[:,10]
+    accz_s=data[:,11]
+    
+
+
+
+    # Build error ------------------------------------------------
+    
+    errx_s = (accx_s - accx_u )/abs(accx_u) 
+    erry_s = (accy_s - accy_u )/abs(accy_u) 
+    errz_s = (accz_s - accz_u )/abs(accz_u) 
+    
+    e_errx_s = errx_s#[abs(errx_s) > 0.001]
+    e_erry_s = erry_s#[abs(erry_s) > 0.001]
+    e_errz_s = errz_s#[abs(errz_s) > 0.001]
+    
+    # Statistics
+    meanx_s = mean(errx_s[abs(errx_s) < 0.1])
+    stdx_s = std(errx_s[abs(errx_s) < 0.1])
+    meany_s = mean(erry_s[abs(erry_s) < 0.1])
+    stdy_s = std(erry_s[abs(erry_s) < 0.1])
+    meanz_s = mean(errz_s[abs(errz_s) < 0.1])
+    stdz_s = std(errz_s[abs(errz_s) < 0.1])
+    
+
+
+    # Plot -------------------------------------------------------
+    figure(frameon=True)
+    
+    subplot(311, title="Acceleration along X")
+    #plot(id[abs(errx_s) > 0.001], e_errx_s , 'rx')
+    plot(pos, e_errx_s , 'rx')
+    #text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" )
+    ylim(-0.2, 0.2)
+    #xlim(0, size(id)-1)
+    grid()
+    
+    subplot(312, title="Acceleration along Y")
+    #plot(id[abs(erry_s) > 0.001], e_erry_s , 'rx')
+    plot(pos, e_erry_s , 'rx')
+    #text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" )
+    ylim(-0.2, 0.2)
+    #xlim(0, size(id)-1)
+    
+    grid()
+    
+    subplot(313, title="Acceleration along Z")
+    #plot(id[abs(errz_s) > 0.001], e_errz_s , 'rx', label="Sorted")
+    plot(pos, e_errz_s , 'rx', label="Sorted")
+    #text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanz_bh, stdz_bh, meanz_new, stdz_new), backgroundcolor="w", va="top", ha="right" )
+    legend(loc="upper right")
+    
+    ylim(-0.2, 0.2)
+    #xlim(0, size(id)-1)
+    grid()
+
+    savefig("accelerations_%s.png"% names[orientation])
+    
+    
+    
+    
+    # Plot -------------------------------------------------------
+    bins = linspace(-3, 3, 10000)
+    
+    e_errx_s = errx_s[(abs(errx_s) > 0.001) & (abs(errx_s) < 1.)]
+    e_erry_s = erry_s[(abs(erry_s) > 0.001) & (abs(erry_s) < 1.)]
+    e_errz_s = errz_s[(abs(errz_s) > 0.001) & (abs(errz_s) < 1.)]
+    
+    
+    
+    figure(frameon=True)
+    subplot(311, title="Acceleration along X")#, yscale='log')
+    hist(e_errx_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Sorted")
+    legend(loc="upper right")
+    xlim(-0.12, 0.12)
+    
+    subplot(312, title="Acceleration along Y")
+    hist(e_erry_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
+    xlim(-0.12, 0.12)
+
+    subplot(313, title="Acceleration along Z")
+    hist(e_errz_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
+    xlim(-0.12, 0.12)
+    
+    savefig("histogram_%s.png"% names[orientation])
+
+
+
+
+
+    print "E_error for sorted: ( x= %4.3e"%std(e_errx_s), "y= %4.3e"%std(e_erry_s), "z= %4.3e"%std(e_errz_s), ") Combined YZ = %4.3e"%(( std(e_erry_s) + std(e_errz_s) )/2.)
+
+    #print std(e_errx_s), std(e_erry_s), std(e_errz_s)
+
+
+    #print "Min   for sorted: ( x= %4.3e"%min(errx_s), "y= %4.3e"%min(erry_s), "z= %4.3e"%min(errz_s), ") Combined YZ = %4.3e"%(min( min(erry_s), min(errz_s) ))
+    #print "Max   for sorted: ( x= %4.3e"%max(errx_s), "y= %4.3e"%max(erry_s), "z= %4.3e"%max(errz_s), ") Combined YZ = %4.3e"%(max( max(erry_s), max(errz_s) ))
 
 
 
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index adb7cf7879deffba87d1924a31032d3f58ad9aee..57c36f2585b140dc3f6d6e577c36a561fb2851f9 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -44,7 +44,7 @@
 #define dist_cutoff_ratio 1.5
 #define iact_pair_direct iact_pair_direct_sorted
 
-#define ICHECK 2
+#define ICHECK -1
 #define NO_SANITY_CHECKS
 #define NO_COM_AS_TASK
 #define COUNTERS
@@ -2018,8 +2018,10 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
 /**
  * @brief Creates two neighbouring cells wiht N_parts per celland makes them interact using both the 
  * sorted and unsortde interactions. Outputs then the tow sets of accelerations for accuracy tests.
+ * @param N_parts Number of particles in each cell
+ * @param orientation Orientation of the cells ( 0 == side, 1 == edge, 2 == corner )
  */
-void test_direct_neighbour( int N_parts ) {
+void test_direct_neighbour( int N_parts , int orientation ) {
 
   int k;
   struct part *parts;
@@ -2032,19 +2034,34 @@ void test_direct_neighbour( int N_parts ) {
   /* Create random set of particles in both cells */
   for (k = 0; k < 2*N_parts; k++) {
       parts[k].id = k;
-      parts[k].x[0] = ((double)rand()) / RAND_MAX + (k >= N_parts ? 1 : 0. );
+
+      parts[k].x[0] = ((double)rand()) / RAND_MAX;
+      if ( k >= N_parts )
+	parts[k].x[0] += 1.;
+
       parts[k].x[1] = ((double)rand()) / RAND_MAX;
+      if ( orientation > 0 && k >= N_parts )
+      	parts[k].x[1] += 1.;
+
       parts[k].x[2] = ((double)rand()) / RAND_MAX;
-      parts[k].mass = 1.;//((double)rand()) / RAND_MAX;
+      if ( orientation > 1 && k >= N_parts )
+      	parts[k].x[2] += 1.;
+
+      parts[k].mass = ((double)rand()) / RAND_MAX;
       parts[k].a[0] = 0.0;
       parts[k].a[1] = 0.0;
       parts[k].a[2] = 0.0;
     }
 
   /* Get the cell geometry right */
-  left.loc[0] = 0.; left.loc[1] = 0.; left.loc[2] = 0.;
+  left.loc[0] = 0.; 
+  left.loc[1] = 0.; 
+  left.loc[2] = 0.;
   left.h = 1.;
-  right.loc[0] = 1.; right.loc[1] = 0.; right.loc[2] = 0.;
+  
+  right.loc[0] = 1.; 
+  right.loc[1] = ( orientation > 0 ? 1 : 0 ); 
+  right.loc[2] = ( orientation > 1 ? 1 : 0 );
   right.h = 1.;
 
   /* Put the particles in the cell */
@@ -2107,27 +2124,50 @@ void test_direct_neighbour( int N_parts ) {
     }
   }
 
-  /* Sort the particles along axis to simply interpretation */
-  int compParts(const void* c1, const void* c2) {
-    if      ( ((struct part*)c1)->x[0] <  ((struct part*)c2)->x[0] ) return -1;
-    else if ( ((struct part*)c1)->x[0] == ((struct part*)c2)->x[0] ) return 0;
-    else    return 1;
+  
+  /* Position along the axis */
+  double* position;
+  if ((position = (double *)malloc(sizeof(double) * N_parts * 2)) == NULL)
+    error("Failed to allocate position buffer.");
+
+  for ( k = 0; k < 2*N_parts ; ++k ) {
+   
+    switch( orientation ) {
+
+    case 0:
+      position[k] = parts[k].x[0] - 1.;
+      break;
+
+    case 1:
+      position[k] = sqrt( ( parts[k].x[0] * parts[k].x[0] ) + ( parts[k].x[1] * parts[k].x[1] ) ) - sqrt(2.);
+      break;
+
+    case 2:
+      position[k] = sqrt( ( parts[k].x[0] * parts[k].x[0] ) + ( parts[k].x[1] * parts[k].x[1] ) + ( parts[k].x[2] * parts[k].x[2] ) ) - sqrt(3.);
+      break;
+
+    default:
+      error("Wrong switch statement");
+      break;
+    }
   }
-  //qsort( parts, 2*N_parts, sizeof(struct part), compParts);
 
 
   /* Now, output everything */
-  message( "Writing file 'interaction_dump.dat'" );
-  FILE* file = fopen("interaction_dump.dat", "w");
-  fprintf(file, "# ID m x y z a_u.x   a_u.y    a_u.z    a_s.x    a_s.y    a_s.z\n");
+  char fileName[100];
+  sprintf( fileName, "interaction_dump_%d.dat", orientation );
+  message( "Writing file '%s'", fileName );
+  FILE* file = fopen( fileName , "w" );
+  fprintf(file, "# ID m r x y z a_u.x   a_u.y    a_u.z    a_s.x    a_s.y    a_s.z\n");
   for (k = 0; k < 2*N_parts; k++) {
-    fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e\n", parts[k].id, parts[k].mass,
+    fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id, parts[k].mass, position[k],
             parts[k].x[0], parts[k].x[1], parts[k].x[2], 
 	    parts[k].a_exact[0], parts[k].a_exact[1], parts[k].a_exact[2], 
 	    parts[k].a[0], parts[k].a[1], parts[k].a[2]);
   }
   fclose( file );
 
+
 #ifdef COUNTERS
   message( "Unsorted intereactions:           %d" ,count_direct_unsorted );
   message( "Sorted intereactions PP:          %d" ,count_direct_sorted_pp );
@@ -2211,7 +2251,9 @@ int main(int argc, char *argv[]) {
     message("Interacting 2 neighbouring cells with %d particles per cell", N_parts);
     
     /* Run the test */
-    test_direct_neighbour( N_parts );
+    test_direct_neighbour( N_parts, 0 );
+    test_direct_neighbour( N_parts, 1 );
+    test_direct_neighbour( N_parts, 2 );
     
   }
   else {