Repository: FloatingArrayDesign/MoorDyn (language: python)

Issue:
Pinned bodies feature
This PR includes a new option for bodies, allowing them to behave as a pinned object. It contains a corresponding unit test and documentation updates. It also contains a feature update that prints time steps to the terminal while running, similar to how MDF works at runtime.

The formulation follows the pinned rods orientation closely. Inertial responses and forces have been added to the getStateDeriv and getFnet functions for both bodies and rods. The torque from translational acceleration is:
**F_iner [4:6] = M.bottomLeftCorner * a6[1:3]** 

and the inertial response force to acceleration is given as (contains two terms to avoid counting the externally applied force):
**F_iner[1:3] = - M.topLeftCorner * a6[1:3] - M.topRightCorner * a6[3:6]**

This feature has been verified in MDF using CAD software, analytical solutions, and Subdyn. For the MDC version I verified using MDF, and the two versions match very closely.

The corresponding update will be included in OpenFAST v4. 

This update did require the addition of an acceleration term for each calling time step for the calculation of inertial forces. In MDF acceleration a value passed into the step function. For MDC, I did this the same way velocity is handled: constant over the coupling step. While this isn't the most accurate, it was the simplest approach and in most use cases the coupling time step will be small enough that the effect is negligible (there was never issues when comparing with MDF). 


Discussion / hints:

Fix patch:
diff --git a/docs/inputs.rst b/docs/inputs.rst
index aa6a5a99..35fc2235 100644
--- a/docs/inputs.rst
+++ b/docs/inputs.rst
@@ -654,6 +654,10 @@ Footnotes:
 - The tension on the Line n anchor can be output with the ANCHTEN[n] flag (see examples above)
 - Object indicates output is for whole object, Node indicates output is for node of object
 - Coupled/fixed bodies and points will output acceleration 0 because no forces are calculated
+- When looking at the rotational outputs of coupled pinned rods that are hanging near vertical, 
+  the best approach is to attach a rod to a zero-mass, zero-volume pinned body and output the body 
+  rotations. Hanging pinned rods are inverted (end A over end B) in MoorDyn and the output range 
+  for roll/pitch of rods is +/- 180 degrees. 
 - There are a couple additional outputs left over from OpenFAST conventions that don’t follow the 
   same format: FairTen and AnchTen. FairTen[n] is the same as Line[n]TenB. For example, the 
   fairlead tension of line 1 would be FAIRTEN1 or LINE1TENB.
diff --git a/source/Body.cpp b/source/Body.cpp
index bcf52fb6..836c9b69 100644
--- a/source/Body.cpp
+++ b/source/Body.cpp
@@ -85,6 +85,16 @@ Body::setup(int number_in,
 		bodyI = vec::Zero();
 		bodyCdA = vec6::Zero();
 		bodyCa = vec6::Zero();
+	} else if (type == CPLDPIN){
+		bodyM = M_in;
+		bodyV = V_in;
+
+		body_r6.head<3>() = vec3::Zero();
+		body_r6.tail<3>() = deg2rad * r6_in.tail<3>();
+		body_rCG = rCG_in;
+		bodyI = I_in;
+		bodyCdA = CdA_in;
+		bodyCa = Ca_in;
     } else // coupled bodies have no need for these variables...
 	{
 		bodyM = 0.0;
@@ -158,13 +168,13 @@ Body::addRod(Rod* rod, vec6 coords)
 };
 
 void
-Body::initializeUnfreeBody(vec6 r6_in, vec6 v6_in)
+Body::initializeUnfreeBody(vec6 r6_in, vec6 v6_in, vec6 a6_in)
 {
 	if (type == FREE) {
 		LOGERR << "Invalid initializator for a FREE body" << endl;
 		throw moordyn::invalid_value_error("Invalid body type");
 	}
-	initiateStep(r6_in, v6_in);
+	initiateStep(r6_in, v6_in, a6_in);
 	updateFairlead(0.0);
 
 	// If any Rod is fixed to the body (not pinned), initialize it now because
@@ -181,7 +191,7 @@ Body::initializeUnfreeBody(vec6 r6_in, vec6 v6_in)
 std::pair<XYZQuat, vec6>
 Body::initialize()
 {
-	if (type != FREE) {
+	if ((type != FREE) && (type != CPLDPIN)) {
 		LOGERR << "Invalid initializator for a non FREE body ("
 		       << TypeName(type) << ")" << endl;
 		throw moordyn::invalid_value_error("Invalid body type");
@@ -191,15 +201,17 @@ Body::initialize()
 	// initialized)
 	setDependentStates();
 
-	// If any Rod is fixed to the body (not pinned), initialize it now because
-	// otherwise it won't be initialized
-	for (auto attached : attachedR)
-		if (attached->type == Rod::FIXED)
+	if (type == FREE) { // Attached objects already initalized for coupled pinned body
+		// If any Rod is fixed to the body (not pinned), initialize it now because
+		// otherwise it won't be initialized
+		for (auto attached : attachedR)
+			if (attached->type == Rod::FIXED)
+				attached->initialize();
+		// If there's an attached Point, initialize it now because it won't be
+		// initialized otherwise
+		for (auto attached : attachedP)
 			attached->initialize();
-	// If there's an attached Point, initialize it now because it won't be
-	// initialized otherwise
-	for (auto attached : attachedP)
-		attached->initialize();
+	}
 
 	// create output file for writing output (and write channel header and units
 	// lines) if applicable
@@ -313,7 +325,13 @@ real
 Body::GetBodyOutput(OutChanProps outChan)
 {
 
-	vec3 rotations = rad2deg * Quat2Euler(r7.quat);
+	vec3 rotations;
+	vec6 Fout;
+
+	if ((outChan.QType == RX)||(outChan.QType == RY)||(outChan.QType == RZ))
+		rotations = rad2deg * Quat2Euler(r7.quat);
+	if ((outChan.QType == FX)||(outChan.QType == FY)||(outChan.QType == FZ)||(outChan.QType == MX)||(outChan.QType == MY)||(outChan.QType == MZ))
+		Fout = getFnet();
 
 	if (outChan.QType == PosX)
 		return r7.pos.x();
@@ -355,17 +373,17 @@ Body::GetBodyOutput(OutChanProps outChan)
 		return sqrt(F6net[0] * F6net[0] + F6net[1] * F6net[1] +
 		            F6net[2] * F6net[2]);
 	else if (outChan.QType == FX)
-		return F6net[0];
+		return Fout[0];
 	else if (outChan.QType == FY)
-		return F6net[1];
+		return Fout[1];
 	else if (outChan.QType == FZ)
-		return F6net[2];
+		return Fout[2];
 	else if (outChan.QType == MX)
-		return F6net[3];
+		return Fout[3];
 	else if (outChan.QType == MY)
-		return F6net[4];
+		return Fout[4];
 	else if (outChan.QType == MZ)
-		return F6net[5];
+		return Fout[5];
 	else {
 		LOGWRN << "Unrecognized output channel " << outChan.QType << endl;
 		return 0.0;
@@ -375,12 +393,20 @@ Body::GetBodyOutput(OutChanProps outChan)
 // called at the beginning of each coupling step to update the boundary
 // conditions (body kinematics) for the proceeding time steps
 void
-Body::initiateStep(vec6 r, vec6 rd)
+Body::initiateStep(vec6 r, vec6 rd, vec6 rdd)
 {
 	if (type == COUPLED) // if coupled, update boundary conditions
 	{
 		r_ves = r;
 		rd_ves = rd;
+		rdd_ves = rdd;
+		return;
+	}
+	if (type == CPLDPIN) // if coupled, update boundary conditions
+	{
+		r_ves(Eigen::seqN(0, 3)) = r(Eigen::seqN(0, 3));
+		rd_ves(Eigen::seqN(0, 3)) = rd(Eigen::seqN(0, 3));
+		rdd_ves(Eigen::seqN(0, 3)) = rdd(Eigen::seqN(0, 3));
 		return;
 	}
 	if (type == FIXED) // if fixed body, set the BCs to stationary
@@ -388,6 +414,7 @@ Body::initiateStep(vec6 r, vec6 rd)
 		if (bodyId == 0) r_ves = vec6::Zero(); // special ground body case
 		else r_ves = r;
 		rd_ves = vec6::Zero();
+		rdd_ves = vec6::Zero();
 		return;
 	}
 	LOGERR << "Body " << number << "is not of type COUPLED or FIXED." << endl;
@@ -402,6 +429,7 @@ Body::updateFairlead(real time)
 		// set Body kinematics based on BCs (linear model for now)
 		r7 = XYZQuat::fromVec6(r_ves + rd_ves * time);
 		v6 = rd_ves;
+		a6 = rdd_ves;
 
 		// calculate orientation matrix based on latest angles
 		OrMat = r7.quat.toRotationMatrix();
@@ -411,6 +439,14 @@ Body::updateFairlead(real time)
 
 		return;
 	}
+	if (type == CPLDPIN) {
+		// set Body translational kinematics based on BCs (linear model for now)
+		r7.pos = r_ves(Eigen::seqN(0, 3)) + rd_ves(Eigen::seqN(0, 3)) * time;
+		v6(Eigen::seqN(0, 3)) = rd_ves(Eigen::seqN(0, 3));
+		a6(Eigen::seqN(0, 3)) = rdd_ves(Eigen::seqN(0, 3));
+
+		return;
+	}
 	LOGERR << "The body is not a coupled/fixed one" << endl;
 	throw moordyn::invalid_value_error("Invalid body type");
 }
@@ -418,9 +454,17 @@ Body::updateFairlead(real time)
 void
 Body::setState(XYZQuat pos, vec6 vel)
 {
-	// set position and velocity vectors from state vector
-	r7 = pos;
-	v6 = vel;
+	if (type == FREE) {
+		// set position and velocity vectors from state vector
+		r7 = pos;
+		v6 = vel;
+	} else if (type == CPLDPIN) {
+		r7.quat = pos.quat;
+		v6(Eigen::seqN(3, 3)) = vel(Eigen::seqN(3, 3));
+	} else {
+		LOGERR << "Invalid body type: " << TypeName(type) << endl;
+		throw moordyn::invalid_value_error("Invalid body type");
+	}
 
 	// calculate orientation matrix based on latest angles
 	OrMat = r7.quat.toRotationMatrix();
@@ -432,26 +476,65 @@ Body::setState(XYZQuat pos, vec6 vel)
 std::pair<XYZQuat, vec6>
 Body::getStateDeriv()
 {
-	if (type != FREE) {
-		LOGERR << "The body is not a free one" << endl;
+	if ((type != FREE) && (type != CPLDPIN)) {
+		LOGERR << "getStateDeriv called for non-free body" << endl;
 		throw moordyn::invalid_value_error("Invalid body type");
 	}
 
 	// Get contributions from attached points (and lines attached to
 	// them) and store in FNet:
 	doRHS();
-
-	// solve for accelerations in [M]{a}={f}
-	a6 = solveMat6(M, F6net);
-
-	// NOTE; is the above still valid even though it includes rotational DOFs?
-	dPos.pos = v6.head<3>();
-	// this assumes that the angular velocity is about the global coordinates
-	// which is true for bodies
-	dPos.quat = 0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
+	if (type == FREE){
+		// solve for accelerations in [M]{a}={f}
+		a6 = solveMat6(M, F6net);
+
+		// NOTE; is the above still valid even though it includes rotational DOFs?
+		dPos.pos = v6.head<3>();
+		// this assumes that the angular velocity is about the global coordinates
+		// which is true for bodies
+		dPos.quat = 0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
+	} else { // Pinned body
+
+		// account for moment in response to acceleration due to inertial coupling (off-diagonal sub-matrix terms)
+		const vec Fnet_out3 = F6net(Eigen::seqN(3, 3)) - (M.bottomLeftCorner<3,3>() * a6(Eigen::seqN(0, 3)));
+
+		// For small systems it is usually faster to compute the inverse
+		// of the matrix. See
+		// https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
+		const mat M_out3 = M(Eigen::seqN(3, 3), Eigen::seqN(3, 3));
+		a6(Eigen::seqN(3, 3)) = M_out3.inverse() * Fnet_out3;
+
+		// dxdt = V   (velocities)
+		dPos.pos = vec::Zero();
+		dPos.quat = 0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
+	}
 	return std::make_pair(dPos, a6);
 };
 
+vec6
+Body::getFnet()
+{
+	vec6 F6_iner = vec6::Zero();
+	vec6 Fnet_out = vec6::Zero();
+
+	// this assumes doRHS() has already been called
+	if (type == COUPLED) { // if coupled rigidly
+		F6_iner = - M * a6; // Inertial terms
+		Fnet_out = F6net + F6_iner;
+	} else if (type == CPLDPIN) { // if coupled pinned
+		F6_iner(Eigen::seqN(0,3)) = (- M.topLeftCorner<3,3>() * a6(Eigen::seqN(0,3))) - (M.topRightCorner<3,3>() * a6(Eigen::seqN(3,3))); // Inertial term
+		Fnet_out(Eigen::seqN(0,3)) = F6net(Eigen::seqN(0,3)) + F6_iner(Eigen::seqN(0,3));
+		Fnet_out(Eigen::seqN(3,3)) = vec3::Zero();
+	} else {
+		// LOGERR << "Invalid body type: " << TypeName(type) << endl;
+		// throw moordyn::invalid_value_error("Invalid body type");
+
+		Fnet_out = F6net;
+	}
+
+	return Fnet_out;
+}
+
 //  this is the big function that calculates the forces on the body
 void
 Body::doRHS()
diff --git a/source/Body.hpp b/source/Body.hpp
index c7140811..914d600c 100644
--- a/source/Body.hpp
+++ b/source/Body.hpp
@@ -134,6 +134,8 @@ class Body final : public io::IO
 	vec6 r_ves;
 	/// fairlead velocity for coupled bodies [x/y/z]
 	vec6 rd_ves;
+	/// fairlead acceleration for coupled bodies [x/y/z]
+	vec6 rdd_ves;
 
 	/// total force and moment vector on body
 	vec6 F6net;
@@ -161,6 +163,8 @@ class Body final : public io::IO
 		FREE = 0,
 		/// Is fixed, either to a location or to another moving entity
 		FIXED = 1,
+		/// Is coupled pinned, i.e. translational dof are controlled by the user
+		CPLDPIN = 2,
 		// Some aliases
 		VESSEL = COUPLED,
 		ANCHOR = FIXED,
@@ -175,6 +179,8 @@ class Body final : public io::IO
 		switch (t) {
 			case COUPLED:
 				return "COUPLED";
+			case CPLDPIN:
+				return "CPLDPIN";	
 			case FREE:
 				return "FREE";
 			case FIXED:
@@ -237,10 +243,11 @@ class Body final : public io::IO
 	 * moordyn::Body::FIXED
 	 * @param r The position (6 dof)
 	 * @param rd The velocity (6 dof)
+	 * @param rdd The acceleration (6 dof)
 	 * @throw moordyn::invalid_value_error If the body is of type
 	 * moordyn::Body::FREE
 	 */
-	void initializeUnfreeBody(vec6 r = vec6::Zero(), vec6 rd = vec6::Zero());
+	void initializeUnfreeBody(vec6 r = vec6::Zero(), vec6 rd = vec6::Zero(), vec6 rdd = vec6::Zero());
 
 	/** @brief Initialize the FREE point state
 	 * @return The 6-dof position (first) and the 6-dof velocity (second)
@@ -317,7 +324,7 @@ class Body final : public io::IO
 	/** @brief Get the forces and moments exerted over the body
 	 * @return The net force
 	 */
-	inline vec6 getFnet() const { return F6net; }
+	vec6 getFnet();
 
 	/** @brief Get the mass and intertia matrix
 	 * @return The mass and inertia matrix
@@ -339,12 +346,13 @@ class Body final : public io::IO
 
 	/** @brief Called at the beginning of each coupling step to update the
 	 * boundary conditions (body kinematics) for the proceeding time steps
-	 * @param r The output position
-	 * @param rd The output velocity
+	 * @param r The input position
+	 * @param rd The input velocity
+	 * @param rdd The input velocity
 	 * @throw moordyn::invalid_value_error If the body is not of type
 	 * moordyn::Body::COUPLED or moordyn::Body::FIXED
 	 */
-	void initiateStep(vec6 r, vec6 rd);
+	void initiateStep(vec6 r_in, vec6 rd_in, vec6 rdd_in);
 
 	/** @brief Sets the kinematics based on the position and velocity of the
 	 * fairlead.
diff --git a/source/Line.hpp b/source/Line.hpp
index e278ebb0..93c1e445 100644
--- a/source/Line.hpp
+++ b/source/Line.hpp
@@ -409,6 +409,14 @@ class Line final : public io::IO
 			       << ", which only has " << N + 1 << " nodes" << std::endl;
 			throw moordyn::invalid_value_error("Invalid node index");
 		}
+		if (isnan(r[i].sum())) {
+			stringstream s;
+			s << "NaN detected" << endl
+			  << "Line " << number << " node positions:" << endl;
+			for (unsigned int j = 0; j <= N; j++)
+				s << j << " : " << r[j] << ";" << endl;
+			throw moordyn::nan_error(s.str().c_str());
+		}
 		return r[i];
 	}
 
diff --git a/source/Misc.hpp b/source/Misc.hpp
index 1fab4ea5..8d725e19 100644
--- a/source/Misc.hpp
+++ b/source/Misc.hpp
@@ -672,7 +672,9 @@ split(const string& str, const char sep);
 inline vector<string>
 split(const string& s)
 {
-	return split(s, ' ');
+	vector<string> sout = split(s, ' ');
+	if (sout.size() == 1) return split(sout[0], '	');
+	else return sout;
 }
 
 /**
diff --git a/source/MoorDyn2.cpp b/source/MoorDyn2.cpp
index 152ca6d3..ac36ffa6 100644
--- a/source/MoorDyn2.cpp
+++ b/source/MoorDyn2.cpp
@@ -214,35 +214,62 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic)
 	int ix = 0;
 
 	for (auto l : CpldBodyIs) {
-		LOGMSG << "Initializing coupled Body " << l << " in " << x[ix] << ", "
+		LOGMSG << "Initializing coupled Body " << l+1 << " at " << x[ix] << ", "
 		       << x[ix + 1] << ", " << x[ix + 2] << "..." << endl;
-		// this calls initiateStep and updateFairlead, then initializes
-		// dependent Rods
+
 		// BUG: These conversions will not be needed in the future
-		vec6 r, rd;
-		moordyn::array2vec6(x + ix, r);
-		moordyn::array2vec6(xd + ix, rd);
-		BodyList[l]->initializeUnfreeBody(r, rd);
-		ix += 6;
+		vec6 r, rd, rdd;
+		if (BodyList[l]->type == Body::COUPLED){
+			moordyn::array2vec6(x + ix, r);
+			moordyn::array2vec6(xd + ix, rd);
+			// determine acceleration 
+			rdd = (rd - rd_b) / dtM0;
+			rd_b = rd;
+
+			ix += 6;
+		} else {
+			// for pinned body 3 entries will be taken
+			vec3 r3, rd3, rdd3;
+			moordyn::array2vec(x + ix, r3);
+			r(Eigen::seqN(0, 3)) = r3;
+			moordyn::array2vec(xd + ix, rd3);
+			rd(Eigen::seqN(0, 3)) = rd3;
+			// determine acceleration 
+		    rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_b) / dtM0;
+			rd3_b = rd3;
+
+			ix += 3;
+		}
+		BodyList[l]->initializeUnfreeBody(r, rd, rdd);
 	}
 
 	for (auto l : CpldRodIs) {
-		LOGMSG << "Initializing coupled Rod " << l << " in " << x[ix] << ", "
+		LOGMSG << "Initializing coupled Rod " << l+1 << " at " << x[ix] << ", "
 		       << x[ix + 1] << ", " << x[ix + 2] << "..." << endl;
-		vec6 r, rd;
+		vec6 r, rd, rdd;
 		if (RodList[l]->type == Rod::COUPLED) {
+			// for cantilevered rods 6 entries will be taken
 			moordyn::array2vec6(x + ix, r);
 			moordyn::array2vec6(xd + ix, rd);
-			ix += 6; // for cantilevered rods 6 entries will be taken
+			// determine acceleration 
+		    rdd = (rd - rd_r) / dtM0;
+			rd_r = rd;
+
+			ix += 6;
 		} else {
-			vec3 r3, rd3;
+			// for pinned rods 3 entries will be taken
+			vec3 r3, rd3, rdd3;
 			moordyn::array2vec(x + ix, r3);
 			r(Eigen::seqN(0, 3)) = r3;
 			moordyn::array2vec(xd + ix, rd3);
 			rd(Eigen::seqN(0, 3)) = rd3;
-			ix += 3; // for pinned rods 3 entries will be taken
+			// determine acceleration 
+		    rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_r) / dtM0;
+			rd3_r = rd3;
+
+			ix += 3;
 		}
-		RodList[l]->initiateStep(r, rd);
+		RodList[l]->initiateStep(r, rd, rdd);
 		RodList[l]->updateFairlead(0.0);
 		// call this just to set up the output file header
 		RodList[l]->initialize();
@@ -250,7 +277,7 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic)
 
 	for (auto l : CpldPointIs) {
 		LOGMSG << "Initializing coupled Point " << l+1 << " at " << x[ix] << ", "
-		       << x[ix + 1] << ", " << x[ix + 2] << "..." << endl;
+		       << x[ix + 1] << ", " << x[ix + 2] << endl;
 		vec r, rd;
 		moordyn::array2vec(x + ix, r);
 		moordyn::array2vec(xd + ix, rd);
@@ -485,6 +512,9 @@ moordyn::MoorDyn::Step(const double* x,
 	// should check if wave kinematics have been set up if expected!
 	LOGDBG << "t = " << t << "s     \r";
 
+	cout << fixed << setprecision(1);
+	cout << "\rt = " << t << " " << flush;
+
 	if (dt <= 0) {
 		// Nothing to do, just recover the forces if there are coupled DOFs
 		if (NCoupledDOF())
@@ -505,29 +535,55 @@ moordyn::MoorDyn::Step(const double* x,
 	// used later for extrapolating motions
 	for (auto l : CpldBodyIs) {
 		// BUG: These conversions will not be needed in the future
-		vec6 r, rd;
-		moordyn::array2vec6(x + ix, r);
-		moordyn::array2vec6(xd + ix, rd);
-		BodyList[l]->initiateStep(r, rd);
-		ix += 6;
+		vec6 r, rd, rdd;
+		if (BodyList[l]->type == Body::COUPLED){
+			moordyn::array2vec6(x + ix, r);
+			moordyn::array2vec6(xd + ix, rd);
+			// determine acceleration 
+			rdd = (rd - rd_b) / dtM0;
+			rd_b = rd;
+
+			ix += 6;
+		} else {
+			// for pinned body 3 entries will be taken
+			vec3 r3, rd3, rdd3;
+			moordyn::array2vec(x + ix, r3);
+			r(Eigen::seqN(0, 3)) = r3;
+			moordyn::array2vec(xd + ix, rd3);
+			rd(Eigen::seqN(0, 3)) = rd3;
+			// determine acceleration 
+		    rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_b) / dt;
+			rd3_b = rd3;
+
+			ix += 3;
+		}
+		BodyList[l]->initiateStep(r, rd, rdd); // acceleration required for inertial terms
 	}
 	for (auto l : CpldRodIs) {
-		vec6 r, rd;
+		vec6 r, rd, rdd;
 		if (RodList[l]->type == Rod::COUPLED) {
 			// for cantilevered rods 6 entries will be taken
 			moordyn::array2vec6(x + ix, r);
 			moordyn::array2vec6(xd + ix, rd);
+			// determine acceleration 
+		    rdd = (rd - rd_r) / dt;
+			rd_r = rd;
+
 			ix += 6;
 		} else {
 			// for pinned rods 3 entries will be taken
-			vec3 r3, rd3;
+			vec3 r3, rd3, rdd3;
 			moordyn::array2vec(x + ix, r3);
 			r(Eigen::seqN(0, 3)) = r3;
 			moordyn::array2vec(xd + ix, rd3);
 			rd(Eigen::seqN(0, 3)) = rd3;
+			// determine acceleration 
+		    rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_r) / dt;
+			rd3_r = rd3;
+
 			ix += 3;
 		}
-		RodList[l]->initiateStep(r, rd);
+		RodList[l]->initiateStep(r, rd, rdd); // acceleration required for inertial terms
 	}
 	for (auto l : CpldPointIs) {
 		vec r, rd;
@@ -1344,7 +1400,7 @@ moordyn::MoorDyn::ReadInFile()
 					} else if (let3 == "AX") {
 						dummy.QType = AccX;
 						dummy.Units = moordyn::UnitList[AccX];
-					} else if (let3 == "Ay") {
+					} else if (let3 == "AY") {
 						dummy.QType = AccY;
 						dummy.Units = moordyn::UnitList[AccY];
 					} else if (let3 == "AZ") {
@@ -1745,6 +1801,16 @@ moordyn::MoorDyn::readBody(string inputText)
 		// it is coupled - controlled from outside
 		type = Body::COUPLED;
 		CpldBodyIs.push_back(ui_size(BodyList));
+	} else if (str::isOneOf(let1, { "VESSELPINNED", "VESPIN", "CPLDPIN", "COUPLEDPINNED" })) {
+		// if a pinned fairlead, add to list and add
+		type = Body::CPLDPIN;
+		CpldBodyIs.push_back(
+		    ui_size(BodyList)); // index of fairlead in BodyList vector
+		FreeBodyIs.push_back(
+		    ui_size(BodyList));    // also add this pinned body to the free
+		                          // list because it is half free
+		BodyStateIs.push_back(nX); // assign start index of this body's states
+		nX += 6;                  // add 6 state variables for each pinned Body
 	} else {
 		// it is free - controlled by MoorDyn
 		type = Body::FREE;
@@ -1836,7 +1902,7 @@ moordyn::MoorDyn::readRod(string inputText)
 		type = Rod::COUPLED;
 		CpldRodIs.push_back(
 		    ui_size(RodList)); // index of fairlead in RodList vector
-	} else if (str::isOneOf(let1, { "VESPIN", "CPLDPIN" })) {
+	} else if (str::isOneOf(let1, { "VESSELPINNED", "VESPIN", "CPLDPIN", "COUPLEDPINNED" })) {
 		// if a pinned fairlead, add to list and add
 		type = Rod::CPLDPIN;
 		CpldRodIs.push_back(
diff --git a/source/MoorDyn2.hpp b/source/MoorDyn2.hpp
index 55290e44..ff4f249c 100644
--- a/source/MoorDyn2.hpp
+++ b/source/MoorDyn2.hpp
@@ -388,8 +388,13 @@ class MoorDyn final : public io::IO
 		for (auto l : CpldBodyIs) {
 			// BUG: These conversions will not be needed in the future
 			const vec6 f_i = BodyList[l]->getFnet();
-			moordyn::vec62array(f_i, f + ix);
-			ix += 6;
+			if (BodyList[l]->type == Body::COUPLED) {
+				moordyn::vec62array(f_i, f + ix);
+				ix += 6; // for coupled bodies 6 entries will be taken
+			} else {
+				moordyn::vec2array(f_i(Eigen::seqN(0, 3)), f + ix);
+				ix += 3; // for pinned bodies 3 entries will be taken
+			}
 		}
 		for (auto l : CpldRodIs) {
 			const vec6 f_i = RodList[l]->getFnet();
@@ -502,6 +507,12 @@ class MoorDyn final : public io::IO
 	/// (if using env.WaveKin=1)
 	unsigned int npW;
 
+	/// Previous time step velocity for calculating coupled point acceleration (first time step assumed stationary)
+	vec6 rd_b = vec6::Zero(); // body
+	vec6 rd_r = vec6::Zero(); // coupled rod
+	vec3 rd3_r = vec3::Zero(); // coupled pinned rod
+	vec3 rd3_b = vec3::Zero(); // coupled pinned body
+
 	/// main output file
 	ofstream outfileMain;
 
diff --git a/source/Rod.cpp b/source/Rod.cpp
index fc352245..24afae24 100644
--- a/source/Rod.cpp
+++ b/source/Rod.cpp
@@ -148,12 +148,13 @@ Rod::setup(int number_in,
 	if (type == FREE) {
 		// For an independent rod, set the position right off the bat
 		r7.pos = endCoords.head<3>();
-		r7.quat = quaternion::Identity();
+		r7.quat = quaternion::FromTwoVectors(q0, endCoords.tail<3>());
 		v6 = vec6::Zero();
 	} else if ((type == PINNED) || (type == CPLDPIN)) {
 		// for a pinned rod, just set the orientation (position will be set
-		// later by parent object)
-		r7 = XYZQuat::Zero();
+		// later by parent object) 
+		r7.pos = vec3::Zero(); 
+		r7.quat = quaternion::FromTwoVectors(q0, endCoords.tail<3>()); // TODO: Check this is right
 		v6(Eigen::seqN(3, 3)) = vec::Zero();
 	}
 	// otherwise (for a fixed rod) the positions will be set by the parent body
@@ -337,6 +338,12 @@ Rod::initialize()
 real
 Rod::GetRodOutput(OutChanProps outChan)
 {
+	
+	vec6 Fout;
+	if ((outChan.QType == FX)||(outChan.QType == FY)||(outChan.QType == FZ)||(outChan.QType == MX)||(outChan.QType == MY)||(outChan.QType == MZ))
+		Fout = getFnet();
+
+
 	if (outChan.NodeID == -1) {
 		if (outChan.QType == PosX)
 			return r7.pos.x();
@@ -373,17 +380,17 @@ Rod::GetRodOutput(OutChanProps outChan)
 		else if (outChan.QType == TenB)
 			return FextB.norm();
 		else if (outChan.QType == FX)
-			return F6net[0];
+			return Fout[0];
 		else if (outChan.QType == FY)
-			return F6net[1];
+			return Fout[1];
 		else if (outChan.QType == FZ)
-			return F6net[2];
+			return Fout[2];
 		else if (outChan.QType == MX)
-			return F6net[3];
+			return Fout[3];
 		else if (outChan.QType == MY)
-			return F6net[4];
+			return Fout[4];
 		else if (outChan.QType == MZ)
-			return F6net[5];
+			return Fout[5];
 		else if (outChan.QType == Sub) {
 			real VOFsum = 0.0;
 			for (unsigned int i = 0; i <= N; i++)
@@ -446,20 +453,22 @@ Rod::setState(XYZQuat pos, vec6 vel)
 }
 
 void
-Rod::initiateStep(vec6 rFairIn, vec6 rdFairIn)
+Rod::initiateStep(vec6 r_in, vec6 rd_in, vec6 rdd_in)
 {
 	if (type == COUPLED) {
 		// set Rod kinematics based on BCs (linear model for now)
-		r_ves = rFairIn;
-		rd_ves = rdFairIn;
+		r_ves = r_in;
+		rd_ves = rd_in;
+		rdd_ves = rdd_in;
 
 		// since this rod has no states and all DOFs have been set, pass its
 		// kinematics to dependent Lines
 		setDependentStates();
 	} else if (type == CPLDPIN) {
 		// set Rod *end A only* kinematics based on BCs (linear model for now)
-		r_ves(Eigen::seqN(0, 3)) = rFairIn(Eigen::seqN(0, 3));
-		rd_ves(Eigen::seqN(0, 3)) = rdFairIn(Eigen::seqN(0, 3));
+		r_ves(Eigen::seqN(0, 3)) = r_in(Eigen::seqN(0, 3));
+		rd_ves(Eigen::seqN(0, 3)) = rd_in(Eigen::seqN(0, 3));
+		rdd_ves(Eigen::seqN(0, 3)) = rdd_in(Eigen::seqN(0, 3));
 	} else {
 		LOGERR << "Invalid rod type: " << TypeName(type) << endl;
 		throw moordyn::invalid_value_error("Invalid rod type");
@@ -475,6 +484,7 @@ Rod::updateFairlead(real time)
 		// set Rod kinematics based on BCs (linear model for now)
 		r7 = XYZQuat::fromVec6(r_ves + rd_ves * time);
 		v6 = rd_ves;
+		acc6 = rdd_ves;
 
 		// since this rod has no states and all DOFs have been set, pass its
 		// kinematics to dependent Lines
@@ -483,6 +493,7 @@ Rod::updateFairlead(real time)
 		// set Rod *end A only* kinematics based on BCs (linear model for now)
 		r7.pos = r_ves(Eigen::seqN(0, 3)) + rd_ves(Eigen::seqN(0, 3)) * time;
 		v6(Eigen::seqN(0, 3)) = rd_ves(Eigen::seqN(0, 3));
+		acc6(Eigen::seqN(0, 3)) = rdd_ves(Eigen::seqN(0, 3));
 
 		// Rod is pinned so only end A is specified, rotations are left alone
 		// and will be handled, along with passing kinematics to dependent
@@ -505,7 +516,7 @@ Rod::setKinematics(vec6 r_in, vec6 rd_in)
 		// since this rod has no states and all DOFs have been set, pass its
 		// kinematics to dependent Lines
 		setDependentStates();
-	} else if (type == PINNED) // rod end A pinned to a body
+	} else if ((type == PINNED) || (type == CPLDPIN))// rod end A pinned to a body
 	{
 		// set Rod *end A only* kinematics based on BCs (linear model for now)
 		r7.pos = r_in(Eigen::seqN(0, 3));
@@ -650,22 +661,21 @@ Rod::getStateDeriv()
 			acc6 = solveMat6(M_out6, Fnet_out);
 		}
 	} else {
-		// Pinned rod, where the position rate of change is null (Dirichlet BC)
+		// Pinned rod
+
+		// account for moment in response to end A acceleration due to inertial coupling (off-diagonal sub-matrix terms)
+		const vec Fnet_out3 = Fnet_out(Eigen::seqN(3, 3)) - (M_out6.bottomLeftCorner<3,3>() * acc6(Eigen::seqN(0, 3)));
 
 		// For small systems it is usually faster to compute the inverse
 		// of the matrix. See
 		// https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
-		const vec Fnet_out3 = Fnet_out(Eigen::seqN(3, 3));
 		const mat M_out3 = M_out6(Eigen::seqN(3, 3), Eigen::seqN(3, 3));
-		const vec acc = M_out3.inverse() * Fnet_out3;
+		acc6(Eigen::seqN(3, 3)) = M_out3.inverse() * Fnet_out3;
 
 		// dxdt = V   (velocities)
 		vel7.pos = vec::Zero();
 		vel7.quat =
 		    0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
-		// dVdt = a   (accelerations)
-		acc6(Eigen::seqN(0, 3)) = vec::Zero();
-		acc6(Eigen::seqN(3, 3)) = acc;
 	}
 
 	return std::make_pair(vel7, acc6);
@@ -674,48 +684,50 @@ Rod::getStateDeriv()
 vec6
 Rod::getFnet()
 {
-	return F6net;
-	// >>>>>>>>>>>>> do I want to leverage getNetForceAndMass or a saved global
-	// to save comp time and code?  >>>>>>>> this function likely not used
-	// >>>>>>>>>>>
-
-	// Fnet_out is assumed to point to a size-3 array if the rod is pinned and
-	// size-6 array if the rod is fixed
+	// return F6net;
+	// // >>>>>>>>>>>>> do I want to leverage getNetForceAndMass or a saved global
+	// // to save comp time and code?  >>>>>>>> this function likely not used
+	// // >>>>>>>>>>>
 
-	int nDOF = 0;
+	// // Fnet_out is assumed to point to a size-3 array if the rod is pinned and
+	// // size-6 array if the rod is fixed
 
-	if (type == CPLDPIN) // if coupled pinned
-		nDOF = 3;
-	else if (type == COUPLED) // if coupled rigidly
-		nDOF = 6;
-	else {
-		LOGERR << "Invalid rod type: " << TypeName(type) << endl;
-		throw moordyn::invalid_value_error("Invalid rod type");
-	}
+	vec6 F6_iner = vec6::Zero();
+	vec6 Fnet_out = vec6::Zero();
 
 	// this assumes doRHS() has already been called
+	if (type == COUPLED) { // if coupled rigidly
+		F6_iner = - M6net * acc6; // Inertial terms
+		Fnet_out = F6net + F6_iner;
+	} else if (type == CPLDPIN) { // if coupled pinned
+		F6_iner(Eigen::seqN(0,3)) = (- M6net.topLeftCorner<3,3>() * acc6(Eigen::seqN(0,3))) - (M6net.topRightCorner<3,3>() * acc6(Eigen::seqN(3,3))); // Inertial term
+		Fnet_out(Eigen::seqN(0,3)) = F6net(Eigen::seqN(0,3)) + F6_iner(Eigen::seqN(0,3));
+		Fnet_out(Eigen::seqN(3,3)) = vec3::Zero();
+	} else {
+		// LOGERR << "Invalid rod type: " << TypeName(type) << endl;
+		// throw moordyn::invalid_value_error("Invalid rod type");
 
-	// make sure Fnet_out is zeroed first
-	vec6 Fnet_out = vec6::Zero();
+		Fnet_out = F6net;
+	}
 
-	// now go through each node's contributions, put them in body ref frame, and
-	// sum them
-	for (unsigned int i = 0; i <= N; i++) {
-		// position of a given node relative to the body reference point (global
-		// orientation frame)
-		vec rRel = r[i] - r[0];
+	// // now go through each node's contributions, put them in body ref frame, and
+	// // sum them
+	// for (unsigned int i = 0; i <= N; i++) {
+	// 	// position of a given node relative to the body reference point (global
+	// 	// orientation frame)
+	// 	vec rRel = r[i] - r[0];
 
-		vec6 F6_i; // 6dof force-moment from rod about body ref point (but
-		           // global orientation frame of course)
-		// convert segment net force into 6dof force about body ref point
-		F6_i.head<3>() = Fnet[i];
-		F6_i.tail<3>() = rRel.cross(Fnet[i]);
+	// 	vec6 F6_i; // 6dof force-moment from rod about body ref point (but
+	// 	           // global orientation frame of course)
+	// 	// convert segment net force into 6dof force about body ref point
+	// 	F6_i.head<3>() = Fnet[i];
+	// 	F6_i.tail<3>() = rRel.cross(Fnet[i]);
 
-		Fnet_out(Eigen::seqN(0, nDOF)) += F6_i(Eigen::seqN(0, nDOF));
-	}
+	// 	Fnet_out(Eigen::seqN(0, nDOF)) += F6_i(Eigen::seqN(0, nDOF));
+	// }
 
-	// add any moments applied from lines at either end (might be zero)
-	Fnet_out(Eigen::seqN(3, 3)) += Mext;
+	// // add any moments applied from lines at either end (might be zero)
+	// Fnet_out(Eigen::seqN(3, 3)) += Mext;
 
 	//  this is where we'd add inertial loads for coupled rods! <<<<<<<<<<<<
 
diff --git a/source/Rod.hpp b/source/Rod.hpp
index 8223cbc4..c33e9a49 100644
--- a/source/Rod.hpp
+++ b/source/Rod.hpp
@@ -203,6 +203,8 @@ class Rod final : public io::IO
 	vec6 r_ves;
 	/// fairlead velocity for coupled rods [x/y/z]
 	vec6 rd_ves;
+	/// fairlead acceleration for coupled rods [x/y/z]
+	vec6 rdd_ves;
 
 	// file stuff
 	/// Pointer to moordyn::MoorDyn::outfileMain
@@ -363,6 +365,14 @@ class Rod final : public io::IO
 			       << ", which only has " << N + 1 << " nodes" << std::endl;
 			throw moordyn::invalid_value_error("Invalid node index");
 		}
+		if (isnan(r[i].sum())) {
+			stringstream s;
+			s << "NaN detected" << endl
+			  << "Rod " << number << " node positions:" << endl;
+			for (unsigned int j = 0; j <= N; j++)
+				s << j << " : " << r[j] << ";" << endl;
+			throw moordyn::nan_error(s.str().c_str());
+		}
 		return r[i];
 	}
 
@@ -442,14 +452,15 @@ class Rod final : public io::IO
 
 	/** @brief Called at the beginning of each coupling step to update the
 	 * boundary conditions (rod kinematics) for the proceeding time steps
-	 * @param r The output position
-	 * @param rd The output velocity
+	 * @param r The input position
+	 * @param rd The input velocity
+	 * @param rdd The input velocity
 	 * @throw moordyn::invalid_value_error If the rod is not of type
 	 * moordyn::Rod::COUPLED or moordyn::Rod::CPLDPIN
 	 * @note If the rod is of type moordyn::Rod::CPLDPIN, then just 3 components
 	 * of @p r and @p rd are considered
 	 */
-	void initiateStep(vec6 r, vec6 rd);
+	void initiateStep(vec6 r_in, vec6 rd_in, vec6 rdd_in);
 
 	/** @brief Sets the kinematics
 	 *
diff --git a/source/Time.cpp b/source/Time.cpp
index acffbc19..12aa4094 100644
--- a/source/Time.cpp
+++ b/source/Time.cpp
@@ -44,7 +44,7 @@ TimeSchemeBase<NSTATE, NDERIV>::Update(real t_local, unsigned int substep)
 
 	t_local += this->t_local;
 	for (auto obj : bodies) {
-		if (obj->type != Body::COUPLED)
+		if ((obj->type != Body::COUPLED) && (obj->type != Body::CPLDPIN))
 			continue;
 		obj->updateFairlead(t_local);
 	}
@@ -63,7 +63,7 @@ TimeSchemeBase<NSTATE, NDERIV>::Update(real t_local, unsigned int substep)
 	}
 
 	for (unsigned int i = 0; i < bodies.size(); i++) {
-		if (bodies[i]->type != Body::FREE)
+		if ((bodies[i]->type != Body::FREE) && (bodies[i]->type != Body::CPLDPIN))
 			continue;
 		bodies[i]->setState(r[substep].bodies[i].pos, r[substep].bodies[i].vel);
 	}
@@ -115,7 +115,7 @@ TimeSchemeBase<NSTATE, NDERIV>::CalcStateDeriv(unsigned int substep)
 	}
 
 	for (unsigned int i = 0; i < bodies.size(); i++) {
-		if (bodies[i]->type != Body::FREE)
+		if ((bodies[i]->type != Body::FREE) && (bodies[i]->type != Body::CPLDPIN))
 			continue;
 		std::tie(rd[substep].bodies[i].vel, rd[substep].bodies[i].acc) =
 		    bodies[i]->getStateDeriv();
@@ -127,7 +127,7 @@ TimeSchemeBase<NSTATE, NDERIV>::CalcStateDeriv(unsigned int substep)
 		obj->doRHS();
 	}
 	for (auto obj : rods) {
-		if ((obj->type != Rod::COUPLED) && (obj->type != Rod::CPLDPIN))
+		if (obj->type != Rod::COUPLED)
 			continue;
 		obj->doRHS();
 	}
diff --git a/source/Time.hpp b/source/Time.hpp
index 70446c9f..f8f2018e 100644
--- a/source/Time.hpp
+++ b/source/Time.hpp
@@ -494,7 +494,7 @@ class TimeSchemeBase : public TimeScheme
 		// integrator, no matter if they are free or not. Thus they can change
 		// types (mutate) without needing to micromanage them in the time scheme
 		for (unsigned int i = 0; i < bodies.size(); i++) {
-			if (bodies[i]->type != Body::FREE)
+			if ((bodies[i]->type != Body::FREE) && (bodies[i]->type != Body::CPLDPIN)) // Only fully coupled bodies are intialized in MD2.cpp
 				continue;
 			std::tie(r[0].bodies[i].pos, r[0].bodies[i].vel) =
 			    bodies[i]->initialize();


Test patch:
diff --git a/tests/Mooring/body_tests/pinnedBody.txt b/tests/Mooring/body_tests/pinnedBody.txt
new file mode 100644
index 00000000..3eb0cdbe
--- /dev/null
+++ b/tests/Mooring/body_tests/pinnedBody.txt
@@ -0,0 +1,43 @@
+The points and line are dummy objects just to allow MD to run w/o error
+------------------------- LINE TYPES --------------------------------------------------
+LineType  Diam    MassDenInAir    EA       BA/-zeta    EI      Can   Cat  Cdn   Cdt 
+(-)       (m)       (kg/m)        (N)      (Pa-s/-)  (n-m^2)   (-)   (-)  (-)   (-)  
+nylon      0.124      35.0    5515288.0     -0.8      0.0     1.0   0.0  1.6   0.05
+---------------------- ROD TYPES ------------------------------------
+TypeName      Diam     Mass/m    Cd     Ca      CdEnd    CaEnd
+(name)        (m)      (kg/m)    (-)    (-)     (-)      (-)
+Can          2       5e3     0.6    1.0     1.2      1.0
+----------------------- BODIES ------------------------------------------------------
+ID   Attachment         X0      Y0     Z0     r0      p0     y0     Mass          CG*          I*      Volume   CdA*   Ca*
+(#)     (-)             (m)     (m)    (m)   (deg)   (deg)  (deg)   (kg)          (m)         (kg-m^2)  (m^3)   (m^2)  (-)
+1     CoupledPinned    0.00    0.00  -5.00   0.00   0.00   0.00   0.00        0.00        0.000e+00  0.00     0.00  0.00
+---------------------- RODS ----------------------------------------
+ID   RodType  Attachment     Xa    Ya    Za    Xb    Yb     Zb    NumSegs  RodOutputs
+(#)  (name)    (#/key)      (m)   (m)   (m)   (m)   (m)    (m)   (-)       (-)
+1      Can     Body1         0     0     -5     0     0     -15     3         p
+----------------------- POINTS ----------------------------------------------
+Node      Type      X        Y         Z        M        V         CdA   CA
+(-)       (-)      (m)      (m)       (m)      (kg)     (m^3)     (m^2)  (-)
+1         Fixed    30    0      -5       0        0          0     0
+2         Fixed    30    0      -10       0        0          0     0
+-------------------------- LINES -------------------------------------------------
+Line     LineType NodeA     NodeB  UnstrLen  NumSegs     Flags/Outputs
+(-)      (-)       (-)       (-)   (m)         (-)          (-)
+1        nylon      1         2    5        1           p
+-------------------------- SOLVER OPTIONS---------------------------------------------------
+2        writeLog     - Write a log file
+0.0     WtrDnsty     - Water Density
+0.001   dtM          - time step to use in mooring integration
+3.0e6    kb           - bottom stiffness
+3.0e5    cb           - bottom damping
+30       WtrDpth      - water depth
+0.015    threshIC     - threshold for IC convergence
+0.0    TmaxIC       - threshold for IC convergence
+0.01     dtIC         - Time lapse between convergence tests (s)
+------------------------------ OUTPUTS -----------------------------------------------------
+BODY1RX
+BODY1RY
+BODY1PX
+BODY1PY
+BODY1AY
+--------------------------- need this line -------------------------------------------------
\ No newline at end of file
diff --git a/tests/bodies.cpp b/tests/bodies.cpp
index c116a75d..a02e1fa6 100644
--- a/tests/bodies.cpp
+++ b/tests/bodies.cpp
@@ -42,6 +42,8 @@
 #include <algorithm>
 #include <cmath>
 #include <limits>
+#include <vector>
+#include <numeric>
 
 #include "util.h"
 
@@ -415,6 +417,122 @@ rotatingBody(SeriesWriter* series_writer)
 
 	return true;
 }
+
+/**
+ * @brief Compares inertial deflection of a pinned body to analytical solution
+ *
+ * 
+ * The coupled pinned body that is massless and volumeless with a rod fixed to it
+ * is moved with constant acceleration in a vaccum (0 water density). 
+ * The resulting avg inertial deflection should match an analytical solution of 
+ * theta = arctan(-accel/gravity)
+ * 
+ * This solution was derived both with Newtonian and LaGrangian mechanics. 
+ * It is the same as the pendulum on an accelerating cart problem. 
+ * 
+ * This only tests the inertial properties of pinned bodies, other tests deal 
+ * with hydrodynamics and general body properties
+ *
+ * 
+ * @param series_writer
+ * @return true
+ * @return false
+ */
+bool
+pinnedBody(SeriesWriter* series_writer)
+{
+	int err;
+	cout << endl << " => " << __PRETTY_FUNC_NAME__ << "..." << endl;
+
+	MoorDyn system = MoorDyn_Create("Mooring/body_tests/pinnedBody.txt");
+	if (!system) {
+		cerr << "Failure Creating the Mooring system" << endl;
+		return false;
+	}
+
+	unsigned int n_dof;
+	err = MoorDyn_NCoupledDOF(system, &n_dof);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure getting NCoupledDOF: " << err << endl;
+		return false;
+	}
+	if (n_dof != 6) { // rotational DOF are ignored by MDC, same as a coupled pinned rods
+		cerr << "Expected 6 DOFs but got " << n_dof << endl;
+		return false;
+	}
+
+	moordyn::vec6 x{ 0, 0, -5, 0, 0, 0};
+	moordyn::vec6 xd{ 0, 0, 0 , 0, 0, 0};
+	double f[6];
+
+	err = MoorDyn_Init(system, x.data(), xd.data());
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure during the mooring initialization: " << err << endl;
+		return false;
+	}
+
+	if (!write_system_vtk(system, 0, series_writer)) {
+		return false;
+	}
+
+	auto body = MoorDyn_GetBody(system, 1);
+	moordyn::vec6 r, rd;
+	vector<double> roll;
+	int i = 0, j = 0;
+	double t = 0.0, dt = 0.01, accel = 0.5;
+	bool local_min_max;
+
+	if (!write_system_vtk(system, t, series_writer)) {
+		return false;
+	}
+
+	while (t < 50.0) {
+
+		x[1] = 0.5*accel*pow(t, 2);
+		xd[1] = accel*t;
+		err = MoorDyn_Step(system, x.data(), xd.data(), f, &t, &dt);
+		if (err != MOORDYN_SUCCESS) {
+			cerr << "Failure during the mooring initialization: " << err
+			     << endl;
+			return false;
+		}
+		if (!write_system_vtk(system, t, series_writer)) {
+			return false;
+		}
+
+		MoorDyn_GetBodyState(body, r.data(), rd.data());
+		roll.push_back(r[3]);
+
+		if (i >= 30) { // after the simulation has run for a few time steps
+			// When local min or max of oscillation, indicates half of an oscialltion has occured
+			local_min_max = (((roll[i]-roll[i-1])/dt) * ((roll[i-1]-roll[i-2])/dt)) < 0; 
+			if (local_min_max) j++;
+		}
+		if (j > 3) break; // after 2 full oscillations
+
+		t = t + dt;
+		i++;
+	}
+	double theta = atan(-accel/9.80665);
+	double average = reduce(roll.begin(), roll.end()) / roll.size();
+	if (abs(average - theta) > 0.001) {
+		cerr << "Pinned body inertial deflection should be "
+		<< theta << " but it is " << average << endl;
+		return false;
+	}
+
+	err = MoorDyn_Close(system);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure closing Moordyn: " << err << endl;
+		return false;
+	}
+
+	cout << setprecision(4) << "Average roll is " << average << endl;
+	cout << setprecision(4) << "Theoretical roll is " << theta << endl;
+
+	return true;
+}
+
 /** @brief Runs all the test
  * @return 0 if the tests have ran just fine. The index of the failing test
  * otherwise
@@ -435,6 +553,19 @@ main(int, char**)
 		return 3;
 	}
 
+	try {
+		// SeriesWriter series_writer;
+		if (!pinnedBody(NULL)) {
+			// series_writer.writeJson("../../vtk_out/");
+			return 3;
+		}
+		// series_writer.writeJson("../../vtk_out/");
+
+	} catch (std::exception& e) {
+		cerr << "pinnedBody failed with exception " << e.what() << endl;
+		return 3;
+	}
+
 	cout << "bodies.cpp passed successfully" << endl;
 	return 0;
 }