 	cudaThreadSynchronize();
 	CUT_CHECK_ERROR("Kernel execution failed: CalcTorqueShare");
 	//******************************************************************************
+	thrust::device_vector<float4> qD1(numRigidBodies);
 	thrust::device_vector<float3> AD1(numRigidBodies);
 	thrust::device_vector<float3> AD2(numRigidBodies);
 	thrust::device_vector<float3> AD3(numRigidBodies);
+	thrust::fill(qD1.begin(), qD1.end(), F4(1, 0, 0, 0));
 	thrust::fill(AD1.begin(), AD1.end(), F3(1, 0, 0));
 	thrust::fill(AD2.begin(), AD2.end(), F3(0, 1, 0));
 	thrust::fill(AD3.begin(), AD3.end(), F3(0, 0, 1));

					 const thrust::device_vector<int> & rigidIdentifierD,
					 const thrust::device_vector<float3> & rigidSPH_MeshPos_LRF_D,
					 const thrust::host_vector<int3> & referenceArray,

					 const thrust::device_vector<float3> & jD1,
					 const thrust::device_vector<float3> & jD2,
					 const thrust::device_vector<float3> & jInvD1,

	velMassRigidD[rigidSphereA] = dummyVelMas;
}
//--------------------------------------------------------------------------------------------------------------------------------
//updates the rigid body Quaternion of Rotation
// A is rotation matrix, A = [AD1; AD2; AD3]
__global__  void UpdateRigidBodyQuaternion_kernel(
									float4 * qD,
									float3 * omegaLRF_D) {
	uint rigidSphereA = blockIdx.x * blockDim.x + threadIdx.x;
	if (rigidSphereA > numRigidBodiesD) {return;}
	float3 omega = omegaLRF_D[rigidSphereA];
	float4 q = qD[rigidSphereA];
	float4 qDot = omega.x * F4(-(q.y), q.x, q.w, -(q.z)) + omega.y * F4(-(q.z), -(q.w), q.x, q.y) + omega.z * F4(-(q.w), q.z, -(q.y), q.x);
	q += .5 * dTD * qDot;
	q *= (1 / length(q));
	qD[rigidSphereA] = q;
}
//--------------------------------------------------------------------------------------------------------------------------------
//updates the rigid body Rotation
// A is rotation matrix, A = [AD1; AD2; AD3], first comp of q is rotation, last 3 components are axis of rot
__global__  void RotationMatirixFromQuaternion_kernel(
									float3 * AD1,
									float3 * AD2,
									float3 * AD3,
									float4 * qD) {
	uint rigidSphereA = blockIdx.x * blockDim.x + threadIdx.x;
	if (rigidSphereA > numRigidBodiesD) {return;}
	float4 q = qD[rigidSphereA];
	AD1[rigidSphereA] = 2 * F3(0.5f - q.z * q.z - q.w * q.w, q.y * q.z - q.x * q.w, q.y * q.w + q.x * q.z);
	AD2[rigidSphereA] = 2 * F3(q.y * q.z + q.x * q.w, 0.5f - q.y * q.y - q.w * q.w,  q.z * q.w - q.x * q.y);
	AD3[rigidSphereA] = 2 * F3(q.y * q.w - q.x * q.z,  q.z * q.w + q.x * q.y, 0.5f - q.y * q.y - q.z * q.z);
}
//--------------------------------------------------------------------------------------------------------------------------------
__global__ void UpdateRigidBodyAngularVelocity_kernel (
