From d9546be2a1763484c768300dfcaa9e40115111cf Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sat, 19 Jun 2021 10:23:56 -0400 Subject: [PATCH] implement orbit model layout changes --- src/main.cpp | 76 +++++++++++++++++++++++++++----------------------- src/orbits.cpp | 69 ++++++++++++++++++++++++++++++++------------- src/orbits.h | 64 ++++++++++++++++++++---------------------- 3 files changed, 120 insertions(+), 89 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index eba061b..5949723 100644 --- a/src/main.cpp +++ b/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); diff --git a/src/orbits.cpp b/src/orbits.cpp index cb45c2a..bc72bc9 100644 --- a/src/orbits.cpp +++ b/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); } diff --git a/src/orbits.h b/src/orbits.h index 22cdc12..3a8231e 100644 --- a/src/orbits.h +++ b/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);