 	
 	float sum_k_grad_Ci = 0;
 	float lambda = 0;
-	
 	float rho = 0;
-	
 	float scorr = 0;
-	for (uint j = 0; j < num_neighbours; j++)
-	{
-		float len = length (particle.position - particles[neighbours[j]].position);
+
+	FOR_EACH_NEIGHBOUR
+	{			
+		// compute rho_i (equation 2)
+		float len = length (particle.position - particles[neighbourid].position);
 		float tmp = Wpoly6 (len);
 		rho += tmp;
-		
+	
+		// compute scorr (equation 13) - currently unused
 		tmp = -tmp / Wpoly6 (0.1);
 		tmp *= tmp;
 		scorr += -0.1 * tmp * tmp;
-	}
-	
-	for (uint k = 0; k < num_neighbours; k++)
-	{
+			
+		// sum gradients of Ci (equation 8 and parts of equation 9)
 		vec3 grad_pk_Ci = vec3 (0, 0, 0);
 		// plugging equation 8 into the sum in equation 9 is actually equivalent
 		// to omitting the values for k=i and instead summing over the values for k=j
 		// with a factor of 2, as direction and sign does not matter for the squared length,
 		// that is summed up in equation 9.
-		grad_pk_Ci = gradWspiky (particle.position - particles[neighbours[k]].position);
+		grad_pk_Ci = gradWspiky (particle.position - particles[neighbourid].position);
 		grad_pk_Ci /= rho_0;
 		sum_k_grad_Ci += 2 * dot (grad_pk_Ci, grad_pk_Ci);
 	}
-
+	END_FOR_EACH_NEIGHBOUR
+	
+	// compute lambda_i (equations 1 and 9)
 	float C_i = rho / rho_0 - 1;
 	lambda = -C_i / (sum_k_grad_Ci + 0.01);
 	

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

//	vec3 flipvelocity = vec3 (greaterThanEqual (particle.position, vec3 (0, 0, 0))) * vec3(lessThanEqual (particle.position, GRID_SIZE));
	particle.velocity = 0.5 * (particle.position - oldposition) / timestep;
//	particle.velocity *= -2 * (flipvelocity - 0.5);
/*
	vec3 v_i = particle.velocity;
	for (uint j = 0; j < num_neighbours; j++)
	{
		vec3 v_j = particles[neighbours[j]].velocity;
		particle.velocity += 0.01 * (v_j - v_i) * Wpoly6 (particle.position - particles[neighbours[j]].position);
	}*/

//	barrier ();
	
	particles[particleid].velocity = particle.velocity;
}

	
	vec3 deltap = vec3 (0, 0, 0);
	
	FOR_EACH_NEIGHBOUR
	{
		// accumulate position corrections (part of equation 12)
		deltap += (lambda + lambdas[neighbourid]) * gradWspiky (particle.position - particles[neighbourid].position);
	}
	END_FOR_EACH_NEIGHBOUR

	particle.position += deltap / rho_0;
	
	// collision detection begin
