 	return (-3 * 4.774648292756860 * tmp * tmp) * r;
 }
 
+// offsets to find neighbouring cells in one dimensional array
+const int gridoffsets[27] = {
+	-GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH - 1,		// -1, -1, -1
+	-GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH,			// -1, -1, 0
+	-GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH + 1,		// -1, -1, 1
+	-GRID_WIDTH * GRID_HEIGHT - 1,					// -1, 0, -1
+	-GRID_WIDTH * GRID_HEIGHT,						// -1, 0, 0
+	-GRID_WIDTH * GRID_HEIGHT + 1,					// -1, 0, 1
+	-GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH - 1,		// -1, 1, -1
+	-GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH,			// -1, 1, 0
+	-GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH + 1,		// -1, 1, 1
+	-GRID_WIDTH - 1,								// 0, -1, -1
+	-GRID_WIDTH,									// 0, -1, 0
+	-GRID_WIDTH + 1,								// 0, -1, 1
+	-1,												// 0, 0, -1
+	0,												// 0, 0, 0
+	1,												// 0, 0, 1
+	GRID_WIDTH - 1,									// 0, 1, -1
+	GRID_WIDTH,										// 0, 1, 0
+	GRID_WIDTH + 1,									// 0, 1, 1
+	GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH - 1,		// 1, -1, -1
+	GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH,			// 1, -1, 0
+	GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH + 1,		// 1, -1, 1
+	GRID_WIDTH * GRID_HEIGHT - 1,					// 1, 0, -1
+	GRID_WIDTH * GRID_HEIGHT,						// 1, 0, 0
+	GRID_WIDTH * GRID_HEIGHT + 1,					// 1, 0, 1
+	GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH - 1,		// 1, 1, -1
+	GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH,			// 1, 1, 0
+	GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH + 1		// 1, 1, 1
+};
+
+// convenience macros to iterate over all neighbouring particles
+// within the loop "neighbourid" contains the index of the current
+// neighbourparticle in the particles[] array
+#define FOR_EACH_NEIGHBOUR 	for (int o = 0; o < 27; o++) {\
+		int ngridid = gridid + gridoffsets[o];\
+		if (ngridid < 0 || ngridid > MAX_GRID_ID)\
+			continue;\
+		for (uint gridparticle = 0; gridparticle < gridcounters[ngridid]; gridparticle++) {\
+			uint neighbourid = gridcells[ngridid].particleids[gridparticle];
+#define END_FOR_EACH_NEIGHBOUR 	}}
+
 void main (void)
 {
 	uint particleid;

	barrier ();
	memoryBarrierBuffer ();
	
	uint neighbours[8*27];
	uint num_neighbours = 0;
	
	const int gridoffsets[27] = {
		-GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH - 1,		// -1, -1, -1
		-GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH,			// -1, -1, 0
		-GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH + 1,		// -1, -1, 1
		-GRID_WIDTH * GRID_HEIGHT - 1,					// -1, 0, -1
		-GRID_WIDTH * GRID_HEIGHT,						// -1, 0, 0
		-GRID_WIDTH * GRID_HEIGHT + 1,					// -1, 0, 1
		-GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH - 1,		// -1, 1, -1
		-GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH,			// -1, 1, 0
		-GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH + 1,		// -1, 1, 1
		-GRID_WIDTH - 1,								// 0, -1, -1
		-GRID_WIDTH,									// 0, -1, 0
		-GRID_WIDTH + 1,								// 0, -1, 1
		-1,												// 0, 0, -1
		0,												// 0, 0, 0
		1,												// 0, 0, 1
		GRID_WIDTH - 1,									// 0, 1, -1
		GRID_WIDTH,										// 0, 1, 0
		GRID_WIDTH + 1,									// 0, 1, 1
		GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH - 1,		// 1, -1, -1
		GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH,			// 1, -1, 0
		GRID_WIDTH * GRID_HEIGHT - GRID_WIDTH + 1,		// 1, -1, 1
		GRID_WIDTH * GRID_HEIGHT - 1,					// 1, 0, -1
		GRID_WIDTH * GRID_HEIGHT,						// 1, 0, 0
		GRID_WIDTH * GRID_HEIGHT + 1,					// 1, 0, 1
		GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH - 1,		// 1, 1, -1
		GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH,			// 1, 1, 0
		GRID_WIDTH * GRID_HEIGHT + GRID_WIDTH + 1		// 1, 1, 1
	};
	
	for (int o = 0; o < 27; o++)
	{
		int ngridid = gridid + gridoffsets[o];
		if (ngridid < 0 || ngridid > MAX_GRID_ID)
			continue;
		for (uint gridparticle = 0; gridparticle < gridcounters[ngridid]; gridparticle++)
		{
			uint neighbourid = gridcells[ngridid].particleids[gridparticle];
			neighbours[num_neighbours] = neighbourid;
			num_neighbours++;
		}
		
	}

	for (uint solver = 0; solver < 1; solver++) {

	barrier ();

	particles[particleid].position = particle.position;
	}

	particles[particleid].velocity = 0.5 * (particle.position - oldposition) / timestep;
}
