implement orbit model layout changes

src/main.cpp

@@ -16,9 +16,7 @@
 
 
 const double SCALING = 0.001;
-static orbital_elements g_orbit = {};
-static ellipse_3d g_ellipse3d;
-static uint g_time_step = 100;
+static system_2body g_sys = {0};
 
 
 simple_mesh*
@@ -64,31 +62,34 @@ createSatelliteMesh()
 	return sm;
 }
 
+// TODO: this should be 2 separate functions. 1 for updating the g_sys.sat true
+//	anomaly, and another for updating the renderer entity
 void
-updateSatellitePosition(entity& satellite)
-{ const static glm::mat4 xform =
+updateSatellitePosition(entity& sat)
+{
+	// TODO: decouple framerate from time_step
+	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);
-	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));
+	satellite& s1 = g_sys.sat;
+	s1.theta = getPropagatedTrueAnomaly(g_sys, s1.theta, g_sys.time_step);
+	s1.gamma = orbitGetFlightPathAngle(g_sys.ep.e, s1.theta);
+	s1.r = orbitGetRadialDistance(g_sys.ep.e, g_sys.ep.p, s1.theta);
+	s1.v = orbitGetVelocity(g_sys.epsilon, g_sys.body.mu, s1.r);
+	glm::vec2 coords = polarToRect(s1.theta, s1.r);
+	glm::vec3 v = s1.position = glm::vec3(coords, 0);
+	entSetWorldPosition(sat, xform * glm::vec4(v.x, v.y, v.z, 1));
 }
 
-// NOTE: use ellipseValidate(orbit.ep) before calling to avoid failing
-//	assertions
+// NOTE: use ellipseValidate(ep) before calling to avoid failing assertions
 void
-updateOrbit(orbital_elements& orbit, ellipse_3d& e3d, entity& ellipse_entity)
+updateOrbit(system_2body sys, entity& ellipse_entity)
 {
-	orbit.ep = ellipseInitAE(orbit.ep.a, orbit.ep.e);
-	ellipse3DUpdate(orbit.ep, e3d);
+	//sys.ep = ellipseInitAE(sys.ep.a, sys.ep.e);
+	ellipse3DUpdate(sys.ep, sys.e3d);
 
-	for (uint i = 0; i < e3d.vert_count; i++)
-		ellipse_entity.mesh->vertices[i] = e3d.vertices[i];
+	for (uint i = 0; i < sys.e3d.vert_count; i++)
+		ellipse_entity.mesh->vertices[i] = sys.e3d.vertices[i];
 
 	entUpdateSimpleMesh(ellipse_entity, ellipse_entity.mesh, GL_LINE_LOOP);
 }
@@ -108,19 +109,21 @@ preFrameCallback(render_state* rs)
 
 	if (is.window_closed || is.escape)
 		rs->running = false;
-
-	if (g_time_step > 0)
-		updateSatellitePosition(rs->render_groups[0].entities[2]);
 }
 
 void
 postFrameCallback(render_state* rs)
 {
-	static orbital_elements c_orb = {};
-	c_orb = g_orbit;
-	gooDraw(rs->handles->window, c_orb, g_time_step);
-	bool update_if_paused = false;
+	//bool update_if_paused = false;
+	//static orbital_elements c_orb = {};
+	//c_orb = g_orbit;
+	gooDraw(rs->handles->window, g_sys);
+	g_sys.ep = ellipseInitAE(g_sys.ep.a, g_sys.ep.e);
+	entity& ellipse_entity = rs->render_groups[0].entities[0];
+	updateOrbit(g_sys, ellipse_entity);
+	updateSatellitePosition(rs->render_groups[0].entities[2]);
 
+#if 0
 	if (c_orb.nu != g_orbit.nu &&
 		((c_orb.nu >= 0 && c_orb.nu <= M_PI) ||
 		(c_orb.nu < 0 && c_orb.nu >= -1 * M_PI)))
@@ -144,7 +147,7 @@ postFrameCallback(render_state* rs)
 			c_orb.ep = ellipseInitAE(c_orb.ep.a, c_orb.ep.e);
 			g_orbit.ep = c_orb.ep;
 			entity& ellipse_entity = rs->render_groups[0].entities[0];
-			updateOrbit(g_orbit, g_ellipse3d, ellipse_entity);
+			updateOrbit(g_sys, ellipse_entity);
 
 			// TODO: would be nice to have an interface for
 			//	updating/constructing an orbit
@@ -154,8 +157,9 @@ postFrameCallback(render_state* rs)
 		}
 	}
 
-	if (g_time_step == 0 && update_if_paused)
+	if (g_sys.time_step == 0 && update_if_paused)
 		updateSatellitePosition(rs->render_groups[0].entities[2]);
+#endif
 }
 
 int
@@ -184,9 +188,11 @@ main()
 		dims.x / dims.y
 	);
 
-	// NOTE: semi-major axis in km, eccentricity, gravitational parameter
-	g_orbit = orbitInit(26564.5, 0.7411, 398601.68);
-	ellipse_parameters ep = ellipseInitAB(g_orbit.ep.a, g_orbit.ep.b);
+	double a = 26564.5;		// NOTE: semi-major axis in km
+	double e = 0.7411;		// NOTE: eccentricity
+	double mu = 398601.68;	// NOTE: gravitational parameter
+	double r = 6378;		// NOTE: body radius in km
+	g_sys = systemInit(gravBodyInit(mu, r), orbitInit(a, e));
 
 	// TODO: (renderer) this needs to be more convenient
 	shader_wrapper sw = { SIMPLE_SHADER, nullptr, rs->simple_shader };
@@ -194,8 +200,8 @@ main()
 	rs->render_group_count = 1;
 
 	entity& ellipse_entity = rs->render_groups[0].entities[0];
-	g_ellipse3d = ellipseInit3D(ep, 256);
-	simple_mesh* sm = constructEllipseMesh(g_ellipse3d, glm::vec3(255, 0, 255));
+	g_sys.e3d = ellipseInit3D(g_sys.ep, 256);
+	simple_mesh* sm = constructEllipseMesh(g_sys.e3d, glm::vec3(255, 0, 255));
 	entInitMesh(ellipse_entity, sm, GL_LINE_LOOP);
 	entRotate(ellipse_entity, (float) M_PI_2, glm::vec3(1, 0, 0));
 
@@ -208,7 +214,7 @@ main()
 	sm = createSatelliteMesh();
 	entInitMesh(satellite_entity, sm, GL_TRIANGLES);
 	entRotate(satellite_entity, (float) M_PI_2, glm::vec3(1, 0, 0));
-	entSetWorldPosition(satellite_entity, g_ellipse3d.vertices[0]);
+	entSetWorldPosition(satellite_entity, g_sys.e3d.vertices[0]);
 
 	renDoRenderLoop(rs, 60 , preFrameCallback, postFrameCallback);
 

src/orbits.cpp

@@ -2,6 +2,41 @@
 #include "orbits.h"
 
 
+//
+// WIP
+//
+const static uint ELLIPSE_VERT_COUNT = 256;
+const static uint DEFAULT_TIME_STEP = 100;
+
+
+grav_body
+gravBodyInit(double mu, double r)
+{
+	grav_body gb = {0};
+	gb.mu = mu;
+	gb.radius = r;
+	return gb;
+}
+
+system_2body
+systemInit(grav_body gb, orbital_elements el)
+{
+	system_2body s = {0};
+	s.body = gb;
+	s.elements = el;
+	s.ep = ellipseInitAE(el.a, el.e);
+	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.time_step = DEFAULT_TIME_STEP;
+
+	return s;
+}
+//
+// WIP
+//
+
+
 ellipse_parameters
 ellipseInitAB(double a, double b)
 {
@@ -51,7 +86,7 @@ ellipseInit3D(ellipse_parameters ep, uint vert_count)
 		ep.a >= ep.b &&
 		vert_count > 0);
 
-	ellipse_3d e3d = { ep, nullptr, vert_count};
+	ellipse_3d e3d = { nullptr, vert_count};
 	// TODO: need to free this allocation at some point
 	e3d.vertices = UTIL_ALLOC(vert_count, glm::vec3);
 	ellipse3DUpdate(ep, e3d);
@@ -59,24 +94,18 @@ ellipseInit3D(ellipse_parameters ep, uint vert_count)
 }
 
 orbital_elements
-orbitInit(double a, double e, double mu)
+orbitInit(double a, double e)
 {
-	orbital_elements o = {};
-	o.mu = mu;
-	orbitUpdate(o, a, e);
+	// TODO: remaining elements: iota, ohm, omega, nu
+	orbital_elements o = {0};
+	o.a = a;
+	o.e = e;
+	o.nu = 0;
 	return o;
 }
 
-void
-orbitUpdate(orbital_elements& o, double a, double e)
-{
-	o.ep = ellipseInitAE(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)
+orbitGetRadialDistance(double e, double p, double true_anom)
 {
 	return p / (1 + e * cos(true_anom));
 }
@@ -183,14 +212,14 @@ getNextTrialValue(double err, double ecc, double test_anom, double mean_anom)
 }
 
 double
-getPropagatedEccAnomaly(orbital_elements orbit,
+getPropagatedEccAnomaly(system_2body sys,
 	double initial_anom,
 	unsigned int time_step)
 {
-	double e = orbit.ep.e;
+	double e = sys.ep.e;
 	double E1 = getEccAnomFromTrueAnom(e, initial_anom);
 	double M1 = getMeanAnomFromEccAnom(E1, e);
-	double n = getMeanMotion(orbit.mu, orbit.ep.a);
+	double n = getMeanMotion(sys.body.mu, sys.elements.a);
 	double M2 = getPropagatedMeanAnom(M1, n, time_step);
 	double E2_1 = getInitialTrialValue(M2, e);
 
@@ -211,11 +240,11 @@ getPropagatedEccAnomaly(orbital_elements orbit,
 }
 
 double
-getPropagatedTrueAnomaly(orbital_elements orbit,
+getPropagatedTrueAnomaly(system_2body sys,
 	double initial_anom,
 	unsigned int time_step)
 {
-	double ecc_anom = getPropagatedEccAnomaly(orbit, initial_anom, time_step);
-	return getTrueAnomFromEccAnom(orbit.ep.e, ecc_anom);
+	double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom, time_step);
+	return getTrueAnomFromEccAnom(sys.ep.e, ecc_anom);
 }
 

src/orbits.h

@@ -23,40 +23,21 @@ struct ellipse_parameters
 struct orbital_elements
 {
 	// NOTE: classical orbital elements
-	// TODO: would be better to have a higher level composite object so we
-	// 	don't have to store this twice
-	ellipse_parameters ep;
-	// semimajor axis, a
-	// eccentricity, e
+	double a;		// NOTE: semimajor axis
+	double e;		// NOTE: eccentricity
 	double iota;	// NOTE: (ι) inclination
+	double ohm;		// NOTE: (Ω) longitude of the ascending node
 	double omega;	// NOTE: (ω) argument of periapsis
-	double nu;		// NOTE: (ν) true anomaly
-	// longitude of the ascending node, Ω
-	// argument of periapsis, ω
-	// true anomaly at epoch, θ0
-
-#if 1
-	double mu;		// NOTE: (μ) gravitational parameter
-	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 v;		// NOTE: velocity magnitute
-	double gamma;	// NOTE: (γ) flight path angle
-#endif
+	double nu;		// NOTE: (ν) true anomaly at T0
 };
 
 struct ellipse_3d
 {
-	// TODO: would be better to have a higher level composite object so we
-	// 	don't have to store this twice
-	ellipse_parameters ep;
 	glm::vec3* vertices;
 	uint vert_count;
 };
 
-struct gravitational_body
+struct grav_body
 {
 	double mu;		// NOTE: (μ) gravitational parameter
 	double radius;	// NOTE: radius of ideal sphere representing the body
@@ -67,18 +48,18 @@ struct satellite
 {
 	glm::vec3 position;
 	glm::vec3 velocity;
-	// double theta;	// NOTE: true anomaly
-	// double r;		// NOTE: radius magnitude
-	// double gamma;	// NOTE: (γ) flight path angle
-	// double v;		// NOTE: velocity magnitute
+	double theta;		// NOTE: true anomaly
+	double r;			// NOTE: radius magnitude at theta
+	double gamma;		// NOTE: (γ) flight path angle
+	double v;			// NOTE: velocity magnitute
 };
 
 // NOTE: top level composite for 2 body system
-struct system
+struct system_2body
 {
-	gravitational_body grav_body;
+	grav_body body;
 	satellite sat;
-	//ellipse_parameters ep;
+	ellipse_parameters ep;
 	ellipse_3d e3d;
 	orbital_elements elements;
 
@@ -86,9 +67,24 @@ struct system
 	double h;			// NOTE: angular momentum
 	double r_apoapsis;	// NOTE: apoapse distance from body center
 	double r_periapsis;	// NOTE: periapsis distance from body center
+
+	double time_step;
 };
 
 
+
+//
+// WIP
+//
+grav_body
+gravBodyInit(double mu, double r);
+
+system_2body
+systemInit(grav_body gb, orbital_elements el);
+//
+// WIP
+//
+
 ellipse_parameters
 ellipseInitAB(double a, double b);
 
@@ -110,7 +106,7 @@ void
 ellipse3DUpdate(ellipse_parameters ep, ellipse_3d& e3d);
 
 orbital_elements
-orbitInit(double a, double e, double mu);
+orbitInit(double a, double e);
 
 double
 orbitGetVelocity(double epsilon, double mu, double r);
@@ -156,13 +152,13 @@ orbitUpdate(orbital_elements& o, double a, double e);
  *			tan(ϴ2/2) = sqrt((1+e) / (1-e)) * tan(E2/2)
  */
 double
-getPropagatedTrueAnomaly(orbital_elements orbit,
+getPropagatedTrueAnomaly(system_2body sys,
 	double initial_anom,
 	unsigned int time_step);
 
 // NOTE: aka) the trajectory equation (eq. 2.45)
 double
-getRadialDistance(double e, double p, double true_anom);
+orbitGetRadialDistance(double e, double p, double true_anom);
 
 double
 orbitGetAngularMomentum(double p, double mu);