update demo to use correct orbit propagation

src/main.cpp

@@ -13,8 +13,9 @@
 #include "orbits.h"
 
 
+const double SCALING = 0.001;
+static orbital_elements g_orbit = {};
 static ellipse_3d g_ellipse;
-static uint g_index = 0;
 
 
 simple_mesh*
@@ -35,6 +36,7 @@ constructCircleMesh(float r, uint vert_count, glm::vec3 vertex_color)
 {
 	simple_mesh* sm = meInitMesh(vert_count);
 	float angle = (float) M_PI * 2 / vert_count;
+	r = r / SCALING;
 
 	for (uint i = 0; i < vert_count; i++) {
 		sm->vertices[i] = glm::vec3(cos(i * angle) * r, sin(i * angle) * r, 0);
@@ -44,21 +46,50 @@ constructCircleMesh(float r, uint vert_count, glm::vec3 vertex_color)
 	return sm;
 }
 
+simple_mesh*
+createSatelliteMesh()
+{
+	simple_mesh* sm = meInitMesh(3);
+	glm::vec3 s_color(255, 0, 0);
+	sm->vertices[0] = glm::vec3(0, 1 / SCALING, 0);
+	sm->vert_colors[0] = s_color;
+	sm->vertices[1] = glm::vec3(-1 / SCALING, -1 / SCALING, 0);
+	sm->vert_colors[1] = s_color;
+	sm->vertices[2] = glm::vec3(1 / SCALING, -1 / SCALING, 0);
+	sm->vert_colors[2] = s_color;
+
+	return sm;
+}
+
+orbital_elements
+initOrbit()
+{
+	orbital_elements orbit = {};
+
+	double a = 26564.5; // kilometers
+	double e = 0.7411;
+	double b = a * sqrt(1 - pow(e, 2.0));
+	orbit.ep = constructEllipse(a, b);
+	orbit.mu = 398601.68; // NOTE: earth's gravitational parameter in km^3/s^2
+
+	return orbit;
+}
+
 void
 doFrameCallback(render_state* rs)
 {
+	const static glm::mat4 xform =
+		glm::rotate(glm::mat4(1.0), (float) M_PI_2, glm::vec3(1, 0, 0));
 	entity& satellite = rs->render_groups[0].entities[2];
 
-	g_index++;
-	if (g_index == g_ellipse.vert_count) g_index = 0;
-	static glm::mat4 xform =
-		glm::rotate(glm::mat4(1.0), (float) M_PI_2, glm::vec3(1, 0, 0));
-	glm::vec3 v = g_ellipse.vertices[g_index];
-	entSetWorldPosition(satellite, xform * glm::vec4(v.x, v.y, v.z, 1));
+	// TODO: decouple framerate from time_step
+	unsigned int time_step = 100; // NOTE: seconds
 
-#if 0
-	LOG(Debug) << "SDL_GetTicks(): " << SDL_GetTicks() << "ms\n";
-#endif
+	g_orbit.nu = getPropagatedTrueAnomaly(g_orbit, g_orbit.nu, time_step);
+	double r2 = getRadialPosition(g_orbit.ep, g_orbit.nu);
+	glm::vec2 coords = polarToRect(g_orbit.nu, r2);
+	glm::vec3 v = glm::vec3(coords, 0);
+	entSetWorldPosition(satellite, xform * glm::vec4(v.x, v.y, v.z, 1));
 }
 
 int
@@ -75,44 +106,40 @@ main()
 
 	cameraInitPerspective(
 		rs->cam,
-		glm::vec3(0, -100, 0),
+		glm::vec3(0, -75 / SCALING, 0),
 		glm::vec3(0, 0, 0),
 		glm::vec3(0,0,1)
 	);
 
+	g_orbit = initOrbit();
+	ellipse_parameters ep = constructEllipse(g_orbit.ep.a, g_orbit.ep.b);
+
 	// TODO: (renderer) this needs to be more convenient
 	shader_wrapper sw = { SIMPLE_SHADER, nullptr, rs->simple_shader };
 	rs->render_groups = renAllocateGroup(3, sw);
 	rs->render_group_count = 1;
 
 	entity& ellipse_entity = rs->render_groups[0].entities[0];
-	ellipse_parameters ep = constructEllipse(50, 30);
 	g_ellipse = constructEllipse3D(ep, 256);
 	simple_mesh* sm = constructEllipseMesh(g_ellipse, glm::vec3(255, 0, 255));
-	entInitMesh(ellipse_entity, sm, GL_LINES);
+	entInitMesh(ellipse_entity, sm, GL_LINE_LOOP);
 	entRotate(ellipse_entity, (float) M_PI_2, glm::vec3(1, 0, 0));
 
 	entity& planet_entity = rs->render_groups[0].entities[1];
 	sm = constructCircleMesh(3, 24, glm::vec3(255, 255, 0));
 	entInitMesh(planet_entity, sm, GL_LINE_LOOP);
 	entRotate(planet_entity, (float) M_PI_2, glm::vec3(1, 0, 0));
-	entTranslate(planet_entity, glm::vec3(ep.f2.x, 0, 0));
 
 	entity& satellite_entity = rs->render_groups[0].entities[2];
-	sm = meInitMesh(3);
-	glm::vec3 s_color(255, 0, 0);
-	sm->vertices[0] = glm::vec3(0, 1, 0);
-	sm->vert_colors[0] = s_color;
-	sm->vertices[1] = glm::vec3(-1, -1, 0);
-	sm->vert_colors[1] = s_color;
-	sm->vertices[2] = glm::vec3(1, -1, 0);
-	sm->vert_colors[2] = s_color;
+	sm = createSatelliteMesh();
 	entInitMesh(satellite_entity, sm, GL_TRIANGLES);
 	entRotate(satellite_entity, (float) M_PI_2, glm::vec3(1, 0, 0));
 	entSetWorldPosition(satellite_entity, g_ellipse.vertices[0]);
 
 	renDoRenderLoop(rs, 60, doFrameCallback);
 
+	// TODO: clean up mesh pointers? don't remember if renderer does that
+	// 	automatically
 	renShutdown(rs);
 	return 0;
 }

src/orbits.cpp

@@ -11,10 +11,10 @@ constructEllipse(double a, double b)
 	ellipse_parameters ep = { a, b };
 	ep.c = sqrt(a * a - b * b);
 	ep.e = ep.c / ep.a;
+	ep.p = ep.a * (1 - pow(ep.e, 2));
 	ep.f1.x = -1 * ep.c;
 	ep.f2.x = ep.c;
-	return ep;
-}
+	return ep; }
 
 ellipse_3d
 constructEllipse3D(ellipse_parameters ep, uint vert_count)
@@ -29,10 +29,9 @@ constructEllipse3D(ellipse_parameters ep, uint 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);
+		// NOTE: solving for distance in polar coordinates relative to focus
+		double r = ep.a * (1 - pow(ep.e, 2)) / (1 + ep.e * cos(a));
+		e3d.vertices[i] = glm::vec3(polarToRect(a, r), 0);
 	}
 
 	return e3d;
@@ -128,3 +127,15 @@ getPropagatedTrueAnomaly(orbital_elements orbit,
 	return getTrueAnomFromEccAnom(orbit.ep.e, ecc_anom);
 }
 
+double
+getRadialPosition(ellipse_parameters ep, double true_anom)
+{
+	return ep.p / (1 + ep.e * cos(true_anom));
+}
+
+glm::vec2
+polarToRect(double angle, double r)
+{
+	return glm::vec2(r * cos(angle), r * sin(angle));
+}
+

src/orbits.h

@@ -14,6 +14,7 @@ struct ellipse_parameters
 	double b;	// NOTE: semi-minor axis
 	double e;	// NOTE: eccentricity
 	double c;	// NOTE: linear eccentricity
+	double p;	// NOTE: semilatus rectum
 	glm::vec2 f1;
 	glm::vec2 f2;
 };
@@ -29,7 +30,7 @@ struct orbital_elements
 
 struct ellipse_3d
 {
-	ellipse_parameters params;
+	ellipse_parameters ep;
 	glm::vec3* vertices;
 	uint vert_count;
 };
@@ -82,3 +83,11 @@ getPropagatedTrueAnomaly(orbital_elements orbit,
 	double initial_anom,
 	unsigned int time_step);
 
+// NOTE: aka) the trajectory equation (eq. 2.45)
+// NOTE: returns radial distance in kilometers
+double
+getRadialPosition(ellipse_parameters ep, double true_anom);
+
+glm::vec2
+polarToRect(double true_anom, double r);
+

tests/orbit_test.cpp

@@ -1,7 +1,7 @@
 
 #include <catch2/catch.hpp>
 
-#include "dumbLog.h"
+#include <glm/glm.hpp>
 
 #include "../src/orbits.cpp"
 
@@ -11,21 +11,28 @@
 using Catch::Matchers::WithinAbs;
 
 
+TEST_CASE("orbit construction", "[orbits]")
+{
+
+}
+
 // NOTE: example 4.6 in "Space Flight Dynamics" by Craig A. Kluever
 TEST_CASE("orbit propagation", "[orbits]")
 {
 	orbital_elements orbit = {};
-	ellipse_parameters ep = {};
-	ep.a = 26564.5; // meters
-	ep.e = 0.7411;
+	// TODO: might want another contructor for ellipse_parameters
+	double a = 26564.5; // kilometers
+	double e = 0.7411;
+	double b = a * sqrt(1 - pow(e, 2.0));
+	ellipse_parameters ep = constructEllipse(a, b);
 	orbit.ep = ep;
+	double initial_anom = 260 * M_PI / 180; // radians
+	unsigned int time_step = 60 * 50; // seconds
 
 	// 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 E1 = getEccAnomFromTrueAnom(ep.e, initial_anom);
 	REQUIRE_THAT(E1, WithinAbs(-0.8615, 1e-4));
@@ -48,5 +55,12 @@ TEST_CASE("orbit propagation", "[orbits]")
 
 	double true_anom = getPropagatedTrueAnomaly(orbit, initial_anom, time_step);
 	REQUIRE_THAT(true_anom, WithinAbs(1.1339, 1e-4));
+
+	double r2 = getRadialPosition(ep, true_anom);
+	REQUIRE_THAT(r2, WithinAbs(9116.1, 0.1));
+
+	glm::vec2 pos = polarToRect(true_anom, r2);
+	REQUIRE_THAT(pos.x, WithinAbs(3856.9, 0.1));
+	REQUIRE_THAT(pos.y, WithinAbs(8259.9, 0.1));
 }