 	float rho = 0;
 	float scorr = 0;
 
-	FOR_EACH_NEIGHBOUR
+	vec3 grad_pi_Ci = vec3 (0, 0, 0);
+	FOR_EACH_NEIGHBOUR(j)
 	{			
 		// compute rho_i (equation 2)
-		float len = length (particle.position - particles[neighbourid].position);
+		float len = length (particle.position - particles[j].position);
 		float tmp = Wpoly6 (len);
 		rho += tmp;
 	

};

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

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

