 			Real Wd = W3(d);
 			Real3 velMasB = FETCH(sortedVelMas, j);
 			//			deltaVDenom += mR4(
-			//					velMasB.w / rhoPresMuB.x * mR3(velMasB) * Wd,
-			//					velMasB.w / rhoPresMuB.x * Wd);
+			//					paramsD.markerMass / rhoPresMuB.x * velMasB * Wd,
+			//					paramsD.markerMass / rhoPresMuB.x * Wd);
 			//			deltaVDenom += mR4(
-			//					mR3(velMasB) * Wd,
+			//					velMasB * Wd,
 			//					Wd);
 
 			if (rhoPresMuB.w < -.1) { // only fluid pressure is used to update BCE pressure see Eq 27 of Adami, 2012 paper

					+ rhoPresMuB.y / (rhoPresMuB.x * rhoPresMuB.x)) * gradW
			+ paramsD.markerMass * (8.0f * multViscosity) * paramsD.mu0
					* pow(rhoPresMuA.x + rhoPresMuB.x, Real(-2))
					* rAB_Dot_GradW_OverDist * velMasA - velMasB;
	Real derivRho = rhoPresMuA.x * paramsD.markerMass / rhoPresMuB.x
			* dot(vel_XSPH_A - vel_XSPH_B, gradW);
	//	Real zeta = 0;//.05;//.1;
	//	Real derivRho = rhoPresMuA.x * velMasB.w * invrhoPresMuBx * (dot(vel_XSPH_A - vel_XSPH_B, gradW)
	//			+ zeta * paramsD.HSML * (10 * paramsD.v_Max) * 2 * (rhoPresMuB.x / rhoPresMuA.x - 1) *
	// rAB_Dot_GradW_OverDist
	//			);

	// write new velocity back to original unsorted location
	//  uint originalIndex = gridMarkerIndex[index];

	// Real newDensity = densityShare + paramsD.markerMass * W3(0); //?$ include the particle in its summation as well
	// if (rhoPreMuA.w < -.1) { rhoPreMuA.x = newDensity; }
	// rhoPreMuA.y = Eos(rhoPreMuA.x, rhoPreMuA.w);
	//   oldRhoPreMu[originalIndex] = rhoPreMuA;
