diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index c5a0fc1d7ba29d91a33b6ba9616ef13be1677061..8c47e73ff2a2d718625220ff00635cec979e10a6 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -22,7 +22,7 @@ params = {'axes.labelsize': 16,
           'figure.subplot.wspace'  : 0.26  ,  # the amount of width reserved for blank space between subplots                              
           'figure.subplot.hspace'  : 0.22  ,  # the amount of height reserved for white space between subplots                             
 
-          'lines.markersize' : 6,
+          'lines.markersize' : 2,
           'lines.linewidth' : 2,
 
 #          'axes.formatter.limits' : (-2, 1),                                                                                              
@@ -36,9 +36,26 @@ import os
 
 print "Plotting..."
 
-names = ["side", "edge", "corner"]
+axis = [
+    1.0,  1.0,  1.0, 
+    1.0,  1.0,  0.0, 
+    1.0,  1.0, -1.0, 
+    1.0,  0.0,  1.0, 
+    1.0,  0.0,  0.0, 
+    1.0,  0.0, -1.0, 
+    1.0, -1.0,  1.0, 
+    1.0, -1.0,  0.0, 
+    1.0, -1.0, -1.0, 
+    0.0,  1.0,  1.0, 
+    0.0,  1.0,  0.0, 
+    0.0,  1.0, -1.0, 
+    0.0,  0.0,  1.0  
+  ]
 
-for orientation in range( 3 ):
+
+#names = ["side", "edge", "corner"]
+
+for orientation in range( 13 ):
 
     # Read Quickshed accelerations
     data=loadtxt( "interaction_dump_%d.dat"%orientation )
@@ -80,34 +97,32 @@ for orientation in range( 3 ):
     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" )
+    #plot(id[abs(errx_s) > 0.001], e_errx_s , 'ro')
+    plot(pos, e_errx_s , 'ro')
+    text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
     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)
-    
+    #plot(id[abs(erry_s) > 0.001], e_erry_s , 'ro')
+    plot(pos, e_erry_s , 'ro')
+    text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+    ylim(-0.2, 0.2)  
     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" )
+    #plot(id[abs(errz_s) > 0.001], e_errz_s , 'ro', label="Sorted")
+    plot(pos, e_errz_s , 'ro', label="Sorted")
+    text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+
     legend(loc="upper right")
     
     ylim(-0.2, 0.2)
     #xlim(0, size(id)-1)
     grid()
 
-    savefig("accelerations_%s.png"% names[orientation])
-    
+    savefig("accelerations_%d.png"%orientation)
+    close()
     
     
     
@@ -134,8 +149,8 @@ for orientation in range( 3 ):
     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])
-
+    savefig("histogram_%d.png"%orientation)
+    close()
 
 
 
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 83c126f3801898c6682bc8540d5e1547f284842a..a9609ba8ea53165d8918abc684e2d092ed17c1be 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -169,6 +169,23 @@ const int axis_num_orth[13] = {
   /* (  0 , -1 ,  1 ) */ 1,
   /* (  0 ,  0 , -1 ) */ 2
 };
+/* const int axis_num_orth[13] = { */
+/*   /\* ( -1 , -1 , -1 ) *\/ 2, */
+/*   /\* ( -1 , -1 ,  0 ) *\/ 2, */
+/*   /\* ( -1 , -1 ,  1 ) *\/ 2, */
+/*   /\* ( -1 ,  0 , -1 ) *\/ 2, */
+/*   /\* ( -1 ,  0 ,  0 ) *\/ 2, */
+/*   /\* ( -1 ,  0 ,  1 ) *\/ 2, */
+/*   /\* ( -1 ,  1 , -1 ) *\/ 2, */
+/*   /\* ( -1 ,  1 ,  0 ) *\/ 2, */
+/*   /\* ( -1 ,  1 ,  1 ) *\/ 2, */
+/*   /\* (  0 , -1 , -1 ) *\/ 2, */
+/*   /\* (  0 , -1 ,  0 ) *\/ 2, */
+/*   /\* (  0 , -1 ,  1 ) *\/ 2, */
+/*   /\* (  0 ,  0 , -1 ) *\/ 2 */
+/* }; */
+
+
 const char axis_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0,
                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
@@ -465,8 +482,11 @@ void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i,
   orth1[3] = -d1;
   orth2[3] = -d2;
 
-  /* message("dx=[%.3e,%.3e,%.3e], axis=[%.3e,%.3e,%.3e], corr=%.3e.",
-    dx[0], dx[1], dx[2], axis[0], axis[1], axis[2], *corr); */
+  message("dx=[%.3e,%.3e,%.3e], axis=[%.3e,%.3e,%.3e]",
+	  dx[0], dx[1], dx[2], axis[0], axis[1], axis[2]);
+  message("N_planes= %d", *num_orth_planes);
+  message("o1=[%.3e,%.3e,%.3e] %.3e", orth1[0], orth1[1], orth1[2], orth1[3]);
+  message("o2=[%.3e,%.3e,%.3e] %.3e", orth2[0], orth2[1], orth2[2], orth2[3]);
 
   /* Make sure the sorts are ok. */
   /* for (int k = 1; k < (*ci)->count; k++)
@@ -2044,8 +2064,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
  * 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 )
+ * @param orientation Orientation of the cells ( 0 <= orientation < 13 )
  */
 void test_direct_neighbour(int N_parts, int orientation) {
 
@@ -2053,6 +2072,31 @@ void test_direct_neighbour(int N_parts, int orientation) {
   struct part *parts;
   struct cell left, right;
 
+  /* All 13 configurations */
+  const float cell_shift[13 * 3] = {
+    1.0,  1.0,  1.0, /* 0 */
+    1.0,  1.0,  0.0, /* 1 */
+    1.0,  1.0, -1.0, /* 2 */
+    1.0,  0.0,  1.0, /* 3 */
+    1.0,  0.0,  0.0, /* 4 */
+    1.0,  0.0, -1.0, /* 5 */
+    1.0, -1.0,  1.0, /* 6 */
+    1.0, -1.0,  0.0, /* 7 */
+    1.0, -1.0, -1.0, /* 8 */
+    0.0,  1.0,  1.0, /* 9 */
+    0.0,  1.0,  0.0, /* 10 */
+    0.0,  1.0, -1.0, /* 11 */
+    0.0,  0.0,  1.0  /* 12 */
+  };
+
+  if ( orientation >= 13 )
+    error( "Wrong orientation !" );
+
+  /* Select configuration */
+  float shift[3];
+  for ( k = 0 ; k < 3 ; ++k )
+    shift[k] = cell_shift[ 3 * orientation + k ];
+
   /* Init and fill the particle array. */
   if ((parts = (struct part *)malloc(sizeof(struct part) * N_parts * 2)) ==
       NULL)
@@ -2063,13 +2107,15 @@ void test_direct_neighbour(int N_parts, int orientation) {
     parts[k].id = k;
 
     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;
-    if (orientation > 1 && k >= N_parts) parts[k].x[2] += 1.;
+
+    /* Shift the second cell */
+    if (k >= N_parts) {
+      parts[k].x[0] += shift[0];
+      parts[k].x[1] += shift[1];
+      parts[k].x[2] += shift[2];
+    }
 
     parts[k].mass = ((double)rand()) / RAND_MAX;
     parts[k].a[0] = 0.0;
@@ -2083,9 +2129,9 @@ void test_direct_neighbour(int N_parts, int orientation) {
   left.loc[2] = 0.;
   left.h = 1.;
 
-  right.loc[0] = 1.;
-  right.loc[1] = (orientation > 0 ? 1 : 0);
-  right.loc[2] = (orientation > 1 ? 1 : 0);
+  right.loc[0] = shift[0];
+  right.loc[1] = shift[1];
+  right.loc[2] = shift[2];
   right.h = 1.;
 
   /* Put the particles in the cell */
@@ -2153,33 +2199,16 @@ void test_direct_neighbour(int N_parts, int orientation) {
   double *position;
   if ((position = (double *)malloc(sizeof(double) * N_parts * 2)) == NULL)
     error("Failed to allocate position buffer.");
+  
+  float midPoint = shift[0]*shift[0] + shift[1]*shift[1] + shift[2]*shift[2];
+  midPoint += ( shift[0] < 0 ? -1 : 0);
+  midPoint += ( shift[1] < 0 ? -1 : 0);
+  midPoint += ( shift[2] < 0 ? -1 : 0);
+  message( "MidPoint: %f", midPoint );
 
-  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;
-    }
-  }
+  for (k = 0; k < 2 * N_parts; ++k)
+    position[k] = parts[k].x[0] * shift[0] + parts[k].x[1] * shift[1] + parts[k].x[2] * shift[2] - midPoint;
+    
 
   /* Now, output everything */
   char fileName[100];
@@ -2218,7 +2247,7 @@ void test_direct_neighbour(int N_parts, int orientation) {
 
 int main(int argc, char *argv[]) {
 
-  int c, nr_threads;
+  int c, k, nr_threads;
   int N = 1000, runs = 1;
   char fileName[100] = {0};
   int N_parts = 0;
@@ -2283,9 +2312,8 @@ int main(int argc, char *argv[]) {
             N_parts);
 
     /* Run the test */
-    test_direct_neighbour(N_parts, 0);
-    test_direct_neighbour(N_parts, 1);
-    test_direct_neighbour(N_parts, 2);
+    for ( k = 0 ; k < 13 ; ++k )
+      test_direct_neighbour(N_parts, k);
 
   } else {