 	GridCell gridcells[];
 };
 
+float Wpoly6 (float r)
+{
+	if (r > 1)
+		return 0;
+	float tmp = 1 - r * r;
+	return 1.56668147106 * tmp * tmp * tmp;
+}
+
+float Wspiky (float r)
+{
+	if (r > 1)
+		return 0;
+	float tmp = 1 - r;
+	return 4.774648292756860 * tmp * tmp * tmp;
+}
+
+float gradWspiky (float r)
+{
+	if (r > 1)
+		return 0;
+	float tmp = 1 - r;
+	return -3 * 4.774648292756860 * tmp * tmp;
+}
+
 void main (void)
 {
 	uint particleid;

	barrier ();
	memoryBarrierBuffer ();
	
	particle.position = clamp (particle.position, vec3 (0, 0, 0), vec3 (GRID_WIDTH, GRID_HEIGHT, GRID_DEPTH));

	ivec3 gridoffset;
	for (gridoffset.x = -1; gridoffset.x <= 1; gridoffset.x++)
	for (gridoffset.y = -1; gridoffset.y <= 1; gridoffset.y++)
	for (gridoffset.z = -1; gridoffset.z <= 1; gridoffset.z++)
	{
		ivec3 ngrid = grid + gridoffset;
		
		int ngridid = ngrid.z * GRID_WIDTH * GRID_HEIGHT + ngrid.y * GRID_WIDTH + ngrid.x;
		
		for (uint i = 0; i < gridcounters[ngridid]; i++)
		{
			if (gridoffset == ivec3 (0, 0, 0) && i == pos)
				continue;
			
			ParticleInfo neighbour = gridcells[ngridid].particleids[i];
			
			vec3 impactdir = neighbour.position - particle.position;
			float impactlen = length (impactdir);
			if (impactlen < 1)
			{
				// collisiion
				impactdir /= impactlen;
				particle.position = 0.5 * (particle.position + neighbour.position - impactdir);
			}
		}
	}
		
	barrier ();
	
	particles[particleid] = particle;
}

	ParticleInfo particle = particles[particleid];
	
	// gravity
	vec3 velocity = (particle.position - particle.oldposition) / lasttimestep;
	velocity += 1000 * vec3 (0, -1, 0) * timestep;
	
	// save old position and predict new position
