update orbitGetTimeofFlight for orbits that pass through periapsis

src/orbits.cpp

@@ -273,21 +273,21 @@ orbitGetTimeOfFlight(system_2body sys, double theta_begin, double theta_end)
 {
 	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);
+	double ecc_begin = getEccAnomFromTrueAnom(sys.ep.e, theta_begin);
+	double ecc_end = getEccAnomFromTrueAnom(sys.ep.e, theta_end);
 
-	// get mean anomaly for theta_1 (M1)
-	double ecc_2 = getEccAnomFromTrueAnom(sys.ep.e, theta_1);
+	// NOTE: test if flight passes through perisapsis
+	if (ecc_begin > ecc_end)
+		ecc_end += 2 * M_PI;
 
-	// 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));
+	double M1 = getMeanAnomFromEccAnom(ecc_begin, e);
+	double M2 = getMeanAnomFromEccAnom(ecc_end, e);
 
-	return tof_2 - tof_1;
+	// NOTE: Kepler's equation for time of flight
+	double tof_begin = 1 / n * M1;
+	double tof_end = 1 / n * M2;
+
+	return tof_end - tof_begin;
 }