diff --git a/Makefile b/Makefile
index 42f1b4c..d134d8f 100644
--- a/Makefile
+++ b/Makefile
@@ -15,6 +15,8 @@ COMMON_OBJECTS := $(patsubst $(SRCDIR)/%.cpp, $(OBJDIR)/%.o, $(COMMON_SOURCES))
 
 TEST_SOURCES := $(TESTDIR)/orbit_test.cpp
 TEST_OBJECTS := $(patsubst $(TESTDIR)/%.cpp, $(OBJDIR)/%.o, $(TEST_SOURCES))
+TEST_MAIN_SOURCE := $(TESTDIR)/tests_main.cpp
+TEST_MAIN_OBJ := $(OBJDIR)/tests_main.o
 TEST_BIN := $(BINDIR)/test_orbit
 
 SOURCES := $(SRCDIR)/main.cpp
@@ -31,14 +33,19 @@ $(BIN): $(COMMON_OBJECTS)
 	$(CXX) $(CXXFLAGS) -c -MMD $(SOURCES) -o $(OBJECTS)
 	$(CXX) -o $(BIN) $(LDFLAGS) $(COMMON_OBJECTS) $(OBJECTS) $(LIB)
 
-tests: $(COMMON_OBJECTS)
-	$(CXX) $(CXXFLAGS) -c -MMD $(TEST_SOURCES) -o $(TEST_OBJECTS)
-	$(CXX) -o $(TEST_BIN) $(LDFLAGS) $(COMMON_OBJECTS) $(TEST_OBJECTS) $(LIB)
+tests: $(TEST_OBJECTS) $(TEST_MAIN_OBJ)
+	$(CXX) -o $(TEST_BIN) $(LDFLAGS) $(TEST_OBJECTS) $(TEST_MAIN_OBJ) $(LIB)
 .PHONY: tests
 
 $(COMMON_OBJECTS): $(COMMON_SOURCES)
 	$(CXX) $(CXXFLAGS) -c -MMD $(SRCDIR)/orbits.cpp -o $@
 
+$(TEST_OBJECTS): $(TEST_SOURCES)
+	$(CXX) $(CXXFLAGS) -c -MMD $(TESTDIR)/orbit_test.cpp -o $(OBJDIR)/orbit_test.o
+
+$(TEST_MAIN_OBJ): $(TEST_MAIN_SOURCE)
+	$(CXX) $(CXXFLAGS) -c -MMD $(TESTDIR)/tests_main.cpp -o $(OBJDIR)/tests_main.o
+
 clean:
 	rm -rf $(OBJDIR)/* $(BIN) bin/test_orbit
 .PHONY: clean
diff --git a/README b/README
index 0e51db7..69ff99d 100644
--- a/README
+++ b/README
@@ -3,6 +3,8 @@ Work in progress orbital shipping game
 
 Dependencies:
 	libTangerine: https://gitlab.com/cinnaboot/tangerine
+	Catch2 (for tests):
+		https://github.com/catchorg/Catch2/blob/devel/docs/Readme.md
 
 TODO:
 	add git submodule for libTangerine
diff --git a/src/orbits.cpp b/src/orbits.cpp
index 947f3a5..43ca544 100644
--- a/src/orbits.cpp
+++ b/src/orbits.cpp
@@ -4,63 +4,6 @@
 #include "orbits.h"
 
 
-double
-getPropagatedPosition(orbital_elements orbit,
-	double initial_anom,
-	unsigned int time_step)
-{
-	ellipse_parameters ep = orbit.ep;
-
-	// TODO: break out steps with descriptive function names
-	// initial eccentric anomaly
-	double E1 = 2 * atan(sqrt((1 - ep.e) / (1 + ep.e)) * tan(initial_anom / 2));
-
-	// initial mean anomaly
-	double M1 = E1 - ep.e * sin(E1);
-
-	// mean motion
-	double n = sqrt(orbit.mu / pow(ep.a, 3));
-
-	// propagated mean anomaly
-	double M2 = M1 + n * (time_step);
-
-	// guess starting value for E2
-	double E2_1 = M2 + ep.e * sin(M2) + ((pow(ep.e, 2) / 2) * sin(2 * M2));
-
-	// test if guess is a solution to kepler's equation
-	const double ACCEPTABLE_ERROR = 0.00000001;
-	double E2_test = E2_1;
-
-	for (uint i = 0; i < 10; i++) {
-		double err = E2_test - ep.e * sin(E2_test) - M2;
-
-		LOG(Debug) << "-------------------------\n";
-		LOG(Debug) << "index: " << i << "\n";
-		LOG(Debug) << "E2_test: " << E2_test << "\n";
-		LOG(Debug) << "err: " << err << "\n";
-		LOG(Debug) << "err: " << err << "\n";
-
-		if (fabs(err) < ACCEPTABLE_ERROR)
-			break;
-
-		// compute derivative of the error function
-		double derr = 1 - ep.e * cos(E2_test);
-		LOG(Debug) << "derr: " << derr << "\n";
-
-		// use Newton's method to compute next trial value of E2
-		E2_test = E2_test - (err / derr);
-	}
-
-	return E2_test;
-}
-
-double
-testTrialEccAnomaly(double trial_value, double e, double M2)
-{
-	double err = trial_value - e * sin(trial_value) - M2;
-	return err;
-}
-
 ellipse_parameters
 constructEllipse(double a, double b)
 {
@@ -95,3 +38,93 @@ constructEllipse3D(ellipse_parameters ep, uint vert_count)
 	return e3d;
 }
 
+double
+getEccAnomFromTrueAnom(double ecc, double true_anom)
+{
+	return 2 * atan(sqrt((1 - ecc) / (1 + ecc)) * tan(true_anom / 2));
+}
+
+double
+getTrueAnomFromEccAnom(double ecc, double ecc_anom)
+{
+	return 2 * atan(sqrt((1 + ecc) / (1 - ecc)) * tan(ecc_anom / 2));
+}
+
+double
+getMeanAnomFromEccAnom(double ecc_anom, double ecc)
+{
+	return ecc_anom - ecc * sin(ecc_anom);
+}
+
+double
+getMeanMotion(double mu, double a)
+{
+	return sqrt(mu / pow(a, 3));
+}
+
+double
+getPropagatedMeanAnom(double mean_anom, double mean_motion, double time_step)
+{
+	return mean_anom + mean_motion * (time_step);
+}
+
+double
+getInitialTrialValue(double mean_anom, double ecc)
+{
+	return mean_anom + ecc * sin(mean_anom)
+		+ ((pow(ecc, 2) / 2) * sin(2 * mean_anom));
+}
+
+double
+getTrialError(double ecc, double test_anom, double mean_anom)
+{
+	return test_anom - ecc * sin(test_anom) - mean_anom;
+}
+
+double
+getNextTrialValue(double err, double ecc, double test_anom, double mean_anom)
+{
+	// compute derivative of the error function
+	double derr = 1 - ecc * cos(test_anom);
+
+	// use Newton's method to compute next trial value of E2
+	return test_anom - (err / derr);
+}
+
+double
+getPropagatedEccAnomaly(orbital_elements orbit,
+	double initial_anom,
+	unsigned int time_step)
+{
+	double e = orbit.ep.e;
+	double E1 = getEccAnomFromTrueAnom(e, initial_anom);
+	double M1 = getMeanAnomFromEccAnom(E1, e);
+	double n = getMeanMotion(orbit.mu, orbit.ep.a);
+	double M2 = getPropagatedMeanAnom(M1, n, time_step);
+	double E2_1 = getInitialTrialValue(M2, e);
+
+	// test if guess is a solution to kepler's equation
+	const double ACCEPTABLE_ERROR = 0.00000001;
+	double E2_test = E2_1;
+
+	for (uint i = 0; i < 10; i++) {
+		double err = getTrialError(e, E2_test, M2);
+
+		if (fabs(err) < ACCEPTABLE_ERROR)
+			break;
+
+		E2_test = getNextTrialValue(err, e, E2_test, M2);
+	}
+
+	return E2_test;
+}
+
+double
+getPropagatedTrueAnomaly(orbital_elements orbit,
+	double initial_anom,
+	unsigned int time_step)
+{
+	double ecc_anom = getPropagatedEccAnomaly(orbit, initial_anom, time_step);
+	return getTrueAnomFromEccAnom(orbit.ep.e, ecc_anom);
+}
+
diff --git a/src/orbits.h b/src/orbits.h
index 067aa95..70b9745 100644
--- a/src/orbits.h
+++ b/src/orbits.h
@@ -78,7 +78,7 @@ constructEllipse3D(ellipse_parameters ep, uint vert_count);
  *			tan(ϴ2/2) = sqrt((1+e) / (1-e)) * tan(E2/2)
  */
 double
-getPropagatedPosition(orbital_elements orbit,
+getPropagatedTrueAnomaly(orbital_elements orbit,
 	double initial_anom,
 	unsigned int time_step);
 
diff --git a/tests/orbit_test.cpp b/tests/orbit_test.cpp
index 7766212..ccfc4b6 100644
--- a/tests/orbit_test.cpp
+++ b/tests/orbit_test.cpp
@@ -1,11 +1,18 @@
 
+#include <catch2/catch.hpp>
+
 #include "dumbLog.h"
 
-#include "../src/orbits.h"
+#include "../src/orbits.cpp"
+
+// NOTE: using WithinAbs instead of Approx because:
+// 	https://github.com/catchorg/Catch2/issues/1962
+// 	https://github.com/catchorg/Catch2/issues/1507
+using Catch::Matchers::WithinAbs;
 
 
-void
-testOrbitPropagation()
+// NOTE: example 4.6 in "Space Flight Dynamics" by Craig A. Kluever
+TEST_CASE("orbit propagation", "[orbits]")
 {
 	orbital_elements orbit = {};
 	ellipse_parameters ep = {};
@@ -19,16 +26,27 @@ testOrbitPropagation()
 
 	double initial_anom = 260 * M_PI / 180; // radians
 	unsigned int time_step = 60 * 50; // seconds
-	double next_pos = getPropagatedPosition(orbit, initial_anom, time_step);
 
-	LOG(Debug) << "E2_1: " << next_pos << "\n";
+	double E1 = getEccAnomFromTrueAnom(ep.e, initial_anom);
+	REQUIRE_THAT(E1, WithinAbs(-0.8615, 1e-4));
 
-	return;
-}
+	double M1 = getMeanAnomFromEccAnom(E1, ep.e);
+	REQUIRE_THAT(M1, WithinAbs(-0.2992, 1e-4));
 
-int
-main()
-{
-	testOrbitPropagation();
-	return 0;
+	double n = getMeanMotion(orbit.mu, orbit.ep.a);
+	REQUIRE_THAT(n, WithinAbs(0.00014582, 1e-8));
+
+	double M2 = getPropagatedMeanAnom(M1, n, time_step);
+	REQUIRE_THAT(M2, WithinAbs(0.1383, 1e-5));
+
+	// TODO: could also test for the other trial values listed in table 4.1
+	double E2_1 = getInitialTrialValue(M2, ep.e);
+	REQUIRE_THAT(E2_1, WithinAbs(0.315452, 1e-5));
+
+	double ecc_anom = getPropagatedEccAnomaly(orbit, initial_anom, time_step);
+	REQUIRE_THAT(ecc_anom, WithinAbs(0.481518, 1e-5));
+
+	double true_anom = getPropagatedTrueAnomaly(orbit, initial_anom, time_step);
+	REQUIRE_THAT(true_anom, WithinAbs(1.1339, 1e-4));
 }
+
diff --git a/tests/tests_main.cpp b/tests/tests_main.cpp
new file mode 100644
index 0000000..995517b
--- /dev/null
+++ b/tests/tests_main.cpp
@@ -0,0 +1,4 @@
+
+#define CATCH_CONFIG_MAIN
+#include <catch2/catch.hpp>
+