 
 	omega3 += omegaDot3 * dTD;
 	omegaLRF_D[rigidSphereA] = omega3;
-}
-//--------------------------------------------------------------------------------------------------------------------------------
-//updates the rigid body particles
-__global__ void UpdateKernelRigid_XZ_Motion(
-									 float4 * totalForces4,
-									 float4 * totalTorque4,
-									 float4 * posRadRigidD,
-									 float4 * velMassRigidD,
-									 float4 * cylinderRotOmegaJD,
-									 float3 * deltaPositionD,
-									 float3 * deltaVelocityD,
-									 float3 * deltaXZTetaOmegaD,
-									 float rigid_SPH_mass) {
-	uint rigidSphereA = blockIdx.x * blockDim.x + threadIdx.x;
-	if (rigidSphereA > numRigidBodiesD) {return;}
-
-	float4 dummyPosRad = posRadRigidD[rigidSphereA];
-	float4 dummyVelMas = velMassRigidD[rigidSphereA];
-	float4 dummyTetaOmegaJ = cylinderRotOmegaJD[rigidSphereA];
-
-	float4 totalForce4 = totalForces4[rigidSphereA];
-	totalForce4.y = 0;		//the only different line from UpdateKernelRigid
-	float3 derivV_SPH = rigid_SPH_mass * F3(totalForce4) / dummyVelMas.w;		//in fact, totalForce4 is originially sum of dV/dt of sph particles and should be multiplied by m to produce force. gravity is applied in the force kernel
-
-	float3 deltaPos = F3(dummyVelMas) * dTD;
-	deltaPos.y = 0;
-		//deltaPos = F3(0);//
-	dummyPosRad += F4(deltaPos, 0);
-
-	float3 deltaVel = derivV_SPH * dTD;
-	deltaVel.y = 0;
-		//deltaVel = F3(0);//
-	dummyVelMas += F4(deltaVel, 0);
-//------------- torque
-	float totalTorque = totalTorque4[rigidSphereA].y;					//for the XZ motion, only y component of the torque is necessary
-	float derivOmega = rigid_SPH_mass * totalTorque / dummyTetaOmegaJ.z;
-	float deltaTeta = dummyTetaOmegaJ.y * dTD;
-		//deltaTeta = 0;
-	float deltaOmega = derivOmega * dTD;
-		//deltaOmega = 0;
-	dummyTetaOmegaJ += F4(deltaTeta, deltaOmega, 0, 0);
-	dummyTetaOmegaJ.x -= (dummyTetaOmegaJ.x > 2 * PI) ? 2 * PI : 0;
-	dummyTetaOmegaJ.x += (dummyTetaOmegaJ.x < 0) ? 2 * PI : 0;
-
-	//float3 boundaryDeltaPos = F3(0);
-	//float3 boundaryDeltaVel = F3(0);
-	//ApplyBoundaryRigid2(dummyPosRad, dummyVelMas, boundaryDeltaPos, boundaryDeltaVel);
-	//deltaPos += boundaryDeltaPos;
-	//deltaVel += boundaryDeltaVel;
-
-	deltaPositionD[rigidSphereA] = deltaPos;
-	deltaVelocityD[rigidSphereA] = deltaVel;
-	deltaXZTetaOmegaD[rigidSphereA] = F3(deltaTeta, deltaOmega, 0);
-
-	posRadRigidD[rigidSphereA] = dummyPosRad;
-	velMassRigidD[rigidSphereA] = dummyVelMas;
-	cylinderRotOmegaJD[rigidSphereA] = dummyTetaOmegaJ;
 }
 //--------------------------------------------------------------------------------------------------------------------------------
 //updates the rigid body particles

		posRadRigidD2.clear();
		velMassRigidD2.clear();
		vel_XSPH_D.clear();
		AD1_2.clear();
		AD2_2.clear();
		AD3_2.clear();

		UpdateFluid(posRadD, velMasD, vel_XSPH_D, rhoPresMuD, derivVelRhoD, referenceArray, delT);
		//UpdateBoundary(posRadD, velMasD, rhoPresMuD, derivVelRhoD, referenceArray, delT);
		UpdateRigidBody(posRadD, velMasD,
				posRadRigidD, velMassRigidD, qD1, AD1, AD2, AD3, omegaLRF_D,
				derivVelRhoD, rigidIdentifierD, rigidSPH_MeshPos_LRF_D, referenceArray, jD1, jD2, jInvD1, jInvD2, paramsH, delT);
		ApplyBoundary(posRadD, rhoPresMuD, mNSpheres, posRadRigidD, velMassRigidD, numRigidBodies);

