diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index 1787c78a5cf7f1d2e6f88d09624080e16edead95..d66fec58514b706357492bb10d21d3222ad7a254 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -135,7 +135,6 @@ for orientation in range( 26 ):
     #    continue
     cell_edge = sqrt(sum(abs(array(this_axis))))# - sum(abs(this_axis[this_axis<0]))
     print this_axis
-    # print 2*cell_edge
     print "actual:", orientation, min(dist), max(dist), min_dist[0], max_dist[0]
     print "---"
     continue
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 2e244dacaf93493cc523218a4d72c8db0194f6a1..895d079680bef6b752791f10bc74cc41d5667201 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -1184,7 +1184,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
   struct part_local *parts_i, *parts_j;
   float xi[3];
   float dx[3], ai[3], mi, mj, r2, w, ir;
-  float min_dist, max_dist;
+  float min_dist, max_dist, length_dist, mid_dist;
 
 #ifdef SANITY_CHECKS
 
@@ -1194,10 +1194,14 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
 
 #endif
 
-  /* Get the sorted indices and stuff. */
+  /* Get the sorted indices and geometry. */
   struct index *ind_i, *ind_j;
-  struct multipole multi;
   get_axis(&ci, &cj, &ind_i, &ind_j, &min_dist, &max_dist);
+  length_dist = max_dist - min_dist;
+  mid_dist = min_dist + 0.5f * length_dist;
+
+  /* Initialise the multipole */
+  struct multipole multi;
   multipole_reset(&multi);
 
   /* Allocate and fill-in the local parts. */
@@ -1210,6 +1214,9 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
            (struct part_local *)alloca(sizeof(struct part_local) * count_j)) ==
           NULL)
     error("Failed to allocate local part arrays.");
+
+  /* Store only the particle position relative to a point close to the cells */
+  /* This point does not have to be in the middle cj->loc is good enough */
   float midpoint[3] = { cj->loc[0] , cj->loc[1], cj->loc[2] };
   for (i = 0; i < count_i; i++) {
     int pid = ind_i[i].ind;
@@ -1251,7 +1258,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
   /* Loop over all particles in ci... */
   for (i = count_i - 1; i >= 0; i--) {
 
-    /* Get a local copy of the distance along the axis. */
+    /* Get a local copy of the maximal distance along the axis. */
     float di = parts_i[i].d + d_max;
 
     /* Init the ith particle's data. */
@@ -1303,7 +1310,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
 
     } /* loop over every other particle. */
 
-    /* Add any remaining particles to the COM. */
+    /* Add any remaining particles to the multipole. */
     for (int jj = j; jj < count_j; jj++) {
       multipole_add(&multi, parts_j[jj].x, parts_j[jj].mass);
     }
@@ -1311,12 +1318,12 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
     /* Shrink count_j to the latest valid particle. */
     count_j = j;
 
-    /* Interact part_i with the center of mass. */
+    /* Interact part_i with the multipole. */
     multipole_iact(&multi, xi, ai);
 
 #if ICHECK >= 0
       if (parts_i[i].id == ICHECK)
-        message("Interaction with multipole");  //, parts_j[j].id );
+        message("Interaction with multipole"); 
 #endif
 
 #ifdef COUNTERS
@@ -1353,12 +1360,24 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
     /* Interact part_j with the COM. */
     multipole_iact(&multi, parts_j[j].x, parts_j[j].a);
 
+#if ICHECK >= 0
+      if (parts_j[j].id == ICHECK)
+        message("Interaction with multipole"); 
+#endif
+
+
 #ifdef COUNTERS
       ++count_direct_sorted_pm_j;
 #endif
     
   }
 
+  /* Apply the correction to cancel the bias */
+
+  for (i = 0; i < count_i; i++) {
+    message("Length = %f Mid = %f dist = %f", length_dist, mid_dist, parts_i[i].d);
+  }
+
   /* Copy the accelerations back to the original particles. */
   for (i = 0; i < count_i; i++) {
     int pid = ind_i[i].ind;