diff --git a/examples/EAGLE_12/check_ngbs.py b/examples/EAGLE_12/check_ngbs.py
new file mode 100644
index 0000000000000000000000000000000000000000..b40dbd7e923b78677321fdeaba65d328e7229144
--- /dev/null
+++ b/examples/EAGLE_12/check_ngbs.py
@@ -0,0 +1,248 @@
+import h5py as h
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import os.path
+
+kernel_gamma = 1.825742
+kernel_gamma2 = kernel_gamma * kernel_gamma
+kernel_gamma_dim = np.power(kernel_gamma,3)
+hydro_dimension_unit_sphere = 4. * np.pi / 3.
+kernel_norm = hydro_dimension_unit_sphere * kernel_gamma_dim
+error = False
+
+inputFile1 = ""
+inputFile2 = ""
+
+# Check list of density neighbours and check that they are correct.
+def check_density_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
+
+ error_val = False
+
+ for k in range(0,num_invalid):
+
+ # Filter neighbour lists for valid particle ids
+ filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
+ filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
+
+ # Check neighbour lists for differences
+ id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
+
+ pid = pids[mask][k]
+
+ # Loop over discrepancies and check if they are actually neighbours
+ for pjd in id_list:
+ pi_pos = pos[np.where(pids == pid)]
+ pj_pos = pos[np.where(pids == pjd)]
+
+ hi = h[np.where(pids == pid)]
+
+ dx = pi_pos[0][0] - pj_pos[0][0]
+ dy = pi_pos[0][1] - pj_pos[0][1]
+ dz = pi_pos[0][2] - pj_pos[0][2]
+
+ r2 = dx*dx + dy*dy + dz*dz
+
+ hig2 = hi*hi*kernel_gamma2
+
+ diff = abs(r2 - hig2)
+
+ print "Particle {} is missing {}, hig2: {}, r2: {}, |r2 - hig2|: {}".format(pid,pjd,hig2, r2, diff)
+
+ if diff < acc * hig2:
+ print "Missing interaction due to precision issue will be ignored."
+ else:
+ error_val = True
+
+ return error_val
+
+# Check list of force neighbours and check that they are correct.
+def check_force_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
+
+ error_val = False
+
+ for k in range(0,num_invalid):
+
+ # Filter neighbour lists for valid particle ids
+ filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
+ filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
+
+ # Check neighbour lists for differences
+ id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
+
+ pid = pids[mask][k]
+
+ # Loop over discrepancies and check if they are actually neighbours
+ for pjd in id_list:
+ pi_pos = pos[np.where(pids == pid)]
+ pj_pos = pos[np.where(pids == pjd)]
+
+ hi = h[np.where(pids == pid)]
+ hj = h[np.where(pids == pjd)]
+
+ dx = pi_pos[0][0] - pj_pos[0][0]
+ dy = pi_pos[0][1] - pj_pos[0][1]
+ dz = pi_pos[0][2] - pj_pos[0][2]
+
+ r2 = dx*dx + dy*dy + dz*dz
+
+ hig2 = hi*hi*kernel_gamma2
+ hjg2 = hj*hj*kernel_gamma2
+
+ diff = abs(r2 - max(hig2, hjg2))
+
+ print "Particle {} is missing {}, hig2: {}, hjg2: {}, r2: {}, |r2 - max(hig2,hjg2)|: {}".format(pid,pjd,hig2, hjg2, r2, diff)
+
+
+ if diff < acc * max(hig2,hjg2):
+ print "Missing interaction due to precision issue will be ignored."
+ else:
+ error_val = True
+
+ return error_val
+
+# Parse command line arguments
+if len(sys.argv) < 3:
+ print "Error: pass input files as arguments"
+ sys.exit()
+else:
+ inputFile1 = sys.argv[1]
+ inputFile2 = sys.argv[2]
+ if os.path.exists(inputFile1) != 1:
+ print "\n{} does not exist!\n".format(inputFile1)
+ sys.exit()
+ if os.path.exists(inputFile2) != 1:
+ print "\n{} does not exist!\n".format(inputFile2)
+ sys.exit()
+
+# Open input files
+file_naive = h.File(inputFile1, "r")
+file_sort = h.File(inputFile2, "r")
+
+# Read input file fields
+ids_naive = file_naive["/PartType0/ParticleIDs"][:]
+ids_sort = file_sort["/PartType0/ParticleIDs"][:]
+
+h_naive = file_naive["/PartType0/SmoothingLength"][:]
+h_sort = file_sort["/PartType0/SmoothingLength"][:]
+
+pos_naive = file_naive["/PartType0/Coordinates"][:,:]
+pos_sort = file_sort["/PartType0/Coordinates"][:,:]
+
+num_density_naive = file_naive["/PartType0/Num_ngb_density"][:]
+num_density_sort = file_sort["/PartType0/Num_ngb_density"][:]
+
+num_force_naive = file_naive["/PartType0/Num_ngb_force"][:]
+num_force_sort = file_sort["/PartType0/Num_ngb_force"][:]
+
+neighbour_ids_density_naive = file_naive["/PartType0/Ids_ngb_density"][:]
+neighbour_ids_density_sort = file_sort["/PartType0/Ids_ngb_density"][:]
+
+neighbour_ids_force_naive = file_naive["/PartType0/Ids_ngb_force"][:]
+neighbour_ids_force_sort = file_sort["/PartType0/Ids_ngb_force"][:]
+
+#wcount_naive = file_naive["/PartType0/Wcount"][:]
+#wcount_sort = file_sort["/PartType0/Wcount"][:]
+#
+#wcount_naive = wcount_naive * np.power(h_naive,3) * kernel_norm
+#wcount_sort = wcount_sort * np.power(h_sort,3) * kernel_norm
+
+# Cross check
+print " Min Mean Max "
+print " ---------------------"
+print "Ngbs density naiv: ", np.min(num_density_naive), np.mean(num_density_naive), np.max(num_density_naive)
+print "Ngbs density sort: ", np.min(num_density_sort), np.mean(num_density_sort), np.max(num_density_sort)
+print "Ngbs force naiv: ", np.min(num_force_naive), np.mean(num_force_naive), np.max(num_force_naive)
+print "Ngbs force sort: ", np.min(num_force_sort), np.mean(num_force_sort), np.max(num_force_sort)
+#print "Wcount naiv: ", np.min(wcount_naive), np.mean(wcount_naive), np.max(wcount_naive)
+#print "Wcount sort: ", np.min(wcount_sort), np.mean(wcount_sort), np.max(wcount_sort)
+
+# Sort
+index_naive = np.argsort(ids_naive)
+index_sort = np.argsort(ids_sort)
+
+num_density_naive = num_density_naive[index_naive]
+num_density_sort = num_density_sort[index_sort]
+num_force_naive = num_force_naive[index_naive]
+num_force_sort = num_force_sort[index_sort]
+ids_naive = ids_naive[index_naive]
+ids_sort = ids_sort[index_sort]
+neighbour_ids_density_naive = neighbour_ids_density_naive[index_naive]
+neighbour_ids_density_sort = neighbour_ids_density_sort[index_sort]
+neighbour_ids_force_naive = neighbour_ids_force_naive[index_naive]
+neighbour_ids_force_sort = neighbour_ids_force_sort[index_sort]
+#wcount_naive = wcount_naive[index_naive]
+#wcount_sort = wcount_sort[index_sort]
+h_naive = h_naive[index_naive]
+h_sort = h_sort[index_sort]
+pos_naive = pos_naive[index_naive]
+pos_sort = pos_sort[index_sort]
+
+# First check
+print "\n Min Max"
+print " ----------"
+print "Differences for density: ", min(num_density_naive - num_density_sort), max(num_density_naive - num_density_sort)
+print "Differences for force: ", min(num_force_naive - num_force_sort), max(num_force_naive - num_force_sort)
+
+# Get the IDs that are different
+mask_density = num_density_naive != num_density_sort
+mask_force = num_force_naive != num_force_sort
+num_invalid_density = np.sum(mask_density)
+num_invalid_force = np.sum(mask_force)
+
+print "\nNum non-zero density: ", num_invalid_density
+print "Num non-zero force: ", num_invalid_force
+
+print "\nParticle IDs with incorrect densities"
+print "----------------------------------------"
+print ids_naive[mask_density]
+
+# Check density neighbour lists
+error += check_density_neighbours(ids_naive, neighbour_ids_density_naive,
+ neighbour_ids_density_sort, mask_density, pos_naive, h_naive,
+ num_invalid_density, 2e-6)
+
+print "Num of density interactions", inputFile1
+print num_density_naive[mask_density]
+
+print "Num of density interactions", inputFile2
+print num_density_sort[mask_density]
+
+print "\nParticle IDs with incorrect forces"
+print "------------------------------------"
+print ids_naive[mask_force]
+
+# Check force neighbour lists
+error += check_force_neighbours(ids_naive, neighbour_ids_force_naive,
+ neighbour_ids_force_sort, mask_force, pos_naive, h_naive,
+ num_invalid_force, 2e-6)
+
+print "Num of force interactions", inputFile1
+print num_force_naive[mask_force]
+
+#print "Smoothing lengths", inputFile1
+#print h_naive[mask_force]
+
+print "Num of force interactions", inputFile2
+print num_force_sort[mask_force]
+
+#print "Smoothing lengths", inputFile2
+#print h_sort[mask_force]
+
+# Statistics of h difference
+h_relative = (h_naive - h_sort) / h_naive
+print "h statistics: {} {} (Min, 1st Percentile)".format(np.min(h_relative), np.percentile(h_relative,1))
+print "h statistics: {} {} (Mean, Median)".format(np.mean(h_relative), np.median(h_relative))
+print "h statistics: {} {} (Max, 99th Percentile)".format(np.max(h_relative), np.percentile(h_relative, 99))
+
+if error:
+ print "\n------------------"
+ print "Differences found."
+ print "------------------"
+ exit(1)
+else:
+ print "\n---------------------"
+ print "No differences found."
+ print "---------------------"
+ exit(0)