 	int3 gridPos = calcGridPos(posRadA);
 
 	real_ densityShare = 0.0f;
+	real_ denominator = 0.0f;
 	// examine neighbouring cells
 	for (int z = -1; z <= 1; z++) {
 		for (int y = -1; y <= 1; y++) {
 			for (int x = -1; x <= 1; x++) {
 				int3 neighbourPos = gridPos + I3(x, y, z);
-				densityShare += collideCellDensityReInit_F1(neighbourPos, index, posRadA, sortedPosRad, sortedVelMas, sortedRhoPreMu, cellStart,
+				collideCellDensityReInit(densityShare, denominator, neighbourPos, index, posRadA, sortedPosRad, sortedVelMas, sortedRhoPreMu, cellStart,
 						cellEnd);
 			}
 		}

	//?c2 printf("grid pos %d %d %d \n", gridPos.x, gridPos.y, gridPos.z);
	uint gridHash = calcGridHash(gridPos);
	// get start of bucket for this cell
	real_ densityShare = 0.0f;

	uint startIndex = FETCH(cellStart, gridHash);
	if (startIndex != 0xffffffff) { // cell is not empty

	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;
