From d52a5b1d693912e3c251e091407b9189a5969552 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Fri, 26 Mar 2021 15:48:44 -0400 Subject: [PATCH] add specific energy, flight path angle, and speed functions --- src/gooey.cpp | 6 ++-- src/main.cpp | 32 +++++++---------- src/orbits.cpp | 93 +++++++++++++++++++++++++++++++------------------- src/orbits.h | 27 +++++++++++---- 4 files changed, 95 insertions(+), 63 deletions(-) diff --git a/src/gooey.cpp b/src/gooey.cpp index 93381af..59f8aa9 100644 --- a/src/gooey.cpp +++ b/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(); diff --git a/src/main.cpp b/src/main.cpp index 6680d9c..46a8978 100644 --- a/src/main.cpp +++ b/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 diff --git a/src/orbits.cpp b/src/orbits.cpp index 7ccadb7..a55f522 100644 --- a/src/orbits.cpp +++ b/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)); -} - diff --git a/src/orbits.h b/src/orbits.h index 851b1ad..3dc2a8e 100644 --- a/src/orbits.h +++ b/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);