 // 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) {
+	if (rigidSphereA >= numRigidBodiesD) {
 		return;
 	}
 	float3 omega = omegaLRF_D[rigidSphereA];

	FILE *outFileMultipleZones;

	int povRayCounter = 0;
	int stepEnd =2.4e6 * (.02 * sizeScale) / delT ; //1.4e6 * (.02 * sizeScale) / delT ;//0.7e6 * (.02 * sizeScale) / delT ;//0.7e6;//2.5e6; //200000;//10000;//50000;//100000;
	printf("stepEnd %d\n", stepEnd);

	float delTOrig = delT;

//local torque = T' = A' * T
__global__ void MapTorqueToLRFKernel(float3 * AD1, float3 * AD2, float3 * AD3, float3 * totalTorque3, float3 * LF_totalTorque3) {
	uint rigidSphereA = blockIdx.x * blockDim.x + threadIdx.x;
	if (rigidSphereA >= numRigidBodiesD) {
		return;
	}
	float3 totalTorqueGRF = totalTorque3[rigidSphereA];
