diff --git a/Makefile b/Makefile
index 16995c5..42f1b4c 100644
--- a/Makefile
+++ b/Makefile
@@ -1,26 +1,44 @@
 
-OBJDIR = build
-SRCDIR = src
-LIBDIR = ext/tangerine
-LIB = $(LIBDIR)/build/libTangerine.a
-CXX = g++
-CXXFLAGS = -std=c++11 -g -ggdb3 -Wall -I$(LIBDIR)/include/
-LDFLAGS = -lSDL2 -lGLEW -lGL -lassimp
-BIN = bin/orbital_shipping
+OBJDIR := build
+SRCDIR := src
+LIBDIR := ext/tangerine
+TESTDIR := tests
+BINDIR := bin
+LIB := $(LIBDIR)/build/libTangerine.a
 
-SOURCES = $(wildcard $(SRCDIR)/*.cpp)
-OBJECTS = $(patsubst $(SRCDIR)/%.cpp, $(OBJDIR)/%.o, $(SOURCES))
+CXX := g++
+CXXFLAGS := -std=c++11 -g -ggdb3 -Wall -I$(LIBDIR)/include/
+LDFLAGS := -lSDL2 -lGLEW -lGL -lassimp
 
+COMMON_SOURCES := $(SRCDIR)/orbits.cpp
+COMMON_OBJECTS := $(patsubst $(SRCDIR)/%.cpp, $(OBJDIR)/%.o, $(COMMON_SOURCES))
 
-all: $(OBJECTS)
-	$(CXX) -o $(BIN) $(LDFLAGS) $(OBJDIR)/*.o $(LIB)
+TEST_SOURCES := $(TESTDIR)/orbit_test.cpp
+TEST_OBJECTS := $(patsubst $(TESTDIR)/%.cpp, $(OBJDIR)/%.o, $(TEST_SOURCES))
+TEST_BIN := $(BINDIR)/test_orbit
+
+SOURCES := $(SRCDIR)/main.cpp
+OBJECTS := $(patsubst $(SRCDIR)/%.cpp, $(OBJDIR)/%.o, $(SOURCES))
+BIN := $(BINDIR)/orbital_shipping
+
+
+all: $(BIN) tests
 .PHONY: all
 
 -include $(OBJDIR)/*.d
 
-$(OBJECTS): $(OBJDIR)/%.o : $(SRCDIR)/%.cpp
-	$(CXX) $(CXXFLAGS) -c -MMD $< -o $@
+$(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)
+.PHONY: tests
+
+$(COMMON_OBJECTS): $(COMMON_SOURCES)
+	$(CXX) $(CXXFLAGS) -c -MMD $(SRCDIR)/orbits.cpp -o $@
 
 clean:
-	rm -rf $(OBJDIR)/* $(BIN)
+	rm -rf $(OBJDIR)/* $(BIN) bin/test_orbit
 .PHONY: clean
diff --git a/src/main.cpp b/src/main.cpp
index 33a3c31..2bdc81e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -10,7 +10,6 @@
 #include "renderer.h"
 #include "util.h"
 
-#define ORBIT_IMPLEMENTATION
 #include "orbits.h"
 
 
diff --git a/src/orbits.cpp b/src/orbits.cpp
new file mode 100644
index 0000000..947f3a5
--- /dev/null
+++ b/src/orbits.cpp
@@ -0,0 +1,97 @@
+
+#include "dumbLog.h"
+
+#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)
+{
+	assert(a > 0 && b > 0 && a >= b);
+	ellipse_parameters ep = { a, b };
+	ep.c = sqrt(a * a - b * b);
+	ep.e = ep.c / ep.a;
+	ep.f1.x = -1 * ep.c;
+	ep.f2.x = ep.c;
+	return ep;
+}
+
+ellipse_3d
+constructEllipse3D(ellipse_parameters ep, uint vert_count)
+{
+	assert(ep.a > 0 && ep.b > 0 &&
+		ep.a >= ep.b &&
+		vert_count > 0);
+
+	ellipse_3d e3d = { ep, nullptr, vert_count};
+	e3d.vertices = UTIL_ALLOC(vert_count, glm::vec3);
+	double angle = 2 * M_PI / vert_count;
+
+	for (uint i = 0; i < vert_count; i++) {
+		double a = angle * i;
+		// NOTE: solving for distance in polar coordinates
+		double r = ep.b / sqrt(1 - (ep.e * pow(cos(a), 2)));
+		// NOTE: converting from polar to rectangular
+		e3d.vertices[i] = glm::vec3(cos(a) *  r, sin(a) * r, 0);
+	}
+
+	return e3d;
+}
+
diff --git a/src/orbits.h b/src/orbits.h
index a5e75dc..067aa95 100644
--- a/src/orbits.h
+++ b/src/orbits.h
@@ -5,6 +5,8 @@
 
 #include <glm/glm.hpp>
 
+#include "util.h"
+
 
 struct ellipse_parameters
 {
@@ -16,13 +18,6 @@ struct ellipse_parameters
 	glm::vec2 f2;
 };
 
-struct ellipse_3d
-{
-	ellipse_parameters params;
-	glm::vec3* vertices;
-	uint vert_count;
-};
-
 struct orbital_elements
 {
 	ellipse_parameters ep;
@@ -32,6 +27,14 @@ struct orbital_elements
 	double nu;		// NOTE: (ν) true anomaly
 };
 
+struct ellipse_3d
+{
+	ellipse_parameters params;
+	glm::vec3* vertices;
+	uint vert_count;
+};
+
+
 ellipse_parameters
 constructEllipse(double a, double b);
 
@@ -40,8 +43,7 @@ constructEllipse(double a, double b);
 ellipse_3d
 constructEllipse3D(ellipse_parameters ep, uint vert_count);
 
-
-/* NOTE: how-to propagate orbit position given inital true anomaly, semimajor
+/* NOTE: how-to propagate orbit position given initial true anomaly, semimajor
  *	axis, mean motion, and eccentricity: ref) section 4.4, Kepler's Problem,
  * 	"Space Flight Dynamics" by Craig A. Kluever
  *
@@ -75,48 +77,8 @@ constructEllipse3D(ellipse_parameters ep, uint vert_count);
  *			anomoly:
  *			tan(ϴ2/2) = sqrt((1+e) / (1-e)) * tan(E2/2)
  */
-
-#ifdef ORBIT_IMPLEMENTATION
-
 double
-getMeanMotion(ellipse_parameters ep)
-{
-	return -1;
-}
-
-ellipse_parameters
-constructEllipse(double a, double b)
-{
-	assert(a > 0 && b > 0 && a >= b);
-	ellipse_parameters ep = { a, b };
-	ep.c = sqrt(a * a - b * b);
-	ep.e = ep.c / ep.a;
-	ep.f1.x = -1 * ep.c;
-	ep.f2.x = ep.c;
-	return ep;
-}
-
-ellipse_3d
-constructEllipse3D(ellipse_parameters ep, uint vert_count)
-{
-	assert(ep.a > 0 && ep.b > 0 &&
-		ep.a >= ep.b &&
-		vert_count > 0);
-
-	ellipse_3d e3d = { ep, nullptr, vert_count};
-	e3d.vertices = UTIL_ALLOC(vert_count, glm::vec3);
-	double angle = 2 * M_PI / vert_count;
-
-	for (uint i = 0; i < vert_count; i++) {
-		double a = angle * i;
-		// NOTE: solving for distance in polar coordinates
-		double r = ep.b / sqrt(1 - (ep.e * pow(cos(a), 2)));
-		// NOTE: converting from polar to rectangular
-		e3d.vertices[i] = glm::vec3(cos(a) *  r, sin(a) * r, 0);
-	}
-
-	return e3d;
-}
-
+getPropagatedPosition(orbital_elements orbit,
+	double initial_anom,
+	unsigned int time_step);
 
-#endif
diff --git a/tests/orbit_test.cpp b/tests/orbit_test.cpp
new file mode 100644
index 0000000..7766212
--- /dev/null
+++ b/tests/orbit_test.cpp
@@ -0,0 +1,34 @@
+
+#include "dumbLog.h"
+
+#include "../src/orbits.h"
+
+
+void
+testOrbitPropagation()
+{
+	orbital_elements orbit = {};
+	ellipse_parameters ep = {};
+	ep.a = 26564.5; // meters
+	ep.e = 0.7411;
+	orbit.ep = ep;
+
+	// NOTE: solve for gravitational parameter from mean motion
+	double mean_motion = 0.00014582; // rad/s
+	orbit.mu = pow(mean_motion, 2) * pow(ep.a, 3);
+
+	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";
+
+	return;
+}
+
+int
+main()
+{
+	testOrbitPropagation();
+	return 0;
+}