 using namespace std;
 
 
-real_ AngleF3F3(real3 a, real3 b) {
-	return acos(dot(a, b));
+real_ AngleF3F3(real3 a3, real3 b3) {
+	real_ cosTheta = dot(a3, b3);
+	if(cosTheta>1){
+		cosTheta=1;
+	} else if(cosTheta<-1){
+		cosTheta=-1;
+	}
+	return acos(cosTheta);
 }
 //*******************************************************************************************************************************
 

	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;
		} else {
			fileRigidParticleCenterVsTimeAndDistance.open("dataRigidCenterVsTimeAndDistance.txt", ios::app);
		}

//////-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	ofstream flexTipPos;
	int numFlexBodiesInOnePeriod = int(ANCF_ReferenceArrayNodesOnBeams.size() / real_(paramsH.nPeriod) + .5);
	int tFlexTipPos = 1000;
	if (tStep % tFlexTipPos == 0) {
		if (tStep / tFlexTipPos == 0) {
			flexTipPos.open("dataFlexTip.txt");
