From 11fc5fa5e6b715ad9291f0675e83eae4bc108a84 Mon Sep 17 00:00:00 2001
From: cinnaboot <cinnaboot@protonmail.com>
Date: Thu, 31 Mar 2022 12:21:18 -0400
Subject: [PATCH] make use of system_2body time_step property

---
 src/main.cpp         |  2 +-
 src/orbits.cpp       | 14 ++++++++------
 src/orbits.h         |  3 +--
 tests/orbit_test.cpp |  8 ++++----
 4 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index 7000685..3cb83da 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -163,7 +163,7 @@ loadScene(RenderState* rs)
 void
 updateSatelliteModel(const system_2body& sys, satellite& sat)
 {
-	sat.theta = getPropagatedTrueAnomaly(sys, sat.theta, sys.time_step);
+	sat.theta = getPropagatedTrueAnomaly(sys, sat.theta);
 	sat.gamma = orbitGetFlightPathAngle(sys.ep.e, sat.theta);
 	sat.r = orbitGetRadialDistance(sys.ep.e, sys.ep.p, sat.theta);
 	sat.v = orbitGetVelocity(sys.epsilon, sys.body.mu, sat.r);
diff --git a/src/orbits.cpp b/src/orbits.cpp
index 17848d4..9bece20 100644
--- a/src/orbits.cpp
+++ b/src/orbits.cpp
@@ -208,14 +208,13 @@ getNextTrialValue(double err, double ecc, double test_anom, double mean_anom)
 
 double
 getPropagatedEccAnomaly(system_2body sys,
-	double initial_anom,
-	unsigned int time_step)
+	double initial_anom)
 {
 	double e = sys.ep.e;
 	double E1 = getEccAnomFromTrueAnom(e, initial_anom);
 	double M1 = getMeanAnomFromEccAnom(E1, e);
 	double n = getMeanMotion(sys.body.mu, sys.elements.a);
-	double M2 = getPropagatedMeanAnom(M1, n, time_step);
+	double M2 = getPropagatedMeanAnom(M1, n, sys.time_step);
 	double E2_1 = getInitialTrialValue(M2, e);
 
 	// test if guess is a solution to kepler's equation
@@ -236,10 +235,13 @@ getPropagatedEccAnomaly(system_2body sys,
 
 double
 getPropagatedTrueAnomaly(system_2body sys,
-	double initial_anom,
-	unsigned int time_step)
+	double initial_anom)
 {
-	double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom, time_step);
+	// NOTE: 'pause' simulation when time_step is set close to 0
+	if (sys.time_step < 1e-8)
+		return initial_anom;
+
+	double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom);
 	return getTrueAnomFromEccAnom(sys.ep.e, ecc_anom);
 }
 
diff --git a/src/orbits.h b/src/orbits.h
index 6410068..ee5129a 100644
--- a/src/orbits.h
+++ b/src/orbits.h
@@ -145,8 +145,7 @@ orbitUpdate(orbital_elements& o, double a, double e);
  */
 double
 getPropagatedTrueAnomaly(system_2body sys,
-	double initial_anom,
-	unsigned int time_step);
+	double initial_anom);
 
 // NOTE: aka) the trajectory equation (eq. 2.45)
 double
diff --git a/tests/orbit_test.cpp b/tests/orbit_test.cpp
index 8247377..007fbde 100644
--- a/tests/orbit_test.cpp
+++ b/tests/orbit_test.cpp
@@ -24,8 +24,8 @@ TEST_CASE("orbit propagation", "[orbits]")
 	double mu = 398601.68;	// NOTE: gravitational parameter
 	double r = 6378;		// NOTE: body radius in km
 	double initial_anom = 260 * M_PI / 180; // NOTE: radians
-	unsigned int time_step = 60 * 50; // NOTE: seconds
 	system_2body sys = systemInit(gravBodyInit(mu, r), orbitInit(a, e));
+	sys.time_step = 60 * 50; // NOTE: seconds
 
 	double E1 = getEccAnomFromTrueAnom(sys.ep.e, initial_anom);
 	REQUIRE_THAT(E1, WithinAbs(-0.8615, 1e-4));
@@ -36,17 +36,17 @@ TEST_CASE("orbit propagation", "[orbits]")
 	double n = getMeanMotion(sys.body.mu, sys.ep.a);
 	REQUIRE_THAT(n, WithinAbs(0.00014582, 1e-8));
 
-	double M2 = getPropagatedMeanAnom(M1, n, time_step);
+	double M2 = getPropagatedMeanAnom(M1, n, sys.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, sys.ep.e);
 	REQUIRE_THAT(E2_1, WithinAbs(0.315452, 1e-5));
 
-	double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom, time_step);
+	double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom);
 	REQUIRE_THAT(ecc_anom, WithinAbs(0.481518, 1e-5));
 
-	double true_anom = getPropagatedTrueAnomaly(sys, initial_anom, time_step);
+	double true_anom = getPropagatedTrueAnomaly(sys, initial_anom);
 	REQUIRE_THAT(true_anom, WithinAbs(1.1339, 1e-4));
 
 	double r2 = orbitGetRadialDistance(sys.ep.e, sys.ep.p, true_anom);