diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index 033e496864a328e26bf2289b819b79ae492f69be..4f7abfae1899a3fd10036f82efb849ce0c7fa128 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -33,6 +33,9 @@ rcParams.update(params)
rc('font', family='serif')
import sys
import os
+from scipy import stats
+
+dist_cutoff_ratio=1.2
print "Plotting..."
@@ -68,19 +71,28 @@ axis = [
#names = ["side", "edge", "corner"]
-for orientation in range( 26 ):
-# for jjj in range(2):
-# if jjj == 0:
-# orientation = 0
-# if jjj == 1:
-# orientation = 8
+#for orientation in range( 26 ):
+for jjj in range(3):
+ if jjj == 0:
+ orientation = 0
+ if jjj == 1:
+ orientation = 1
+ if jjj == 2:
+ orientation = 4
+
# Read Quickshed accelerations
data=loadtxt( "interaction_dump_%d.dat"%orientation )
id = data[:,0]
pos = data[:,2]
pos -= mean(pos)
-
+
+ x = data[:,3]
+ y = data[:,4]
+ z = data[:,5]
+
+ dist = data[:,12]
+
accx_u=data[:,6]
accy_u=data[:,7]
accz_u=data[:,8]
@@ -112,56 +124,95 @@ for orientation in range( 26 ):
stdz_s = std(errz_s[abs(errz_s) < 0.1])
+# sample_pos = pos[dist<0.2]
+# sample_x = e_errx_s[dist<0.2]
+# sample_y = e_erry_s[dist<0.2]
+# sample_z = e_errz_s[dist<0.2]
+
+
+# numBins = 100
+# binEdges = linspace(-1.5-1.5/numBins, 1.5+1.5/numBins, numBins+2)
+# bins = linspace(-1.5, 1.5, numBins+1)
+# corr_x, a, b = stats.binned_statistic(sample_pos, sample_x, statistic='mean', bins=binEdges)
+# corr_y, a, b = stats.binned_statistic(sample_pos, sample_y, statistic='mean', bins=binEdges)
+# corr_z, a, b = stats.binned_statistic(sample_pos, sample_z, statistic='mean', bins=binEdges)
+
+# a,b, sample_bin = stats.binned_statistic(pos, pos, statistic='mean', bins=binEdges)
+# sample_bin -= 2
+
+
+# # for j in range(size(pos)):
+# # e_errx_s /= corr_x[sample_bin[j]]
+# # e_erry_s /= corr_y[sample_bin[j]]
+# # e_errz_s /= corr_z[sample_bin[j]]
+
# Plot -------------------------------------------------------
figure(frameon=True)
+
subplot(311, title="Acceleration along X")
- #plot(id[abs(errx_s) > 0.001], e_errx_s , 'ro')
- plot(pos, e_errx_s , 'ro')
- # for j in range(size(pos)):
- # if ( pos[j-1] != pos[j] ):
- # text(pos[j], e_errx_s[j]-0.005, "%d"%id[j], fontsize=10)
- # else:
- # text(pos[j], e_errx_s[j]-0.015, "%d"%id[j], fontsize=10)
- 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)
+ scatter(pos[dist<1.], e_errx_s[dist<1.] ,c=dist[dist<1.], marker='o', s=1, linewidth=0, cmap='jet')
+ plot([-dist_cutoff_ratio/2., -dist_cutoff_ratio/2.], [-0.01, 0.01], 'k--')
+ plot([dist_cutoff_ratio/2., dist_cutoff_ratio/2.], [-0.01, 0.01], 'k--')
+ text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
- ylim(-0.05, 0.05)
+ ylim(-0.03, 0.03)
grid()
subplot(312, title="Acceleration along Y")
- #plot(id[abs(erry_s) > 0.001], e_erry_s , 'ro')
- plot(pos, e_erry_s , 'ro')
- # for j in range(size(pos)):
- # if ( pos[j-1] != pos[j] ):
- # text(pos[j], e_errx_s[j]-0.005, "%d"%id[j], fontsize=10)
- # else:
- # text(pos[j], e_errx_s[j]-0.015, "%d"%id[j], fontsize=10)
- 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)
+ scatter(pos[dist<1.], e_erry_s[dist<1.] , c=dist[dist<1.], marker='o', s=1, linewidth=0, cmap='jet')
+ plot([-dist_cutoff_ratio/2., -dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
+ plot([dist_cutoff_ratio/2., dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
+ text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
- ylim(-0.05, 0.05)
+ ylim(-0.03, 0.03)
grid()
subplot(313, title="Acceleration along Z")
- #plot(id[abs(errz_s) > 0.001], e_errz_s , 'ro', label="Sorted")
- plot(pos, e_errz_s , 'ro', label="Sorted")
- # for j in range(size(pos)):
- # if ( pos[j-1] != pos[j] ):
- # text(pos[j], e_errx_s[j]-0.005, "%d"%id[j], fontsize=10)
- # else:
- # text(pos[j], e_errx_s[j]-0.015, "%d"%id[j], fontsize=10)
- 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")
+ scatter(pos[dist<1.], e_errz_s[dist<1.] , c=dist[dist<1.], label="Sorted", marker='o', s=1, linewidth=0, cmap='jet')
+ plot([-dist_cutoff_ratio/2., -dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
+ plot([dist_cutoff_ratio/2., dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
+ text( 0., 0.005, "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")
xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
- ylim(-0.05, 0.05)
- #xlim(0, size(id)-1)
+ ylim(-0.03, 0.03)
grid()
- savefig("1quadrupole_accelerations_relative_%d.png"%orientation)
+ savefig("quadrupole_accelerations_relative_%d.png"%orientation)
close()
+
+ figure(frameon=True)
+
+ subplot(311, title="Acceleration along X")
+ scatter(dist, e_errx_s )
+ text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+ #xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
+ #ylim(-0.01, 0.01)
+ grid()
+
+ subplot(312, title="Acceleration along Y")
+ scatter(dist, e_erry_s )
+ text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+ #xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
+ #ylim(-0.01, 0.01)
+ grid()
+
+ subplot(313, title="Acceleration along Z")
+ scatter(dist, e_errz_s , label="Sorted")
+ text( 0., 0.005, "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")
+ #xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
+ #ylim(-0.01, 0.01)
+ grid()
+
+ savefig("error_distance_%d.png"%orientation)
+ close()
+
+
+
# # Plot -------------------------------------------------------
# figure(frameon=True)
@@ -207,22 +258,22 @@ for orientation in range( 26 ):
- # figure(frameon=True)
- # subplot(311, title="Acceleration along X")
- # 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)
+ figure(frameon=True)
+ subplot(311, title="Acceleration along X")
+ 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(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)
+ 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_%d.png"%orientation)
- # close()
+ savefig("histogram_%d.png"%orientation)
+ close()
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 4b8469bf83cc8c5b378ce4a30822b3f60d7f5d50..185377b3b663afb413f1b6592cf361cfd2815e7e 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -190,6 +190,7 @@ static inline void multipole_iact(struct multipole *d, const float *x, float* a)
w = const_G * ir * ir * ir * d->mass;
for (k = 0; k < 3; k++) a[k] -= w * dx[k];
+
/* Compute some helpful terms */
float w1 = const_G * ir * ir * ir * ir * ir;
float w2 = -2.5f * const_G * ir * ir * ir * ir * ir * ir * ir;
@@ -2207,18 +2208,30 @@ void test_direct_neighbour(int N_parts, int orientation) {
}
/* Position along the axis */
- double *position;
+ double *position, *ortho_dist;
if ((position = (double *)malloc(sizeof(double) * N_parts * 2)) == NULL)
error("Failed to allocate position buffer.");
+ if ((ortho_dist = (double *)malloc(sizeof(double) * N_parts * 2)) == NULL)
+ error("Failed to allocate orthogonal distance buffer.");
float w = 1.0f / sqrtf(shift[0]*shift[0] + shift[1]*shift[1] + shift[2]*shift[2]);
shift[0] *= w;
shift[1] *= w;
shift[2] *= w;
- 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];
+ for (k = 0; k < 2 * N_parts; ++k) {
+ float dx = parts[k].x[0] * shift[0];
+ float dy = parts[k].x[1] * shift[1];
+ float dz = parts[k].x[2] * shift[2];
+
+ position[k] = dx + dy + dz;
+
+ dx = (parts[k].x[1] - 0.5f)*shift[2] - (parts[k].x[2] - 0.5f)*shift[1];
+ dy = (parts[k].x[2] - 0.5f)*shift[0] - (parts[k].x[0] - 0.5f)*shift[2];
+ dz = (parts[k].x[0] - 0.5f)*shift[1] - (parts[k].x[1] - 0.5f)*shift[0];
+
+ ortho_dist[k] = sqrtf(dx*dx + dy*dy + dz*dz);
+ }
/* Now, output everything */
char fileName[100];
@@ -2226,12 +2239,13 @@ void test_direct_neighbour(int N_parts, int 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");
+ "# ID m r x y z a_u.x a_u.y a_u.z a_s.x a_s.y a_s.z ortho\n");
for (k = 0; k < 2 * N_parts; k++) {
- fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id,
+ fprintf(file, "%d %e %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]);
+ parts[k].a_exact[2], parts[k].a[0], parts[k].a[1], parts[k].a[2],
+ ortho_dist[k]);
}
fclose(file);
@@ -2490,11 +2504,11 @@ int main(int argc, char *argv[]) {
N_parts);
/* Run the test */
- for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k);
- /* k++; */
- /* test_direct_neighbour(N_parts, 0); */
- /* test_direct_neighbour(N_parts, 8); */
-
+ //for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k);
+ test_direct_neighbour(N_parts, 0);
+ test_direct_neighbour(N_parts, 1);
+ test_direct_neighbour(N_parts, 4);
+ k++;
/* Dump arguments */
message("Interacting one cell with %d particles with its 27 neighbours"