 	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;
-}
+	particles[particleid].velocity = 0.5 * (particle.position - oldposition) / timestep;
+}

	return (-3 * 4.774648292756860 * tmp * tmp) * r;
}

void main (void)
{
	uint particleid;

	
	float sum_k_grad_Ci = 0;
	float lambda = 0;
	float rho = 0;
	float scorr = 0;

	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;
			
		// 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[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);
	
