diff --git a/src/orbits.cpp b/src/orbits.cpp
index fbf848a..7c9b48b 100644
--- a/src/orbits.cpp
+++ b/src/orbits.cpp
@@ -263,3 +263,27 @@ orbitClampAngle(double theta)
 	double revs = floor(theta / (2 * M_PI));
 	return theta - revs * 2 * M_PI;
 }
+
+double
+orbitGetTimeOfFlight(system_2body sys, double theta_0, double theta_1)
+{
+	double e = sys.ep.e;
+
+	// get mean motion (n)
+	double n = getMeanMotion(sys.body.mu, sys.ep.a);
+
+	// get mean anomaly for theta_0 (M0)
+	double ecc_1 = getEccAnomFromTrueAnom(sys.ep.e, theta_0);
+
+	// get mean anomaly for theta_1 (M1)
+	double ecc_2 = getEccAnomFromTrueAnom(sys.ep.e, theta_1);
+
+	// test if orbit passes through periapsis, and if so, add 2_PI
+	// 		TODO: could also check for multiple revolutions, but that would be
+	// 			pretty unlikely if we're updating the simulation every 33ms
+	// use Kepler's equation for time of flight
+	double tof_1 = 1 / n * (ecc_1 - e * sin(ecc_1));
+	double tof_2 = 1 / n * (ecc_2 - e * sin(ecc_2));
+
+	return tof_2 - tof_1;
+}
diff --git a/tests/orbit_test.cpp b/tests/orbit_test.cpp
index ac778d8..576a6b5 100644
--- a/tests/orbit_test.cpp
+++ b/tests/orbit_test.cpp
@@ -19,10 +19,10 @@ TEST_CASE("orbit construction", "[orbits]")
 // NOTE: example 4.6 in "Space Flight Dynamics" by Craig A. Kluever
 TEST_CASE("orbit propagation", "[orbits]")
 {
-	double a = 26564.5;		// NOTE: semi-major axis in km
-	double e = 0.7411;		// NOTE: eccentricity
-	double mu = 398601.68;	// NOTE: gravitational parameter
-	double r = 6378;		// NOTE: body radius in km
+	double a = 26564.5;		// semi-major axis in km
+	double e = 0.7411;		// eccentricity
+	double mu = 398601.68;	// gravitational parameter
+	double r = 6378;		// body radius in km
 	double initial_anom = 260 * M_PI / 180; // NOTE: radians
 	double time_step = 60 * 50; // NOTE: seconds
 	system_2body sys = systemInit(gravBodyInit(mu, r), orbitInit(a, e));