 	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;
+//		rhoPreMuA.x = newDensity; // old version
+		rhoPreMuA.x = newDensity/newDenominator; // correct version
 	}
 	rhoPreMuA.y = Eos(rhoPreMuA.x, rhoPreMuA.w);
 	oldRhoPreMu[originalIndex] = rhoPreMuA;

	int3 gridPos = calcGridPos(posRadA);

	real_ densityShare = 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,
						cellEnd);
			}
		}

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