From 1923e7a296d76752932da0215fa23574ea46370f Mon Sep 17 00:00:00 2001
From: cinnaboot <cinnaboot@protonmail.com>
Date: Fri, 1 Apr 2022 11:26:09 -0400
Subject: [PATCH] add orbital period property, and test case

---
 src/gooey.cpp        |  1 +
 src/orbits.cpp       |  7 +++++++
 src/orbits.h         | 12 ++++++++----
 tests/orbit_test.cpp |  9 +++++++++
 4 files changed, 25 insertions(+), 4 deletions(-)

diff --git a/src/gooey.cpp b/src/gooey.cpp
index 5a396c9..9aeea30 100644
--- a/src/gooey.cpp
+++ b/src/gooey.cpp
@@ -139,6 +139,7 @@ void drawSystemWindow(system_2body& sys, bool& running)
 		ImGui::SameLine();
 		ImGui::InputDouble("##", &sys.time_step, 5);
 
+		ImGui::Text("orbital period:             %f", sys.orbital_period);
 		ImGui::Text("epsilon, spec. orb. energy: %f", sys.epsilon);
 		ImGui::Text("h, angular momentum:        %f", sys.h);
 		ImGui::Text("r_apoapsis:                 %f", sys.r_apoapsis);
diff --git a/src/orbits.cpp b/src/orbits.cpp
index daf2ba1..45ec5c0 100644
--- a/src/orbits.cpp
+++ b/src/orbits.cpp
@@ -16,6 +16,7 @@ systemInit(grav_body gb, orbital_elements el)
 	s.e3d = ellipseInit3D(s.ep, ELLIPSE_VERT_COUNT);
 	s.epsilon = orbitGetSpecificEnergy(s.ep.a, gb.mu);
 	s.h = orbitGetAngularMomentum(s.ep.p, gb.mu);
+	s.orbital_period = orbitGetPeriod(s.ep.a, gb.mu);
 	s.r_apoapsis = s.ep.a - s.ep.c;
 	s.r_periapsis = 2 * s.ep.a - s.r_apoapsis;
 	s.time_step = DEFAULT_TIME_STEP;
@@ -129,6 +130,12 @@ orbitGetFlightPathAngle(double e, double true_anom)
 	return atan(e * sin(true_anom) / (1 + e * cos(true_anom)));
 }
 
+double
+orbitGetPeriod(double a, double mu)
+{
+	return 2 * M_PI / sqrt(mu) * sqrt(pow(a, 3));
+}
+
 glm::vec2
 polarToRect(double angle, double r)
 {
diff --git a/src/orbits.h b/src/orbits.h
index ee5129a..293f0a0 100644
--- a/src/orbits.h
+++ b/src/orbits.h
@@ -62,10 +62,11 @@ struct system_2body
 	ellipse_3d e3d;
 	orbital_elements elements;
 
-	double epsilon;		// NOTE: (ε) specific orbital energy, MJ/kg
-	double h;			// NOTE: angular momentum
-	double r_apoapsis;	// NOTE: apoapse distance from body center
-	double r_periapsis;	// NOTE: periapsis distance from body center
+	double epsilon;			// NOTE: (ε) specific orbital energy, MJ/kg
+	double h;				// NOTE: angular momentum
+	double r_apoapsis;		// NOTE: apoapse distance from body center
+	double r_periapsis;		// NOTE: periapsis distance from body center
+	double orbital_period;	// NOTE: in seconds
 
 	double time_step;
 };
@@ -106,6 +107,9 @@ orbitGetVelocity(double epsilon, double mu, double r);
 double
 orbitGetFlightPathAngle(double e, double true_anom);
 
+double
+orbitGetPeriod(double a, double mu);
+
 void
 orbitUpdate(orbital_elements& o, double a, double e);
 
diff --git a/tests/orbit_test.cpp b/tests/orbit_test.cpp
index 007fbde..b04355c 100644
--- a/tests/orbit_test.cpp
+++ b/tests/orbit_test.cpp
@@ -57,3 +57,12 @@ TEST_CASE("orbit propagation", "[orbits]")
 	REQUIRE_THAT(pos.y, WithinAbs(8259.9, 0.1));
 }
 
+// NOTE: example 2.5 - c
+TEST_CASE("orbital period", "[orbits]")
+{
+	double a = 24371;		// NOTE: semi-major axis in km
+	double mu = 398601.68;	// NOTE: gravitational parameter
+
+	double T = orbitGetPeriod(a, mu);
+	REQUIRE_THAT(T, WithinAbs(37863.5, 50e-3));
+}