add specific energy, flight path angle, and speed functions

src/gooey.cpp

@@ -91,8 +91,10 @@ gooDraw(SDL_Window* window, orbital_elements& orbit, uint& time_step)
 	ImGui::Text("	x: %f km", orbit.pos.x);
 	ImGui::Text("	y: %f km", orbit.pos.y);
 	ImGui::Text("	z: %f km", orbit.pos.z);
-	ImGui::Text("radial velocity:");
-	ImGui::Text("flight path angle:");
+	ImGui::Text("velocity, km/s:             %f", orbit.v);
+	ImGui::Text("flight path angle:          %f", orbit.gamma);
+	ImGui::Text("epsilon, spec. orb. energy: %f", orbit.epsilon);
+	ImGui::Text("h, angular momentum:        %f", orbit.h);
 
 	ImGui::End();
 	ImGui::Render();

src/main.cpp

@@ -64,28 +64,17 @@ createSatelliteMesh()
 	return sm;
 }
 
-orbital_elements
-initOrbit()
-{
-	orbital_elements orbit = {};
-
-	double a = 26564.5; // kilometers
-	double e = 0.7411;
-	orbit.ep = constructEllipseAE(a, e);
-	orbit.mu = 398601.68; // NOTE: earth's gravitational parameter in km^3/s^2
-	return orbit;
-}
-
 void
 updateSatellitePosition(entity& satellite)
-{
-	const static glm::mat4 xform =
+{ const static glm::mat4 xform =
 		glm::rotate(glm::mat4(1.0), (float) M_PI_2, glm::vec3(1, 0, 0));
 
 	// TODO: decouple framerate from time_step
 	g_orbit.nu = getPropagatedTrueAnomaly(g_orbit, g_orbit.nu, g_time_step);
-	double r2 = getRadialDistance(g_orbit.ep.e, g_orbit.ep.p, g_orbit.nu);
-	glm::vec2 coords = polarToRect(g_orbit.nu, r2);
+	g_orbit.gamma = orbitGetFlightPathAngle(g_orbit.ep.e, g_orbit.nu);
+	g_orbit.r = getRadialDistance(g_orbit.ep.e, g_orbit.ep.p, g_orbit.nu);
+	g_orbit.v = orbitGetVelocity(g_orbit.epsilon, g_orbit.mu, g_orbit.r);
+	glm::vec2 coords = polarToRect(g_orbit.nu, g_orbit.r);
 	glm::vec3 v = g_orbit.pos = glm::vec3(coords, 0);
 	entSetWorldPosition(satellite, xform * glm::vec4(v.x, v.y, v.z, 1));
 }
@@ -127,7 +116,7 @@ void
 postFrameCallback(render_state* rs)
 {
 	static orbital_elements c_orb = {};
-	orbitCopy(g_orbit, c_orb);
+	c_orb = g_orbit;
 	gooDraw(rs->handles->window, c_orb, g_time_step);
 
 	// TODO: will need more validation from GUI as added
@@ -142,9 +131,14 @@ postFrameCallback(render_state* rs)
 
 		if (validateEllipse(c_orb.ep)) {
 			c_orb.ep = constructEllipseAE(c_orb.ep.a, c_orb.ep.e);
-			ellipseCopy(c_orb.ep, g_orbit.ep);
+			g_orbit.ep = c_orb.ep;
 			entity& ellipse_entity = rs->render_groups[0].entities[0];
 			updateOrbit(g_orbit, g_ellipse3d, ellipse_entity);
+
+			// TODO: would be nice to have an interface for
+			//	updating/constructing an orbit
+			g_orbit.epsilon = orbitGetSpecificEnergy(g_orbit.ep.a, g_orbit.mu);
+			g_orbit.h = orbitGetAngularMomentum(g_orbit.ep.p, g_orbit.mu);
 		}
 	}
 
@@ -182,7 +176,7 @@ main()
 		dims.x / dims.y
 	);
 
-	g_orbit = initOrbit();
+	g_orbit = orbitInit(26564.5, 0.7411, 398601.68);
 	ellipse_parameters ep = constructEllipseAB(g_orbit.ep.a, g_orbit.ep.b);
 
 	// TODO: (renderer) this needs to be more convenient

src/orbits.cpp

@@ -44,29 +44,6 @@ ellipsesEqual(ellipse_parameters& e1, ellipse_parameters& e2)
 		e1.e == e2.e);
 }
 
-void
-ellipseCopy(const ellipse_parameters& e1, ellipse_parameters& e2)
-{
-	e2.a = e1.a;
-	e2.b = e1.b;
-	e2.e = e1.e;
-	e2.c = e1.c;
-	e2.p = e1.p;
-	e2.f1 = e1.f1;
-	e2.f2 = e1.f2;
-}
-
-void
-orbitCopy(const orbital_elements& o1, orbital_elements& o2)
-{
-	ellipseCopy(o1.ep, o2.ep);
-	o2.iota = o1.iota;
-	o2.omega = o1.omega;
-	o2.mu = o1.mu;
-	o2.nu = o1.nu;
-	o2.pos = o1.pos;
-}
-
 ellipse_3d
 constructEllipse3D(ellipse_parameters ep, uint vert_count)
 {
@@ -81,6 +58,59 @@ constructEllipse3D(ellipse_parameters ep, uint vert_count)
 	return e3d;
 }
 
+orbital_elements
+orbitInit(double a, double e, double mu)
+{
+	orbital_elements o = {};
+	o.mu = mu;
+	orbitUpdate(o, a, e);
+	return o;
+}
+
+void
+orbitUpdate(orbital_elements& o, double a, double e)
+{
+	o.ep = constructEllipseAE(a, e);
+	o.epsilon = orbitGetSpecificEnergy(o.ep.a, o.mu);
+	o.h = orbitGetAngularMomentum(o.ep.p, o.mu);
+}
+
+double
+getRadialDistance(double e, double p, double true_anom)
+{
+	return p / (1 + e * cos(true_anom));
+}
+
+double
+orbitGetAngularMomentum(double p, double mu)
+{
+	return sqrt(mu * p);
+}
+
+double
+orbitGetSpecificEnergy(double a, double mu)
+{
+	return -1 * mu / (2 * a);
+}
+
+double
+orbitGetVelocity(double epsilon, double mu, double r)
+{
+	return sqrt(2 * (epsilon + mu / r));
+}
+
+double
+orbitGetFlightPathAngle(double e, double true_anom)
+{
+	return atan(e * sin(true_anom) / (1 + e * cos(true_anom)));
+}
+
+glm::vec2
+polarToRect(double angle, double r)
+{
+	return glm::vec2(r * cos(angle), r * sin(angle));
+}
+
 void
 ellipse3DUpdate(ellipse_parameters ep, ellipse_3d& e3d)
 {
@@ -94,6 +124,11 @@ ellipse3DUpdate(ellipse_parameters ep, ellipse_3d& e3d)
 	}
 }
 
+
+//
+// NOTE: propagate anomaly functions:
+//
+
 double
 getEccAnomFromTrueAnom(double ecc, double true_anom)
 {
@@ -184,15 +219,3 @@ getPropagatedTrueAnomaly(orbital_elements orbit,
 	return getTrueAnomFromEccAnom(orbit.ep.e, ecc_anom);
 }
 
-double
-getRadialDistance(double e, double p, double true_anom)
-{
-	return p / (1 + e * cos(true_anom));
-}
-
-glm::vec2
-polarToRect(double angle, double r)
-{
-	return glm::vec2(r * cos(angle), r * sin(angle));
-}
-

src/orbits.h

@@ -30,8 +30,9 @@ struct orbital_elements
 	double epsilon;	// NOTE: (ε) specific orbital energy, MJ/kg
 	double h;		// NOTE: angular momentum
 
+	double r;		// NOTE: radial distance
 	glm::vec3 pos;	// NOTE: true anomaly in 3d with ep.f2 at origin
-	double vr;		// NOTE: radial velocity
+	double v;		// NOTE: velocity magnitute
 	double gamma;	// NOTE: (γ) flight path angle
 };
 
@@ -55,12 +56,6 @@ validateEllipse(const ellipse_parameters& ep);
 bool
 ellipsesEqual(ellipse_parameters& e1, ellipse_parameters& e2);
 
-void
-ellipseCopy(const ellipse_parameters& e1, ellipse_parameters& e2);
-
-void
-orbitCopy(const orbital_elements& o1, orbital_elements& o2);
-
 // NOTE: create vertices for a 3d ellipse
 // NOTE: all vertices are in the x/y plane with z = 0
 ellipse_3d
@@ -69,6 +64,18 @@ constructEllipse3D(ellipse_parameters ep, uint vert_count);
 void
 ellipse3DUpdate(ellipse_parameters ep, ellipse_3d& e3d);
 
+orbital_elements
+orbitInit(double a, double e, double mu);
+
+double
+orbitGetVelocity(double epsilon, double mu, double r);
+
+double
+orbitGetFlightPathAngle(double e, double true_anom);
+
+void
+orbitUpdate(orbital_elements& o, double a, double e);
+
 /* NOTE: how-to propagate orbit position given initial true anomaly, semimajor
  *	axis, mean motion, and eccentricity: ref) section 4.4, Kepler's Problem,
  * 	"Space Flight Dynamics" by Craig A. Kluever
@@ -112,6 +119,12 @@ getPropagatedTrueAnomaly(orbital_elements orbit,
 double
 getRadialDistance(double e, double p, double true_anom);
 
+double
+orbitGetAngularMomentum(double p, double mu);
+
+double
+orbitGetSpecificEnergy(double a, double mu);
+
 glm::vec2
 polarToRect(double true_anom, double r);