diff --git a/CMakeLists.txt b/CMakeLists.txt
index 39338b4..b69d4aa 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -114,6 +114,7 @@ if(BUILD_TESTS)
   if(BUILD_TESTS_WITH_EIGEN)
     string(REPLACE "Conversions" "ConversionsEigen" TESTS_SRC_EIGEN "${TEST_SRC}")
     string(REPLACE "Model" "ModelEigen" TESTS_SRC_EIGEN "${TEST_SRC}")
+    string(REPLACE "Motion" "MotionEigen" TESTS_SRC_EIGEN "${TEST_SRC}")
     add_executable(${TEST_NAME}_eigen ${TESTS_SRC_EIGEN})
     add_test(NAME ${TEST_NAME}_eigen COMMAND "${TEST_PATH}/${TEST_NAME}_eigen")
   endif(BUILD_TESTS_WITH_EIGEN)
diff --git a/ProjectFiles.cmake b/ProjectFiles.cmake
index ce50e5d..4a0bca6 100644
--- a/ProjectFiles.cmake
+++ b/ProjectFiles.cmake
@@ -7,6 +7,7 @@ set(TEST_SRC
   "${TEST_SRC_PATH}/testAstro.cpp"
   "${TEST_SRC_PATH}/testConstants.cpp"
   "${TEST_SRC_PATH}/testOrbitalElementConversions.cpp"
+  "${TEST_SRC_PATH}/testRelativeMotion.cpp"
   "${TEST_SRC_PATH}/testStateVectorIndices.cpp"
   "${TEST_SRC_PATH}/testSolarRadiationPressureAccelerationModel.cpp"
   "${TEST_SRC_PATH}/testTwoBodyMethods.cpp"
diff --git a/include/Astro/astro.hpp b/include/Astro/astro.hpp
index e07278b..c8ea0d8 100644
--- a/include/Astro/astro.hpp
+++ b/include/Astro/astro.hpp
@@ -9,6 +9,7 @@
 
 #include "Astro/constants.hpp"
 #include "Astro/orbitalElementConversions.hpp"
+#include "Astro/relativeMotion.hpp"
 #include "Astro/twoBodyMethods.hpp"
 #include "Astro/solarRadiationPressureAccelerationModel.hpp"
 #include "Astro/stateVectorIndices.hpp"
diff --git a/include/Astro/relativeMotion.hpp b/include/Astro/relativeMotion.hpp
new file mode 100644
index 0000000..d4450a1
--- /dev/null
+++ b/include/Astro/relativeMotion.hpp
@@ -0,0 +1,147 @@
+/*
+ * Copyright (c) 2016 Kartik Kumar, Dinamica Srl (me@kartikkumar.com)
+ * Distributed under the MIT License.
+ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
+ */
+
+#ifndef ASTRO_RELATIVE_MOTION_HPP
+#define ASTRO_RELATIVE_MOTION_HPP
+
+#include <cmath>
+
+#include "Astro/stateVectorIndices.hpp"
+
+namespace astro
+{
+
+//! Propagate constant-force Clohessy-Wiltshire solution.
+/*!
+ * Propagates the constant-force Clohessy-Wiltshire solution to the Hill equations. The
+ * Clohessy-Wiltshire solution is a closed-form solution to the equations of motion describing the
+ * relative motion between a target and a chaser. The solution is predicated on the assumption that
+ * the target is on a circular orbit. The solution is expressed in terms of the Hill frame. The
+ * reader is referred to Fehse (2003) for a concise statement of the problem and an in-depth
+ * explanation of the derivation of this solution.
+ *
+ * The Clohessy-Wiltshire solution can be evaluated at a specified final time (\f$t_{2}\f$), given
+ * an initial state (\f$\bar{x}_{0}\f$) in the Hill frame, to obtain the chaser's final state. A
+ * constant thrust acceleration can be specified, i.e., a constant acceleration is added to the
+ * right-hand side of the equations of motion.
+ *
+ * It should be noted that the definition of the Hill frame conforms to the norm outlined by Fehse
+ * (2003). The user should be careful to adhere to this when executing this function.
+ *
+ * @sa computeUnperturbedHillEquations
+ * @tparam  Real                   Real type
+ * @tparam  Vector6                6-Vector type
+ * @tparam  Vector3                3-Vector type
+ * @param   initialState           Cartesian state at initial epoch in Hill frame
+ * @param   finalTime              Final epoch
+ * @param   targetMeanMotion       Mean motion of target orbit (circular)
+ * @param   thrustAcceleration     Thrust acceleration to be applied to chaser, in Hill frame
+ * @return                         Final state in Hill frame
+ */
+template< typename Real, typename Vector6, typename Vector3 >
+Vector6 propagateClohessyWiltshireSolution( const Vector6& initialState,
+                                            const Real finalTime,
+                                            const Real targetMeanMotion,
+                                            const Vector3& thrustAcceleration )
+{
+    Vector6 finalState = initialState;
+
+    finalState[ astro::xPositionIndex ]
+        = // unperturbed terms
+          ( 4.0 / targetMeanMotion * initialState[ astro::xVelocityIndex ]
+            - 6.0 * initialState[ astro::zPositionIndex ] )
+            * std::sin( targetMeanMotion * finalTime )
+          - 2.0 * initialState[ astro::zVelocityIndex ] / targetMeanMotion
+            * std::cos( targetMeanMotion * finalTime )
+          + ( 6.0 * targetMeanMotion * initialState[ astro::zPositionIndex ]
+              - 3.0 * initialState[ astro::xVelocityIndex ] ) * finalTime
+          + ( initialState[ astro::xPositionIndex ]
+              + 2.0 * initialState[ astro::zVelocityIndex ] / targetMeanMotion )
+          // constant force terms
+          + 2.0 / ( targetMeanMotion * targetMeanMotion )
+            * thrustAcceleration[ astro::zPositionIndex ]
+            * ( targetMeanMotion * finalTime - std::sin( targetMeanMotion * finalTime ) )
+          + thrustAcceleration[ astro::xPositionIndex ]
+            * ( 4.0 / ( targetMeanMotion * targetMeanMotion )
+                * ( 1.0 - std::cos( targetMeanMotion * finalTime ) )
+                - 1.5 * finalTime * finalTime );
+
+    finalState[ astro::yPositionIndex ]
+        = // unperturbed terms
+          initialState[ astro::yPositionIndex ] * std::cos( targetMeanMotion * finalTime )
+          + initialState[ astro::yVelocityIndex ] / targetMeanMotion
+            * std::sin( targetMeanMotion * finalTime )
+          // constant force terms
+          + thrustAcceleration[ astro::yPositionIndex ] / ( targetMeanMotion * targetMeanMotion )
+            * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+
+    finalState[ astro::zPositionIndex ]
+        = // unperturbed terms
+          ( 2.0 * initialState[ astro::xVelocityIndex ] / targetMeanMotion
+            - 3.0 * initialState[ astro::zPositionIndex ] )
+            * std::cos( targetMeanMotion * finalTime )
+          + initialState[ astro::zVelocityIndex ] / targetMeanMotion
+          * std::sin( targetMeanMotion * finalTime )
+          + ( 4.0 * initialState[ astro::zPositionIndex ]
+            - 2.0 * initialState[ astro::xVelocityIndex ] / targetMeanMotion )
+          // constant force terms
+          + 2.0 / ( targetMeanMotion * targetMeanMotion )
+            * thrustAcceleration[ astro::xPositionIndex ]
+            * ( std::sin( targetMeanMotion * finalTime ) - targetMeanMotion * finalTime )
+          + thrustAcceleration[ astro::zPositionIndex ] / ( targetMeanMotion * targetMeanMotion )
+            * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+
+    finalState[ astro::xVelocityIndex ]
+        = // unperturbed terms
+          targetMeanMotion * ( 4.0 / targetMeanMotion * initialState[ astro::xVelocityIndex ]
+            - 6.0 * initialState[ astro::zPositionIndex ] )
+            * std::cos( targetMeanMotion * finalTime )
+          + 2.0 * initialState[ astro::zVelocityIndex ]
+            * std::sin( targetMeanMotion * finalTime )
+          + ( 6.0 * targetMeanMotion * initialState[ astro::zPositionIndex]
+              - 3.0 * initialState[ astro::xVelocityIndex ] )
+          // constant force terms
+          + thrustAcceleration[ astro::zPositionIndex ]
+            * 2.0 / ( targetMeanMotion * targetMeanMotion )
+            * ( targetMeanMotion - targetMeanMotion * std::cos( targetMeanMotion * finalTime ) )
+          + thrustAcceleration[ astro::xPositionIndex ]
+            * ( 4.0 / targetMeanMotion * std::sin( targetMeanMotion * finalTime )
+                - 3.0 * finalTime );
+
+    finalState[ astro::yVelocityIndex ]
+        = // unperturbed terms
+          -targetMeanMotion * initialState[ astro::yPositionIndex ]
+            * std::sin( targetMeanMotion * finalTime )
+          + initialState[ astro::yVelocityIndex ] * std::cos( targetMeanMotion * finalTime )
+          // constant force terms
+          + thrustAcceleration[ astro::yPositionIndex ] / targetMeanMotion
+            * std::sin( targetMeanMotion * finalTime );
+
+    finalState[ astro::zVelocityIndex ]
+        = // unperturbed terms
+          -targetMeanMotion * ( 2.0 * initialState[ astro::xVelocityIndex ] / targetMeanMotion
+            - 3.0 * initialState[ astro::zPositionIndex ] )
+            * std::sin( targetMeanMotion * finalTime )
+          + initialState[ astro::zVelocityIndex ] * std::cos( targetMeanMotion * finalTime )
+          // constant force terms
+          + thrustAcceleration[ astro::xPositionIndex ]
+            * 2.0 / ( targetMeanMotion * targetMeanMotion )
+            * ( targetMeanMotion * std::cos( targetMeanMotion * finalTime ) - targetMeanMotion )
+          + thrustAcceleration[ astro::zPositionIndex ] / targetMeanMotion
+            * std::sin( targetMeanMotion * finalTime );
+
+    return finalState;
+}
+
+} // namespace astro
+
+#endif // ASTRO_RELATIVE_MOTION_HPP
+
+/*!
+ * References
+ *  Fehse, W. (2003) Automated Rendezvous and Docking of Spacecraft, Cambridge Aerospace Series 16,
+ *    Cambridge University Press, Cambridge, United Kingdom.
+ */
diff --git a/test/testRelativeMotion.cpp b/test/testRelativeMotion.cpp
new file mode 100644
index 0000000..9172c9e
--- /dev/null
+++ b/test/testRelativeMotion.cpp
@@ -0,0 +1,367 @@
+/*
+ * Copyright (c) 2014-2016 Kartik Kumar, Dinamica Srl (me@kartikkumar.com)
+ * Distributed under the MIT License.
+ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
+ */
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+#include <catch.hpp>
+
+#include "Astro/relativeMotion.hpp"
+
+namespace astro
+{
+namespace tests
+{
+
+typedef double Real;
+typedef std::vector< Real > Vector;
+static const Real tolerance = 100.0 * std::numeric_limits< Real >::epsilon( );
+
+TEST_CASE( "Test Clohessy-Wiltshire solution", "[clohessy-wiltshire,relative-motion]" )
+{
+    // Scenarios taken from chapter 3, Fehse (2003).
+
+    const Real earthGravitationalParameter = 398600.14e9;       // m^3 s^-2
+    const Real targetSemiMajorAxis = 7200.0e3;                  // m
+    const Real targetMeanMotion
+        = std::sqrt( earthGravitationalParameter
+            / ( targetSemiMajorAxis * targetSemiMajorAxis * targetSemiMajorAxis ) );
+    const Real finalTime = 100.0;                               // s
+
+    SECTION( "Test free motion" )
+    {
+        Vector thrustAcceleration( 3 );
+        thrustAcceleration[ 0 ] = 0.0;
+        thrustAcceleration[ 1 ] = 0.0;
+        thrustAcceleration[ 2 ] = 0.0;
+
+        SECTION( "Test different altitude" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 0.0;
+            initialState[ 1 ] = 0.0;
+            initialState[ 2 ] = 10.0;
+            initialState[ 3 ] = 1.5 * targetMeanMotion * 10.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 1.5 * targetMeanMotion * 10.0 * finalTime;
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ] = 10.0;
+            expectedFinalState[ 3 ] = 1.5 * targetMeanMotion * 10.0;
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test different altitude with target's velocity" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 0.0;
+            initialState[ 1 ] = 0.0;
+            initialState[ 2 ] = 10.0;
+            initialState[ 3 ] = 0.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 6.0 * 10.0 * ( targetMeanMotion * finalTime
+                - std::sin( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ]
+                = 10.0 * ( 4.0 - 3.0 * std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 3 ] = 6.0 * 10.0 * targetMeanMotion
+                * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ]
+                = 3.0 * 10.0 * targetMeanMotion * std::sin( targetMeanMotion * finalTime );
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test out-of-plane free drift" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 0.0;
+            initialState[ 1 ] = 20.0;
+            initialState[ 2 ] = 0.0;
+            initialState[ 3 ] = 0.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 0.0;
+            expectedFinalState[ 1 ] = 20.0 * std::cos( targetMeanMotion * finalTime );
+            expectedFinalState[ 2 ] = 0.0;
+            expectedFinalState[ 3 ] = 0.0;
+            expectedFinalState[ 4 ]
+                = -20.0 * targetMeanMotion * std::sin( targetMeanMotion * finalTime );
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test along-track free drift" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 5.0;
+            initialState[ 1 ] = 0.0;
+            initialState[ 2 ] = 0.0;
+            initialState[ 3 ] = 0.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 5.0;
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ] = 0.0;
+            expectedFinalState[ 3 ] = 0.0;
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+    }
+
+    SECTION( "Test constant-force motion applied at the target" )
+    {
+        Vector initialState( 6 );
+        initialState[ 0 ] = 0.0;
+        initialState[ 1 ] = 0.0;
+        initialState[ 2 ] = 0.0;
+        initialState[ 3 ] = 0.0;
+        initialState[ 4 ] = 0.0;
+        initialState[ 5 ] = 0.0;
+
+        SECTION( "Test along-track force" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 2.161;
+            thrustAcceleration[ 1 ] = 0.0;
+            thrustAcceleration[ 2 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ]
+                = thrustAcceleration[ 0 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 4.0 * ( 1.0 - std::cos( targetMeanMotion * finalTime ) )
+                      - 1.5 * targetMeanMotion * targetMeanMotion * finalTime * finalTime );
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ]
+                = 2.0 * thrustAcceleration[ 0 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( std::sin( targetMeanMotion * finalTime ) - targetMeanMotion * finalTime );
+            expectedFinalState[ 3 ]
+                = thrustAcceleration[ 0 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 4.0 * targetMeanMotion * std::sin( targetMeanMotion * finalTime )
+                      - 3.0 * targetMeanMotion * targetMeanMotion * finalTime );
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ]
+                = 2.0 * thrustAcceleration[ 0 ] / targetMeanMotion
+                  * ( std::cos( targetMeanMotion * finalTime ) - 1.0 );
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test cross-track force" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 0.0;
+            thrustAcceleration[ 1 ] = -4.915;
+            thrustAcceleration[ 2 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 0.0;
+            expectedFinalState[ 1 ]
+                = thrustAcceleration[ 1 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 2 ] = 0.0;
+            expectedFinalState[ 3 ] = 0.0;
+            expectedFinalState[ 4 ]
+                = thrustAcceleration[ 1 ] / targetMeanMotion
+                  * std::sin( targetMeanMotion * finalTime );
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test radial force" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 0.0;
+            thrustAcceleration[ 1 ] = 0.0;
+            thrustAcceleration[ 2 ] = 5.842;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ]
+                = 2.0 * thrustAcceleration[ 2 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( targetMeanMotion * finalTime - std::sin( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ]
+                = thrustAcceleration [ 2 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 3 ]
+                = 2.0 * thrustAcceleration [ 2 ] / targetMeanMotion
+                  * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ]
+                = thrustAcceleration[ 2 ] / targetMeanMotion
+                  * std::sin( targetMeanMotion * finalTime );
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+    }
+
+    SECTION( "Test internal snapshot" )
+    {
+        // The tests presented here are only to confirm internal consistency by testing against a
+        // snapshot of the results obtained at the time of writing the test.
+
+        SECTION( "Test arbitrary free motion" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 0.0;
+            thrustAcceleration[ 1 ] = 0.0;
+            thrustAcceleration[ 2 ] = 0.0;
+
+            Vector initialState( 6 );
+            initialState[ 0 ] = 15.613;
+            initialState[ 1 ] = -1.6136;
+            initialState[ 2 ] = 43.123;
+            initialState[ 3 ] = -1.35;
+            initialState[ 4 ] = 0.612;
+            initialState[ 5 ] = -5.699;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = -177.220096793708;
+            expectedFinalState[ 1 ] = 59.4861383364061;
+            expectedFinalState[ 2 ] = -511.134488585295;
+            expectedFinalState[ 3 ] = -2.49554339092689;
+            expectedFinalState[ 4 ] = 0.608907076143407;
+            expectedFinalState[ 5 ] = -5.3762829430049;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test arbitrary constant-force motion" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 3.415;
+            thrustAcceleration[ 1 ] = -1.556;
+            thrustAcceleration[ 2 ] = 8.821;
+
+            Vector initialState( 6 );
+            initialState[ 0 ] = 15.613;
+            initialState[ 1 ] = -1.6136;
+            initialState[ 2 ] = 43.123;
+            initialState[ 3 ] = -1.35;
+            initialState[ 4 ] = 0.612;
+            initialState[ 5 ] = -5.699;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 19873.9479718838;
+            expectedFinalState[ 1 ] = -7713.59262479076;
+            expectedFinalState[ 2 ] = 42378.8990414098;
+            expectedFinalState[ 3 ] = 427.649888486457;
+            expectedFinalState[ 4 ] = -154.714292723343;
+            expectedFinalState[ 5 ] = 839.895191983715;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+    }
+}
+
+} // namespace tests
+} // namespace astro
+
+/*!
+ * References
+ *  Fehse, W. (2003) Automated Rendezvous and Docking of Spacecraft, Cambridge Aerospace Series 16,
+ *    Cambridge University Press, Cambridge, United Kingdom.
+ */
diff --git a/test/testRelativeMotionEigen.cpp b/test/testRelativeMotionEigen.cpp
new file mode 100644
index 0000000..6cfd1d0
--- /dev/null
+++ b/test/testRelativeMotionEigen.cpp
@@ -0,0 +1,368 @@
+/*
+ * Copyright (c) 2014-2016 Kartik Kumar, Dinamica Srl (me@kartikkumar.com)
+ * Distributed under the MIT License.
+ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
+ */
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+#include <catch.hpp>
+
+#include <Eigen/Core>
+
+#include "Astro/relativeMotion.hpp"
+
+namespace astro
+{
+namespace tests
+{
+
+typedef double Real;
+typedef Eigen::Matrix< Real, Eigen::Dynamic, 1 > Vector;
+static const Real tolerance = 100.0 * std::numeric_limits< Real >::epsilon( );
+
+TEST_CASE( "Test Clohessy-Wiltshire solution", "[clohessy-wiltshire,relative-motion]" )
+{
+    // Scenarios taken from chapter 3, Fehse (2003).
+
+    const Real earthGravitationalParameter = 398600.14e9;       // m^3 s^-2
+    const Real targetSemiMajorAxis = 7200.0e3;                  // m
+    const Real targetMeanMotion
+        = std::sqrt( earthGravitationalParameter
+            / ( targetSemiMajorAxis * targetSemiMajorAxis * targetSemiMajorAxis ) );
+    const Real finalTime = 100.0;                               // s
+
+    SECTION( "Test free motion" )
+    {
+        Vector thrustAcceleration( 3 );
+        thrustAcceleration[ 0 ] = 0.0;
+        thrustAcceleration[ 1 ] = 0.0;
+        thrustAcceleration[ 2 ] = 0.0;
+
+        SECTION( "Test different altitude" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 0.0;
+            initialState[ 1 ] = 0.0;
+            initialState[ 2 ] = 10.0;
+            initialState[ 3 ] = 1.5 * targetMeanMotion * 10.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 1.5 * targetMeanMotion * 10.0 * finalTime;
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ] = 10.0;
+            expectedFinalState[ 3 ] = 1.5 * targetMeanMotion * 10.0;
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test different altitude with target's velocity" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 0.0;
+            initialState[ 1 ] = 0.0;
+            initialState[ 2 ] = 10.0;
+            initialState[ 3 ] = 0.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 6.0 * 10.0 * ( targetMeanMotion * finalTime
+                - std::sin( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ]
+                = 10.0 * ( 4.0 - 3.0 * std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 3 ] = 6.0 * 10.0 * targetMeanMotion
+                * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ]
+                = 3.0 * 10.0 * targetMeanMotion * std::sin( targetMeanMotion * finalTime );
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test out-of-plane free drift" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 0.0;
+            initialState[ 1 ] = 20.0;
+            initialState[ 2 ] = 0.0;
+            initialState[ 3 ] = 0.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 0.0;
+            expectedFinalState[ 1 ] = 20.0 * std::cos( targetMeanMotion * finalTime );
+            expectedFinalState[ 2 ] = 0.0;
+            expectedFinalState[ 3 ] = 0.0;
+            expectedFinalState[ 4 ]
+                = -20.0 * targetMeanMotion * std::sin( targetMeanMotion * finalTime );
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test along-track free drift" )
+        {
+            Vector initialState( 6 );
+            initialState[ 0 ] = 5.0;
+            initialState[ 1 ] = 0.0;
+            initialState[ 2 ] = 0.0;
+            initialState[ 3 ] = 0.0;
+            initialState[ 4 ] = 0.0;
+            initialState[ 5 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 5.0;
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ] = 0.0;
+            expectedFinalState[ 3 ] = 0.0;
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+    }
+
+    SECTION( "Test constant-force motion applied at the target" )
+    {
+        Vector initialState( 6 );
+        initialState[ 0 ] = 0.0;
+        initialState[ 1 ] = 0.0;
+        initialState[ 2 ] = 0.0;
+        initialState[ 3 ] = 0.0;
+        initialState[ 4 ] = 0.0;
+        initialState[ 5 ] = 0.0;
+
+        SECTION( "Test along-track force" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 2.161;
+            thrustAcceleration[ 1 ] = 0.0;
+            thrustAcceleration[ 2 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ]
+                = thrustAcceleration[ 0 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 4.0 * ( 1.0 - std::cos( targetMeanMotion * finalTime ) )
+                      - 1.5 * targetMeanMotion * targetMeanMotion * finalTime * finalTime );
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ]
+                = 2.0 * thrustAcceleration[ 0 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( std::sin( targetMeanMotion * finalTime ) - targetMeanMotion * finalTime );
+            expectedFinalState[ 3 ]
+                = thrustAcceleration[ 0 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 4.0 * targetMeanMotion * std::sin( targetMeanMotion * finalTime )
+                      - 3.0 * targetMeanMotion * targetMeanMotion * finalTime );
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ]
+                = 2.0 * thrustAcceleration[ 0 ] / targetMeanMotion
+                  * ( std::cos( targetMeanMotion * finalTime ) - 1.0 );
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test cross-track force" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 0.0;
+            thrustAcceleration[ 1 ] = -4.915;
+            thrustAcceleration[ 2 ] = 0.0;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 0.0;
+            expectedFinalState[ 1 ]
+                = thrustAcceleration[ 1 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 2 ] = 0.0;
+            expectedFinalState[ 3 ] = 0.0;
+            expectedFinalState[ 4 ]
+                = thrustAcceleration[ 1 ] / targetMeanMotion
+                  * std::sin( targetMeanMotion * finalTime );
+            expectedFinalState[ 5 ] = 0.0;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test radial force" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 0.0;
+            thrustAcceleration[ 1 ] = 0.0;
+            thrustAcceleration[ 2 ] = 5.842;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ]
+                = 2.0 * thrustAcceleration[ 2 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( targetMeanMotion * finalTime - std::sin( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 1 ] = 0.0;
+            expectedFinalState[ 2 ]
+                = thrustAcceleration [ 2 ] / ( targetMeanMotion * targetMeanMotion )
+                  * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 3 ]
+                = 2.0 * thrustAcceleration [ 2 ] / targetMeanMotion
+                  * ( 1.0 - std::cos( targetMeanMotion * finalTime ) );
+            expectedFinalState[ 4 ] = 0.0;
+            expectedFinalState[ 5 ]
+                = thrustAcceleration[ 2 ] / targetMeanMotion
+                  * std::sin( targetMeanMotion * finalTime );
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+    }
+
+    SECTION( "Test internal snapshot" )
+    {
+        // The tests presented here are only to confirm internal consistency by testing against a
+        // snapshot of the results obtained at the time of writing the test.
+
+        SECTION( "Test arbitrary free motion" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 0.0;
+            thrustAcceleration[ 1 ] = 0.0;
+            thrustAcceleration[ 2 ] = 0.0;
+
+            Vector initialState( 6 );
+            initialState[ 0 ] = 15.613;
+            initialState[ 1 ] = -1.6136;
+            initialState[ 2 ] = 43.123;
+            initialState[ 3 ] = -1.35;
+            initialState[ 4 ] = 0.612;
+            initialState[ 5 ] = -5.699;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = -177.220096793708;
+            expectedFinalState[ 1 ] = 59.4861383364061;
+            expectedFinalState[ 2 ] = -511.134488585295;
+            expectedFinalState[ 3 ] = -2.49554339092689;
+            expectedFinalState[ 4 ] = 0.608907076143407;
+            expectedFinalState[ 5 ] = -5.3762829430049;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+
+        SECTION( "Test arbitrary constant-force motion" )
+        {
+            Vector thrustAcceleration( 3 );
+            thrustAcceleration[ 0 ] = 3.415;
+            thrustAcceleration[ 1 ] = -1.556;
+            thrustAcceleration[ 2 ] = 8.821;
+
+            Vector initialState( 6 );
+            initialState[ 0 ] = 15.613;
+            initialState[ 1 ] = -1.6136;
+            initialState[ 2 ] = 43.123;
+            initialState[ 3 ] = -1.35;
+            initialState[ 4 ] = 0.612;
+            initialState[ 5 ] = -5.699;
+
+            Vector expectedFinalState( 6 );
+            expectedFinalState[ 0 ] = 19873.9479718838;
+            expectedFinalState[ 1 ] = -7713.59262479076;
+            expectedFinalState[ 2 ] = 42378.8990414098;
+            expectedFinalState[ 3 ] = 427.649888486457;
+            expectedFinalState[ 4 ] = -154.714292723343;
+            expectedFinalState[ 5 ] = 839.895191983715;
+
+            Vector computedFinalState = propagateClohessyWiltshireSolution( initialState,
+                                                                            finalTime,
+                                                                            targetMeanMotion,
+                                                                            thrustAcceleration );
+
+            for ( unsigned int i = 0; i < 6; ++i )
+            {
+                REQUIRE( computedFinalState[ i ]
+                            == Approx( expectedFinalState[ i ] ).epsilon( tolerance ) );
+            }
+        }
+    }
+}
+
+} // namespace tests
+} // namespace astro
+
+/*!
+ * References
+ *  Fehse, W. (2003) Automated Rendezvous and Docking of Spacecraft, Cambridge Aerospace Series 16,
+ *    Cambridge University Press, Cambridge, United Kingdom.
+ */
diff --git a/test/testSolarRadiationPressureAccelerationModelEigen.cpp b/test/testSolarRadiationPressureAccelerationModelEigen.cpp
index aaafc73..f71edde 100644
--- a/test/testSolarRadiationPressureAccelerationModelEigen.cpp
+++ b/test/testSolarRadiationPressureAccelerationModelEigen.cpp
@@ -8,7 +8,6 @@
  */
 
 #include <cmath>
-#include <vector>
 
 #include <catch.hpp>