 	ofstream fileRigidParticleCenterVsTimeAndDistance;
 	int numRigidBodiesInOnePeriod = int(posRigidH.size() / real_(paramsH.nPeriod) + .5);
 	int tStepRigidCenterPos = 1000;
+
+
 	if (tStep % tStepRigidCenterPos == 0) {
 		if (tStep / tStepRigidCenterPos == 0) {
 			fileRigidParticleCenterVsTimeAndDistance.open("dataRigidCenterVsTimeAndDistance.txt");
-			fileRigidParticleCenterVsTimeAndDistance<<"(t, x, dist[or y], dum[or z],   x_cumul, y, z,   vx, vy, vz,   omega_x, omega_y, omega_z) (sequentially for all particles), average flow velocity (x, y, z, magnitude), channel radius\n" << endl;
+			fileRigidParticleCenterVsTimeAndDistance<<"(t, x, dist[or y], dum[or z],   x_cumul, y, z,   vx, vy, vz, thetaRigid_ZAxis_withX, thetaRigid_ZAxis_withY, thetaRigid_ZAxis_withZ, omega_x, omega_y, omega_z) (sequentially for all particles), average flow velocity (x, y, z, magnitude), channel radius\n" << endl;
 		} else {
 			fileRigidParticleCenterVsTimeAndDistance.open("dataRigidCenterVsTimeAndDistance.txt", ios::app);
 		}

//				ssParticleCenterVsTime << realTime << ", " << p_rigid.y << ", " << p_rigid.z;
//				ssParticleCenterVsDistance << p_rigidCumul.x << ", " << p_rigid.y << ", " << p_rigid.z;

//				//***radial distance (tube)
				real2 dist2 = R2(channelCenterYZ.x - p_rigid.y, channelCenterYZ.y - p_rigid.z);
				printf("center %f %f and radius %f and py and pz %f %f\n", channelCenterYZ.x, channelCenterYZ.y, channelRadius, p_rigid.y, p_rigid.z);

				int2 nodesPortioni2 = ANCF_ReferenceArrayNodesOnBeams[j];
				real3 p_tip = ANCF_NodesH[nodesPortioni2.y - 1];
				real3 v_tip = ANCF_NodesVelH[nodesPortioni2.y - 1];
				real3 s_tip = ANCF_SlopesH[nodesPortioni2.y - 1];

				real_ thetaX = AngleF3F3(s_tip, R3(1,0,0));
				real_ thetaY = AngleF3F3(s_tip, R3(0,1,0));
				real_ thetaZ = AngleF3F3(s_tip, R3(0,0,1));


				ssBeamTipVsTime << realTime << ", " <<
						p_tip.x << ", " << p_tip.y << ", " << p_tip.z << ", " <<
						v_tip.x << ", " << v_tip.y << ", " << v_tip.z << ", " <<
						thetaX << ", " << thetaY << ", " << thetaZ << ", " <<
						length(aveVel) << ", ";
			}
			ssBeamTipVsTime<<endl;
