diff --git a/examples/plot_sorted_all.py b/examples/plot_sorted_all.py
index 7defd5aa911e65500c3a132acdfd5d81ac06b4fc..ab9cefc5621ccc674d2037890a1d495e4db1c34b 100644
--- a/examples/plot_sorted_all.py
+++ b/examples/plot_sorted_all.py
@@ -61,14 +61,14 @@ accz_bh=data[:,13]
# Build error ------------------------------------------------
-
-errx_s = (accx_s - accx_u ) / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
-erry_s = (accy_s - accy_u ) / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
-errz_s = (accz_s - accz_u ) / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
+inv_acc_tot = 1.0 / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
+errx_s = (accx_s - accx_u ) * inv_acc_tot
+erry_s = (accy_s - accy_u ) * inv_acc_tot
+errz_s = (accz_s - accz_u ) * inv_acc_tot
-errx_bh = (accx_bh - accx_u ) / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
-erry_bh = (accy_bh - accy_u ) / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
-errz_bh = (accz_bh - accz_u ) / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
+errx_bh = (accx_bh - accx_u ) * inv_acc_tot
+erry_bh = (accy_bh - accy_u ) * inv_acc_tot
+errz_bh = (accz_bh - accz_u ) * inv_acc_tot
e_errx_s = errx_s#[abs(errx_s) > 0.001]
e_erry_s = erry_s#[abs(erry_s) > 0.001]
@@ -94,10 +94,10 @@ figure(frameon=True)
subplot(311, title="Acceleration along X")
#plot(id[abs(errx_s) > 0.001], e_errx_s , 'ro')
-plot(pos_x, e_errx_s , 'ro')
plot(pos_x, e_errx_bh , 'bx', label='B-H')
+plot(pos_x, 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)
+ylim(-0.1, 0.1)
grid()
subplot(312, title="Acceleration along Y")
@@ -105,7 +105,7 @@ subplot(312, title="Acceleration along Y")
plot(pos_x, e_erry_bh , 'bx', label='B-H')
plot(pos_y, 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)
+ylim(-0.1, 0.1)
grid()
subplot(313, title="Acceleration along Z")
@@ -116,7 +116,7 @@ plot(pos_z, e_errz_s , 'ro', label="Sorted")
legend(loc="upper right")
-ylim(-0.2, 0.2)
+ylim(-0.1, 0.1)
#xlim(0, size(id)-1)
grid()
@@ -129,7 +129,7 @@ figure(frameon=True)
subplot(311, title="Acceleration along X")
#plot(id[abs(errx_s) > 0.001], e_errx_s , 'ro')
-plot(pos_x, accx_u, 'bx')
+plot(pos_x, accx_bh, 'bx')
plot(pos_x, accx_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(-250, 250)
@@ -137,9 +137,9 @@ grid()
subplot(312, title="Acceleration along Y")
#plot(id[abs(erry_s) > 0.001], e_erry_s , 'ro')
-plot(pos_y, accy_u , 'gs')
+#plot(pos_y, accy_u , 'gs')
plot(pos_y, accy_s , 'ro')
-plot(pos_z, accz_bh , 'bx', label="B-H")
+plot(pos_y, accz_bh , 'bx', label="B-H")
#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(-70, 70)
ylim(-250, 250)
@@ -149,7 +149,7 @@ subplot(313, title="Acceleration along Z")
#plot(id[abs(errz_s) > 0.001], e_errz_s , 'ro', label="Sorted")
plot(pos_z, accz_bh , 'bx', label="B-H")
plot(pos_z, accz_s , 'ro', label="Sorted")
-plot(pos_z, accz_u , 'gs', label="Unsorted")
+#plot(pos_z, accz_u , 'gs', label="Unsorted")
#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)
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 4dbbe35513cfe8f7d9eb3874972a66186cbecc88..6210df0d0fda1885599e9b7cd6b996314462ad76 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -231,8 +231,6 @@ static inline void multipole_add_matrix(struct multipole *d, const float *x, flo
for ( k = 0; k < 3; ++k)
dx[k] = x[k] - d->x[k] / d->mass;
- // message("x_i: %f %f %f dx: %f %f %f", x[0], x[1], x[2], dx[0], dx[1], dx[2]);
-
r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Now fill in the matrix */
@@ -273,7 +271,7 @@ static inline void multipole_iact(struct multipole *d, const float *x, float* a)
float inv_mass = 1. / d->mass;
- /* Compute the distance. */
+ /* Compute the distance to the CoM. */
for (r2 = 0.0f, k = 0; k < 3; k++) {
dx[k] = x[k] - d->x[k] * inv_mass;
r2 += dx[k] * dx[k];
@@ -286,8 +284,6 @@ static inline void multipole_iact(struct multipole *d, const float *x, float* a)
float a_m[3] = {0., 0., 0.};
for (k = 0; k < 3; k++) a_m[k] -= w * dx[k];
-
- //return;
/* Compute the acceleration from the quadrupole */
float a_q[3] = {0., 0., 0.};
@@ -303,7 +299,8 @@ static inline void multipole_iact(struct multipole *d, const float *x, float* a)
a_q[1] = w2 * sum * dx[1] + w1 * (2.*dx[0]*d->Q[1][0] + 2.*dx[1]*d->Q[1][1] + 2.*dx[2]*d->Q[1][2]);
a_q[2] = w2 * sum * dx[2] + w1 * (2.*dx[0]*d->Q[2][0] + 2.*dx[1]*d->Q[2][1] + 2.*dx[2]*d->Q[2][2]);
-
+
+ /* Get the total acceleration */
float a_t[3] = {0., 0., 0.};
for (k = 0; k < 3; k++) a_t[k] = a_m[k] + a_q[k];
@@ -1861,8 +1858,6 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
#endif
- //message("start");
-
/* Get the sorted indices and stuff. */
struct index *ind_i, *ind_j;
struct multipole dip[4];
@@ -1873,9 +1868,6 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
for (k = 0; k < (1 << num_orth_planes); k++) multipole_init(&dip[k]);
cjh = cj->h;
-
- //message("init OK");
-
/* Allocate and fill-in the local parts. */
count_i = ci->count;
count_j = cj->count;
@@ -1907,8 +1899,6 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
parts_j[j].mass = cj->parts[pjd].mass;
}
- //message("fill in ok");
-
#if ICHECK >= 0
for (k = 0; k < count_i; k++)
if (parts_i[k].id == ICHECK)
@@ -1987,8 +1977,6 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
multipole_reset(&dip[l]);
}
- //message(" i=%d m-pole reset", i)
-
/* Add any remaining particles to the COM. */
for (int jj = j; jj < count_j; jj++) {