 	int indexOfClosestNode = int(s / l) * nNodes;
 	if (indexOfClosestNode == nNodes) indexOfClosestNode--;
 
-	real3 beamPointPos = ANCF_Point_Pos(ANCF_Nodes, ANCF_Slopes, indexOfClosestNode, s, l); //interpolation using ANCF beam, cubic hermit equation
+	real3 beamPointPos = Calc_ANCF_Point_Pos(ANCF_Nodes, ANCF_Slopes, indexOfClosestNode, s, l); //interpolation using ANCF beam, cubic hermit equation
 
 	real3 dist3 = posRadD[absMarkerIndex] - beamPointPos;
 	flexSPH_MeshPos_LRF_D[index] = dist3;
+}
+//--------------------------------------------------------------------------------------------------------------------------------
+__global__ void Populate_FlexSPH_MeshSlope_LRF_kernel(
+		real3* flexSPH_MeshSlope_Initial_D,
+		int* flexIdentifierD,
+		real_* parametricDist,
+		real_* ANCF_Beam_Length,
+		int* ANCF_NumNodes_Per_Beam,
+		real3 * ANCF_Nodes,
+		real3 * ANCF_Slopes) {
+	uint index = blockIdx.x * blockDim.x + threadIdx.x;
+	if (index >= numFlex_SphMarkersD) {
+		return;
+	}
+	uint absMarkerIndex = index + startFlexMarkersD; // updatePortionD = [start, end] index of the update portion
+	real_ s = parametricDist[index];
+	int flexBodyIndex = flexIdentifierD[index];
+	real_ l = ANCF_Beam_Length[flexBodyIndex];
+	int nNodes = ANCF_NumNodes_Per_Beam[flexBodyIndex];
+
+	int indexOfClosestNode = int(s / l) * nNodes;
+	if (indexOfClosestNode == nNodes) indexOfClosestNode--;
+
+	real3 beamPointSlope = Calc_ANCF_Point_Slope(ANCF_Nodes, ANCF_Slopes, indexOfClosestNode, s, l); //interpolation using ANCF beam, cubic hermit equation
+	flexSPH_MeshSlope_Initial_D[index] = beamPointSlope;
 }
 
 //--------------------------------------------------------------------------------------------------------------------------------

	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
// in wikipedia, last quat comp is the angle, in my version, first one is the angle.
__global__ void RotationMatirixFromQuaternion_kernel(real3 * AD1, real3 * AD2, real3 * AD3, real4 * qD) {
	uint rigidSphereA = blockIdx.x * blockDim.x + threadIdx.x;
	if (rigidSphereA >= numRigidBodiesD) {


		//******************************************************************************
	thrust::device_vector<real3> flexSPH_MeshPos_LRF_D(numFlex_SphMarkers);
	thrust::device_vector<real3> flexSPH_MeshSlope_Initial_D(numFlex_SphMarkers);
	uint nBlocks_numFlex_SphMarkers;
	uint nThreads_SphMarkers;
	computeGridSize(numFlex_SphMarkers, 256, nBlocks_numFlex_SphMarkers, nThreads_SphMarkers);
	printf("before first kernel\n");

	Populate_FlexSPH_MeshPos_LRF_kernel<<<nBlocks_numFlex_SphMarkers, nThreads_SphMarkers>>>(R3CAST(flexSPH_MeshPos_LRF_D), R3CAST(posRadD), I1CAST(flexIdentifierD), R1CAST(parametricDist), R1CAST(ANCF_Beam_Length),
			I1CAST(ANCF_NumNodes_Per_Beam), R3CAST(ANCF_Nodes), R3CAST(ANCF_Slopes));
	cudaThreadSynchronize();
		CUT_CHECK_ERROR("Kernel execution failed: Populate_FlexSPH_MeshPos_LRF_kernel");	printf("after first kernel\n");

	Populate_FlexSPH_MeshSlope_LRF_kernel<<<nBlocks_numFlex_SphMarkers, nThreads_SphMarkers>>>(R3CAST(flexSPH_MeshSlope_Initial_D), I1CAST(flexIdentifierD), R1CAST(parametricDist), R1CAST(ANCF_Beam_Length),
				I1CAST(ANCF_NumNodes_Per_Beam), R3CAST(ANCF_Nodes), R3CAST(ANCF_Slopes));
	cudaThreadSynchronize();
	CUT_CHECK_ERROR("Kernel execution failed: Populate_FlexSPH_MeshSlope_LRF_kernel");	printf("after first kernel\n");

	//******************************************************************************
	thrust::device_vector<real4> qD1 = mQuatRot;
