initial work for impulsive maneuvers

src/game.cpp

@@ -164,8 +164,68 @@ removeManeuver(GameOrbit* orbit)
 	orbit->maneuver = {0};
 }
 
+// NOTE: test if the maneuver node would have occured between 2 orbit positions
+// 	eg) between 2 time steps
+bool
+testManeuverStep(const ManeuverNode& maneuver,
+		double previous_theta,
+		double next_theta)
+{
+	if (!maneuver.active)
+		return false;
+
+	double r = maneuver.true_anomaly;
+	previous_theta = orbitClampAngle(previous_theta);
+	next_theta = orbitClampAngle(next_theta);
+
+	// FIXME: 2 unhandled cases:
+	// 	1) sat passes through 0 degrees
+	// 	2) direction is ccw
+	//
+	if (previous_theta <= r && next_theta >= r)
+		return true;
+
+	return false;
 }
 
+void
+applyManeuver(GameOrbit* orbit)
+{
+	LOGF(Debug, "applying maneuver\n");
+
+	// FIXME: assuming prograde impulse vector
+	assert(orbit->maneuver.impulse_type == ManeuverType::PROGRADE);
+
+	// re-calculate state vectors at the maneuver node
+	TwoBodySystem& sys = orbit->system;
+	double theta = orbit->maneuver.true_anomaly;
+	double mu = sys.body.mu;
+	double r = orbitGetRadialDistance(sys.ep.e, sys.ep.p, theta);
+	double v = orbitGetVelocity(sys.epsilon, mu, r);
+	double gamma = orbitGetFlightPathAngle(sys.ep.e, theta);
+
+	// calculate new satellite state vectors
+	v += orbit->maneuver.impulse_delta_v; // prograde only
+
+	// get new angular momentum and specific orbital energy
+	double h = orbitGetAngularMomentumFromStateVectors(r, v, gamma);
+	double epsilon = orbitGetSpecificEnergyFromStateVectors(r, v, mu);
+
+	// update orbit pararmeters
+	double a = orbitGetSemiMajorAxis(epsilon, mu);
+	double p = orbitGetSemiLatusRectum(h, mu);
+	double e = ellipseGetEccentricity(a, p);
+	systemInit(orbit->system, sys.body, orbitInit(a, e));
+
+	// update ellipse3D & GLBuffer vertices
+	ellipse3DUpdate(sys.ep, orbit->e3d);
+	GLBuffer* buf = &orbit->ellipse_entity->meshes[0].vertex_attrib_buffers[0];
+	assert(utilCStrMatch(buf->name, "position"));
+	assert(buf->data_size == orbit->e3d.vert_count * sizeof(vec3));
+	updateGLBuffer(buf, orbit->e3d.vertices);
+
+	removeManeuver(orbit);
+}
 
 // internal
 

src/game.h

@@ -19,10 +19,21 @@ struct Ellipse3D
 	uint vert_count;
 };
 
+enum struct ManeuverType
+{
+	PROGRADE,
+	RETROGRADE,
+	INCLINATION_CHANGE,
+	GENERAL_PLANE_CHANGE,
+	OTHER,
+	VECTOR_TYPE_COUNT
+};
+
 struct ManeuverNode
 {
 	double true_anomaly;
 	vec3 impulse_vector;
+	ManeuverType impulse_type;
 	double impulse_delta_v;
 	// TODO: thrust calculations
 
@@ -79,3 +90,9 @@ bool addManeuver(GameOrbit* orbit, double anomaly, vec3 impulse_vector,
 		double impulse_delta_v);
 
 void removeManeuver(GameOrbit* orbit);
+
+bool testManeuverStep(const ManeuverNode& maneuver,
+		double previous_theta,
+		double next_theta);
+
+void applyManeuver(GameOrbit* orbit);

src/main.cpp

@@ -203,7 +203,8 @@ loadScene(GameState* gs, RenderState* rs)
 
 	double a = 26564.5;		// semi-major axis in km
 	double e = 0.7411;		// eccentricity
-	loadOrbit(gs, rs, body, orbitInit(a, e), rg, "sat_01");
+	GameOrbit* go_1 = loadOrbit(gs, rs, body, orbitInit(a, e), rg, "sat_01");
+	addManeuver(go_1, M_PI / 4, vec3(), 0.1);
 
 	double a_2 = r + 10000;
 	double e_2 = 0;

src/orbits.cpp

@@ -165,14 +165,6 @@ orbitGetPropagatedTrueAnomaly(TwoBodySystem sys,
 
 // FIXME: organize into interface/internal functions
 
-// NOTE: return an equivalent angle between 0 and 2 pi radians
-inline double
-orbitClampAngle(double theta)
-{
-	double revs = floor(theta / (2 * M_PI));
-	return theta - revs * 2 * M_PI;
-}
-
 double
 orbitTimeSincePeriapsis(TwoBodySystem sys, double theta);
 

src/orbits.h

@@ -40,7 +40,7 @@ struct GravBody
 struct Satellite
 {
 	glm::vec3 position;
-	glm::vec3 velocity;
+	// glm::vec3 velocity;
 	double theta;		// true anomaly
 	double r;			// radius magnitude at theta
 	double gamma;		// (γ) flight path angle
@@ -165,6 +165,14 @@ orbitGetRadialDistance(double e, double p, double true_anom)
 	return p / (1 + e * cos(true_anom));
 }
 
+// NOTE: return an equivalent angle between 0 and 2 pi radians
+inline double
+orbitClampAngle(double theta)
+{
+	double revs = floor(theta / (2 * M_PI));
+	return theta - revs * 2 * M_PI;
+}
+
 inline glm::vec2
 polarToRect(double angle, double r)
 {