diff --git a/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C b/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C
index e543f3c2a..f098e0ed4 100644
--- a/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C
+++ b/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C
@@ -51,6 +51,9 @@ int main(int argc, char *argv[])
     sixDOFBodies structure(runTime);
     OFstream of(runTime.path()/"motion.dat");
 
+    // Write header for output file
+    of << "# Time, CoG_x, CoG_y, CoG_z, omega_x, omega_y, omega_z" << endl;
+
     Info<< "\nStarting time loop\n" << endl;
 
     for (runTime++; !runTime.end(); runTime++)
@@ -73,11 +76,14 @@ int main(int argc, char *argv[])
                 << structure()[bodyI].omegaAverage().value() << nl
                 << "Current omega in time instant = "
                 << structure()[bodyI].omega().value()  << nl
-                << "Average omegaAbs in time step = "
-                << structure()[bodyI].omegaAverageAbsolute().value() << nl
                 << endl;
 
-            of  << structure()[bodyI].X().value().x() << tab;
+            of  << structure()[bodyI].X().value().x() << tab
+                << structure()[bodyI].X().value().y() << tab
+                << structure()[bodyI].X().value().z() << tab
+                << structure()[bodyI].omega().value().x() << tab
+                << structure()[bodyI].omega().value().y() << tab
+                << structure()[bodyI].omega().value().z() << tab;
         }
 
         of << endl;
diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.C b/src/ODE/sixDOF/finiteRotation/finiteRotation.C
index df12bdb4c..aa407d91f 100644
--- a/src/ODE/sixDOF/finiteRotation/finiteRotation.C
+++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.C
@@ -40,9 +40,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
 {
     vector ur = - *( inv(I + rotT) & (I - rotT) );
 
-    // Scaling to a unit vector.  HJ, problems with round-off
-    // HJ, 4/Aug/2008
-
+    // Scaling to a unit vector. Problems with round-off. HJ, 4/Aug/2008
     if (mag(ur) > SMALL)
     {
         return ur/(mag(ur) + SMALL);
@@ -50,8 +48,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
     else
     {
         // Rotation vector is undertermined at zero rotation
-        // Returning arbitrary unit vector
-        // HJ, 4/Mar/2015
+        // Returning arbitrary unit vector. HJ, 4/Mar/2015
         return vector(0, 0, 1);
     }
 }
diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
index 507c1d69b..f70ad4908 100644
--- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
+++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
@@ -170,7 +170,7 @@ Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
 {
     tensor R;
 
-    // Calculate the magnitude of the rotation increment vector
+    // Calculate the magnitude of the rotational increment vector
     const scalar magRotInc = mag(rotInc);
 
     if (magRotInc < rotIncTensorTol_)
@@ -202,7 +202,7 @@ Foam::vector Foam::geometricSixDOF::dexpMap
 {
     vector rotIncDot;
 
-    // Calculate the magnitude of the rotation increment vector
+    // Calculate the magnitude of the rotational increment vector
     const scalar magRotInc = mag(rotInc);
 
     if (magRotInc < rotIncRateTol_)
@@ -241,11 +241,14 @@ Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
 
     Xrel_(dict().lookup("Xrel")),
     U_(dict().lookup("U")),
-    Uaverage_(U_),
+    Uaverage_("Uaverage", U_),
     rotation_(tensor(dict().lookup("rotationTensor"))),
+    rotIncrement_
+    (
+        dict().lookupOrDefault<tensor>("rotationIncrementTensor", tensor::zero)
+    ),
     omega_(dict().lookup("omega")),
-    omegaAverage_(omega_),
-    omegaAverageAbsolute_(omega_),
+    omegaAverage_("omegaAverage", omega_),
 
     nEqns_(),
     coeffs_(),
@@ -326,11 +329,6 @@ Foam::geometricSixDOF::geometricSixDOF
     rotation_(gsd.rotation_),
     omega_(gsd.omega_.name(), gsd.omega_),
     omegaAverage_(gsd.omegaAverage_.name(), gsd.omegaAverage_),
-    omegaAverageAbsolute_
-    (
-        gsd.omegaAverageAbsolute_.name(),
-        gsd.omegaAverageAbsolute_
-    ),
 
     nEqns_(gsd.nEqns_),
     coeffs_(gsd.coeffs_),
@@ -401,80 +399,6 @@ const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
 }
 
 
-const Foam::finiteRotation& Foam::geometricSixDOF::rotation() const
-{
-    return rotation_;
-}
-
-
-Foam::vector Foam::geometricSixDOF::rotVector() const
-{
-    return rotation_.rotVector();
-}
-
-
-Foam::dimensionedScalar Foam::geometricSixDOF::rotAngle() const
-{
-    return dimensionedScalar("rotAngle", dimless, rotation_.rotAngle());
-}
-
-
-void Foam::geometricSixDOF::setXrel(const vector& x)
-{
-    Xrel_.value() = x;
-
-    // Set X in coefficients
-    coeffs_[0] = x.x();
-    coeffs_[1] = x.y();
-    coeffs_[2] = x.z();
-}
-
-
-void Foam::geometricSixDOF::setU(const vector& U)
-{
-    U_.value() = U;
-    Uaverage_ = U_;
-
-    // Set U in coefficients
-    coeffs_[3] = U.x();
-    coeffs_[4] = U.y();
-    coeffs_[5] = U.z();
-}
-
-
-void Foam::geometricSixDOF::setRotation(const HamiltonRodriguezRot& rot)
-{
-    // Set increment rotation vector to zero
-    coeffs_[9] = 0;
-    coeffs_[10] = 0;
-    coeffs_[11] = 0;
-}
-
-
-void Foam::geometricSixDOF::setOmega(const vector& omega)
-{
-    omega_.value() = omega;
-    omegaAverage_ = omega_;
-    omegaAverageAbsolute_ = omega_;
-
-    // Set omega in coefficients
-    coeffs_[6] = omega.x();
-    coeffs_[7] = omega.y();
-    coeffs_[8] = omega.z();
-}
-
-
-void Foam::geometricSixDOF::setReferentRotation
-(
-    const HamiltonRodriguezRot& rot
-)
-{
-    referentRotation_ = rot;
-
-    referentMotionConstraints_ = true;
-}
-
-
 void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
 {
     // First set the state in base class subobject
@@ -490,7 +414,6 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
     rotation_ = gsd.rotation_;
     omega_ = gsd.omega_;
     omegaAverage_ = gsd.omegaAverage_;
-    omegaAverageAbsolute_ = gsd.omegaAverageAbsolute_;
 
     // Copy ODE coefficients: this carries actual ODE state
     // HJ, 23/Mar/2015
@@ -507,44 +430,27 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
 }
 
 
-Foam::vector Foam::geometricSixDOF::rotVectorAverage() const
-{
-    return rotation_.rotVectorAverage();
-}
-
-
 const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const
 {
     return omegaAverage_;
 }
 
 
-const Foam::dimensionedVector&
-Foam::geometricSixDOF::omegaAverageAbsolute() const
-{
-    return omegaAverageAbsolute_;
-}
-
-
 Foam::tensor Foam::geometricSixDOF::toRelative() const
 {
-    // Note: using rotation tensor directly since we are integrating rotational
-    // increment vector
-    return rotation_.rotTensor();
+    return rotation_;
 }
 
 
 Foam::tensor Foam::geometricSixDOF::toAbsolute() const
 {
-    // Note: using rotation tensor directly since we are integrating rotational
-    // increment vector
-    return rotation_.rotTensor().T();
+    return rotation_.T();
 }
 
 
 const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const
 {
-    return rotation_.rotIncrementTensor();
+    return rotIncrement_;
 }
 
 
@@ -569,7 +475,7 @@ void Foam::geometricSixDOF::derivatives
     const vector rotIncrementVector(y[9], y[10], y[11]);
 
     // Calculate current rotation tensor obtained with exponential map
-    const tensor curRot = (expMap(rotIncrementVector) & rotation_.rotTensor());
+    const tensor curRot = (expMap(rotIncrementVector) & rotation_);
 
     // Calculate translational acceleration using current rotation
     const vector accel = A(curX, curU, curRot).value();
@@ -612,7 +518,7 @@ void Foam::geometricSixDOF::update(const scalar delta)
     // Translation
 
     // Update displacement
-    vector Xold = Xrel_.value();
+    const vector Xold = Xrel_.value();
 
     vector& Xval = Xrel_.value();
 
@@ -638,27 +544,28 @@ void Foam::geometricSixDOF::update(const scalar delta)
     // Rotation
 
     // Update omega
+    const vector omegaOld = omega_.value();
+
     vector& omegaVal = omega_.value();
 
     omegaVal.x() = coeffs_[6];
     omegaVal.y() = coeffs_[7];
     omegaVal.z() = coeffs_[8];
 
-    // Update rotation with final increment vector
-    rotation_.updateRotationWithIncrement
-    (
-        expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]))
-    );
+    // Update rotational increment tensor
+    rotIncrement_ = expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]));
+
+    // Update rotational tensor
+    rotation_ = (rotIncrement_ & rotation_);
 
     // Reset increment vector in coefficients for the next step
     coeffs_[9] = 0;
     coeffs_[10] = 0;
     coeffs_[11] = 0;
 
-    omegaAverage_.value() = rotation_.omegaAverage(delta);
-
-    // Calculate omegaAverageAbsolute
-    omegaAverageAbsolute_.value() = rotation_.omegaAverageAbsolute(delta);
+    // Consider calculating average omega using rotational tensor and rotational
+    // increment tensors
+    omegaAverage_.value() = 0.5*(omegaVal + omegaOld);
 }
 
 
@@ -671,26 +578,28 @@ bool Foam::geometricSixDOF::writeData(Ostream& os) const
     os.writeKeyword("type") << tab << type() << token::END_STATEMENT << endl;
 
     // Write data
-    os.writeKeyword("Xrel") << tab << this->Xrel()
+    os.writeKeyword("Xrel") << tab << Xrel_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl;
-    os.writeKeyword("rotationTensor") << tab << this->toRelative()
+    os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationTensor") << tab << rotation_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("omega") << tab << this->omega()
+    os.writeKeyword("rotationIncrementTensor") << tab << rotIncrement_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("omega") << tab << omega_
         << token::END_STATEMENT << nl << nl;
 
-    os.writeKeyword("fixedSurge") << tab << this->fixedSurge_ <<
-        token::END_STATEMENT << endl;
-    os.writeKeyword("fixedSway") << tab << this->fixedSway_ <<
-        token::END_STATEMENT << endl;
-    os.writeKeyword("fixedHeave") << tab << this->fixedHeave_ <<
-        token::END_STATEMENT << endl;
-    os.writeKeyword("fixedRoll") << tab << this->fixedRoll_ <<
-        token::END_STATEMENT << endl;
-    os.writeKeyword("fixedPitch") << tab << this->fixedPitch_ <<
-        token::END_STATEMENT << endl;
-    os.writeKeyword("fixedYaw") << tab << this->fixedYaw_ <<
-        token::END_STATEMENT << endl;
+    os.writeKeyword("fixedSurge") << tab << fixedSurge_ <<
+        token::END_STATEMENT << nl;
+    os.writeKeyword("fixedSway") << tab << fixedSway_ <<
+        token::END_STATEMENT << nl;
+    os.writeKeyword("fixedHeave") << tab << fixedHeave_ <<
+        token::END_STATEMENT << nl;
+    os.writeKeyword("fixedRoll") << tab << fixedRoll_ <<
+        token::END_STATEMENT << nl;
+    os.writeKeyword("fixedPitch") << tab << fixedPitch_ <<
+        token::END_STATEMENT << nl;
+    os.writeKeyword("fixedYaw") << tab << fixedYaw_ <<
+        token::END_STATEMENT << nl << endl;
 
     return os.good();
 }
diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
index 18777f38a..64a26e9e5 100644
--- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
+++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
@@ -79,24 +79,22 @@ class geometricSixDOF
             //- Velocity of mass centroid
             dimensionedVector U_;
 
-            //- Average velocity of mass centroid at previous time-step
+            //- Average velocity of mass centroid (evaluated at midstep)
             dimensionedVector Uaverage_;
 
-            //- Finite rotation. Note: changed during solution process
-            mutable finiteRotation rotation_;
+            //- Rotation tensor
+            tensor rotation_;
+
+            //- Rotational increment tensor
+            tensor rotIncrement_;
 
             //- Rotational velocity about mass centroid
             dimensionedVector omega_;
 
-
-        // Average variables that need to be stored
-
             //- Average rotational velocity in relative coordinate system
+            //  (evaluated at midstep)
             dimensionedVector omegaAverage_;
 
-            //- Average rotational velocity in absolute coordinate system
-            dimensionedVector omegaAverageAbsolute_;
-
 
         // ODE controls
 
@@ -107,7 +105,7 @@ class geometricSixDOF
             scalarField coeffs_;
 
 
-        //- Motion constraints (given as fixed motion components)
+        // Motion constraints (given as fixed motion components)
 
             //- Fixed surge (x-translation)
             Switch fixedSurge_;
@@ -236,7 +234,7 @@ public:
                 virtual const dimensionedVector& omega() const;
 
 
-            // Displacement and rotation in the absolute coordinate system
+            // Displacement and velocity in the absolute coordinate system
 
                 //- Return position of origin in absolute coordinate system
                 virtual dimensionedVector X() const;
@@ -244,56 +242,22 @@ public:
                 //- Return velocity of origin
                 virtual const dimensionedVector& U() const;
 
-                //- Return average velocity of origin for the previous time-step
+                //- Return average velocity of origin (evaluated at midstep)
                 virtual const dimensionedVector& Uaverage() const;
 
-                //- Return finite rotation
-                virtual const finiteRotation& rotation() const;
-
-                //- Return rotational vector of body
-                virtual vector rotVector() const;
-
-                //- Return rotation angle of body
-                virtual dimensionedScalar rotAngle() const;
-
 
             // Non-access control for setting state variables
 
-                //- Set position of origin
-                virtual void setXrel(const vector& x);
-
-                //- Set velocity of origin
-                virtual void setU(const vector& u);
-
-                //- Set rotational angle in relative coordinate system
-                virtual void setRotation(const HamiltonRodriguezRot& rot);
-
-                //- Set rotational velocity in relative coordinate system
-                virtual void setOmega(const vector& omega);
-
-                //- Set referent coordinate system to apply constraints
-                virtual void setReferentRotation
-                (
-                    const HamiltonRodriguezRot& rot
-                );
-
                 //- Set ODE parameters from another ODE
                 virtual void setState(const sixDOFODE&);
 
 
             // Average motion per time-step
 
-                //- Return average rotational vector of body
-                virtual vector rotVectorAverage() const;
-
                 //- Return average rotational velocity in relative coordinate
-                //  system
+                //  system (evaluated at midstep)
                 virtual const dimensionedVector& omegaAverage() const;
 
-                //- Return average rotational velocity in absolute coordinate
-                //  system
-                virtual const dimensionedVector& omegaAverageAbsolute() const;
-
 
             // Rotations
 
@@ -360,6 +324,7 @@ public:
                 //- Update ODE after the solution, advancing by delta
                 virtual void update(const scalar delta);
 
+
         // Write controls
 
             //- WriteData member function required by regIOobject
diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C
index f414f7aea..45143c562 100644
--- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C
+++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C
@@ -198,15 +198,15 @@ Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
 
     Xrel_(dict().lookup("Xrel")),
     U_(dict().lookup("U")),
-    Uaverage_(U_),
+    Uaverage_("Uaverage", U_),
     rotation_
     (
         vector(dict().lookup("rotationVector")),
         dimensionedScalar(dict().lookup("rotationAngle")).value()
     ),
     omega_(dict().lookup("omega")),
-    omegaAverage_(omega_),
-    omegaAverageAbsolute_(omega_),
+    omegaAverage_("omegaAverage", omega_),
+    omegaAverageAbsolute_("omegaAverageAbsolute", omega_),
 
     coeffs_(13, 0.0),
 
@@ -347,89 +347,13 @@ const Foam::dimensionedVector& Foam::quaternionSixDOF::U() const
     return U_;
 }
 
+
 const Foam::dimensionedVector& Foam::quaternionSixDOF::Uaverage() const
 {
     return Uaverage_;
 }
 
 
-const Foam::finiteRotation& Foam::quaternionSixDOF::rotation() const
-{
-    return rotation_;
-}
-
-
-Foam::vector Foam::quaternionSixDOF::rotVector() const
-{
-    return rotation_.rotVector();
-}
-
-
-Foam::dimensionedScalar Foam::quaternionSixDOF::rotAngle() const
-{
-    return dimensionedScalar("rotAngle", dimless, rotation_.rotAngle());
-}
-
-
-void Foam::quaternionSixDOF::setXrel(const vector& x)
-{
-    Xrel_.value() = x;
-
-    // Set X in coefficients
-    coeffs_[0] = x.x();
-    coeffs_[1] = x.y();
-    coeffs_[2] = x.z();
-}
-
-
-void Foam::quaternionSixDOF::setU(const vector& U)
-{
-    U_.value() = U;
-    Uaverage_ = U_;
-
-    // Set U in coefficients
-    coeffs_[3] = U.x();
-    coeffs_[4] = U.y();
-    coeffs_[5] = U.z();
-}
-
-
-void Foam::quaternionSixDOF::setRotation(const HamiltonRodriguezRot& rot)
-{
-    // Cannot call update rotation: simply set up coefficients
-    // HJ, 23/Mar/2015
-
-    coeffs_[9] = rot.e0();
-    coeffs_[10] = rot.e1();
-    coeffs_[11] = rot.e2();
-    coeffs_[12] = rot.e3();
-}
-
-
-void Foam::quaternionSixDOF::setOmega(const vector& omega)
-{
-    omega_.value() = omega;
-    omegaAverage_ = omega_;
-    omegaAverageAbsolute_ = omega_;
-
-    // Set omega in coefficients
-    coeffs_[6] = omega.x();
-    coeffs_[7] = omega.y();
-    coeffs_[8] = omega.z();
-}
-
-
-void Foam::quaternionSixDOF::setReferentRotation
-(
-    const HamiltonRodriguezRot& rot
-)
-{
-    referentRotation_ = rot;
-
-    referentMotionConstraints_ = true;
-}
-
-
 void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
 {
     // First set the state in base class subobject
@@ -462,25 +386,12 @@ void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
 }
 
 
-Foam::vector Foam::quaternionSixDOF::rotVectorAverage() const
-{
-    return rotation_.rotVectorAverage();
-}
-
-
 const Foam::dimensionedVector& Foam::quaternionSixDOF::omegaAverage() const
 {
     return omegaAverage_;
 }
 
 
-const Foam::dimensionedVector&
-Foam::quaternionSixDOF::omegaAverageAbsolute() const
-{
-    return omegaAverageAbsolute_;
-}
-
-
 Foam::tensor Foam::quaternionSixDOF::toRelative() const
 {
     return rotation_.eCurrent().R();
@@ -642,14 +553,15 @@ bool Foam::quaternionSixDOF::writeData(Ostream& os) const
     os.writeKeyword("type") << tab << type() << token::END_STATEMENT << endl;
 
     // Write data
-    os.writeKeyword("Xrel") << tab << this->Xrel()
+    os.writeKeyword("Xrel") << tab << Xrel_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl;
-    os.writeKeyword("rotationVector") << tab << this->rotVector()
+    os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationVector") << tab << rotation_.rotVector()
         << token::END_STATEMENT << nl;
-    os.writeKeyword("rotationAngle") << tab << this->rotAngle()
+    os.writeKeyword("rotationAngle") << tab
+        << dimensionedScalar("rotAngle", dimless, rotation_.rotAngle())
         << token::END_STATEMENT << nl;
-    os.writeKeyword("omega") << tab << this->omega()
+    os.writeKeyword("omega") << tab << omega_
         << token::END_STATEMENT << nl << nl;
 
     os.writeKeyword("fixedSurge") << tab << fixedSurge_ <<
diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H
index 648b49921..9eed2941b 100644
--- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H
+++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H
@@ -58,7 +58,7 @@ class quaternionSixDOF
 {
     // Private data
 
-        // Initial body state variables
+        // Body state variables
 
             //- Displacement relative to spring equilibrium
             dimensionedVector Xrel_;
@@ -66,7 +66,7 @@ class quaternionSixDOF
             //- Velocity of mass centroid
             dimensionedVector U_;
 
-            //- Average velocity of mass centroid at previous time-step
+            //- Average velocity of mass centroid (evaluated at midstep)
             dimensionedVector Uaverage_;
 
             //- Finite rotation
@@ -75,10 +75,8 @@ class quaternionSixDOF
             //- Rotational velocity about mass centroid
             dimensionedVector omega_;
 
-
-        // Average variables that need to be stored
-
             //- Average rotational velocity in relative coordinate system
+            //  (evaluated at midstep)
             dimensionedVector omegaAverage_;
 
             //- Average rotational velocity in absolute coordinate system
@@ -89,7 +87,7 @@ class quaternionSixDOF
         scalarField coeffs_;
 
 
-        //- Motion constraints (given as fixed motion components)
+        // Motion constraints (given as fixed motion components)
 
             //- Fixed surge (x-translation)
             Switch fixedSurge_;
@@ -201,7 +199,7 @@ public:
                 virtual const dimensionedVector& omega() const;
 
 
-            // Displacement and rotation in the absolute coordinate system
+            // Displacement and velocity in the absolute coordinate system
 
                 //- Return position of origin in absolute coordinate system
                 virtual dimensionedVector X() const;
@@ -209,56 +207,22 @@ public:
                 //- Return velocity of origin
                 virtual const dimensionedVector& U() const;
 
-                //- Return average velocity of origin for the previous time-step
+                //- Return average velocity of origin (evaluated at midstep)
                 virtual const dimensionedVector& Uaverage() const;
 
-                //- Return finite rotation
-                virtual const finiteRotation& rotation() const;
-
-                //- Return rotational vector of body
-                virtual vector rotVector() const;
-
-                //- Return rotation angle of body
-                virtual dimensionedScalar rotAngle() const;
-
 
             // Non-access control for setting state variables
 
-                //- Set position of origin
-                virtual void setXrel(const vector& x);
-
-                //- Set velocity of origin
-                virtual void setU(const vector& u);
-
-                //- Set rotational angle in relative coordinate system
-                virtual void setRotation(const HamiltonRodriguezRot& rot);
-
-                //- Set rotational velocity in relative coordinate system
-                virtual void setOmega(const vector& omega);
-
-                //- Set referent coordinate system to apply constraints
-                virtual void setReferentRotation
-                (
-                    const HamiltonRodriguezRot& rot
-                );
-
                 //- Set ODE parameters from another ODE
                 virtual void setState(const sixDOFODE&);
 
 
             // Average motion per time-step
 
-                //- Return average rotational vector of body
-                virtual vector rotVectorAverage() const;
-
                 //- Return average rotational velocity in relative coordinate
-                //  system
+                //  system (evaluated at midstep)
                 virtual const dimensionedVector& omegaAverage() const;
 
-                //- Return average rotational velocity in absolute coordinate
-                //  system
-                virtual const dimensionedVector& omegaAverageAbsolute() const;
-
 
             // Rotations
 
@@ -325,6 +289,7 @@ public:
                 //- Update ODE after the solution, advancing by delta
                 virtual void update(const scalar delta);
 
+
         // Write controls
 
             //- WriteData member function required by regIOobject
diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C
index 02c734c53..e0c70b045 100644
--- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C
+++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C
@@ -48,7 +48,7 @@ defineRunTimeSelectionTable(sixDOFODE, dictionary);
 
 }
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
 
 void Foam::sixDOFODE::aitkensRelaxation
 (
@@ -173,7 +173,7 @@ void Foam::sixDOFODE::relaxAcceleration
 {
     if (mag(minRelFactor - maxRelFactor) < SMALL)
     {
-        // Fixed relaxation
+       // Fixed relaxation
        relaxFactorT_ = minRelFactor;
        relaxFactorR_ = minRelFactor;
     }
diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
index c9a233793..70ab14f38 100644
--- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
+++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
@@ -135,19 +135,19 @@ class sixDOFODE
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct
-        sixDOFODE(const sixDOFODE&);
+        // Copy control
 
-        //- Disallow default bitwise assignment
-        void operator=(const sixDOFODE&);
+            //- Disallow default bitwise copy construct
+            sixDOFODE(const sixDOFODE&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const sixDOFODE&);
 
 
-protected:
+        // Aitkens relaxation helper function
 
-    // Protected Member Functions
-
-        //- Update Aitkens relaxation parameters
-        void aitkensRelaxation(const scalar min, const scalar max);
+            //- Update Aitkens relaxation parameters
+            void aitkensRelaxation(const scalar min, const scalar max);
 
 
 public:
@@ -182,6 +182,7 @@ public:
         //- Return autoPtr to the selected sixDOFODE
         static autoPtr<sixDOFODE> New(const IOobject& io);
 
+
     // Destructor
 
         virtual ~sixDOFODE();
@@ -271,7 +272,7 @@ public:
                 virtual const dimensionedVector& omega() const = 0;
 
 
-            // Displacement and rotation in the absolute coordinate system
+            // Displacement and velocity in the absolute coordinate system
 
                 //- Return position of origin in absolute coordinate system
                 virtual dimensionedVector X() const = 0;
@@ -279,57 +280,22 @@ public:
                 //- Return velocity of origin
                 virtual const dimensionedVector& U() const = 0;
 
-                //- Return average velocity of origin for the previous time-step
+                //- Return average velocity of origin (evaluated at midstep)
                 virtual const dimensionedVector& Uaverage() const = 0;
 
-                //- Return finite rotation
-                virtual const finiteRotation& rotation() const = 0;
-
-                //- Return rotational vector of body
-                virtual vector rotVector() const = 0;
-
-                //- Return rotation angle of body
-                virtual dimensionedScalar rotAngle() const = 0;
-
 
             // Non-access control for setting state variables
 
-                //- Set position of origin
-                virtual void setXrel(const vector& x) = 0;
-
-                //- Set velocity of origin
-                virtual void setU(const vector& u) = 0;
-
-                //- Set rotational angle in relative coordinate system
-                virtual void setRotation(const HamiltonRodriguezRot& rot) = 0;
-
-                //- Set rotational velocity in relative coordinate system
-                virtual void setOmega(const vector& omega) = 0;
-
-                //- Set referent coordinate system to apply constraints
-                virtual void setReferentRotation
-                (
-                    const HamiltonRodriguezRot& rot
-                ) = 0;
-
                 //- Set ODE parameters from another ODE
                 virtual void setState(const sixDOFODE&);
 
 
             // Average motion per time-step
 
-                //- Return average rotational vector of body
-                virtual vector rotVectorAverage() const = 0;
-
                 //- Return average rotational velocity in relative coordinate
-                //  system
+                //  system (evaluated at midstep)
                 virtual const dimensionedVector& omegaAverage() const = 0;
 
-                //- Return average rotational velocity in absolute coordinate
-                //  system
-                virtual const dimensionedVector&
-                omegaAverageAbsolute() const = 0;
-
 
             // Rotations