 //--------------------------------------------------------------------------------------------------------------------------------
 // collide a particle against all other particles in a given cell
 __device__
-real_ collideCellDensityReInit_F1(
+void collideCellDensityReInit(
+		real_ & densityShare,
+		real_ & denominator,
 		int3 gridPos,
 		uint index,
 		real3 posRadA,

				real3 dist3 = Distance(posRadA, posRadB);
				real_ d = length(dist3);
				if (d > RESOLUTION_LENGTH_MULT * paramsD.HSML) continue;
				densityShare += velMasB.w * W3(d); //optimize it ?$
				//if (fabs(W3(d)) < .00000001) {printf("good evening, distance %f %f %f\n", dist3.x, dist3.y, dist3.z);
				//printf("posRadA %f %f %f, posRadB, %f %f %f\n", posRadA.x, posRadA.y, posRadA.z, posRadB.x, posRadB.y, posRadB.z);
				//}
			}
		}
	}
	return densityShare;
}
//--------------------------------------------------------------------------------------------------------------------------------
// collide a particle against all other particles in a given cell

	uint originalIndex = gridMarkerIndex[index];

	real_ newDensity = densityShare + velMasA.w * W3(0); //?$ include the particle in its summation as well
	real_ newDenominator = denominator + velMasA.w * W3(0) / rhoPreMuA.x;
	if (rhoPreMuA.w < 0) {
//		rhoPreMuA.x = newDensity; // old version
		rhoPreMuA.x = newDensity/newDenominator; // correct version
	}
	rhoPreMuA.y = Eos(rhoPreMuA.x, rhoPreMuA.w);
	oldRhoPreMu[originalIndex] = rhoPreMuA;
