From 9b3137d27bbe9704516e6469aba6d2c1fe84ccc8 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Wed, 3 Mar 2021 16:41:43 -0500 Subject: [PATCH] update demo to use correct orbit propagation --- src/main.cpp | 71 ++++++++++++++++++++++++++++++-------------- src/orbits.cpp | 23 ++++++++++---- src/orbits.h | 11 ++++++- tests/orbit_test.cpp | 26 ++++++++++++---- 4 files changed, 96 insertions(+), 35 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 2bdc81e..006b6c0 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,8 +13,9 @@ #include "orbits.h" +const double SCALING = 0.001; +static orbital_elements g_orbit = {}; static ellipse_3d g_ellipse; -static uint g_index = 0; simple_mesh* @@ -35,6 +36,7 @@ constructCircleMesh(float r, uint vert_count, glm::vec3 vertex_color) { simple_mesh* sm = meInitMesh(vert_count); float angle = (float) M_PI * 2 / vert_count; + r = r / SCALING; for (uint i = 0; i < vert_count; i++) { sm->vertices[i] = glm::vec3(cos(i * angle) * r, sin(i * angle) * r, 0); @@ -44,21 +46,50 @@ constructCircleMesh(float r, uint vert_count, glm::vec3 vertex_color) return sm; } +simple_mesh* +createSatelliteMesh() +{ + simple_mesh* sm = meInitMesh(3); + glm::vec3 s_color(255, 0, 0); + sm->vertices[0] = glm::vec3(0, 1 / SCALING, 0); + sm->vert_colors[0] = s_color; + sm->vertices[1] = glm::vec3(-1 / SCALING, -1 / SCALING, 0); + sm->vert_colors[1] = s_color; + sm->vertices[2] = glm::vec3(1 / SCALING, -1 / SCALING, 0); + sm->vert_colors[2] = s_color; + + return sm; +} + +orbital_elements +initOrbit() +{ + orbital_elements orbit = {}; + + double a = 26564.5; // kilometers + double e = 0.7411; + double b = a * sqrt(1 - pow(e, 2.0)); + orbit.ep = constructEllipse(a, b); + orbit.mu = 398601.68; // NOTE: earth's gravitational parameter in km^3/s^2 + + return orbit; +} + void doFrameCallback(render_state* rs) { + const static glm::mat4 xform = + glm::rotate(glm::mat4(1.0), (float) M_PI_2, glm::vec3(1, 0, 0)); entity& satellite = rs->render_groups[0].entities[2]; - g_index++; - if (g_index == g_ellipse.vert_count) g_index = 0; - static glm::mat4 xform = - glm::rotate(glm::mat4(1.0), (float) M_PI_2, glm::vec3(1, 0, 0)); - glm::vec3 v = g_ellipse.vertices[g_index]; - entSetWorldPosition(satellite, xform * glm::vec4(v.x, v.y, v.z, 1)); + // TODO: decouple framerate from time_step + unsigned int time_step = 100; // NOTE: seconds -#if 0 - LOG(Debug) << "SDL_GetTicks(): " << SDL_GetTicks() << "ms\n"; -#endif + g_orbit.nu = getPropagatedTrueAnomaly(g_orbit, g_orbit.nu, time_step); + double r2 = getRadialPosition(g_orbit.ep, g_orbit.nu); + glm::vec2 coords = polarToRect(g_orbit.nu, r2); + glm::vec3 v = glm::vec3(coords, 0); + entSetWorldPosition(satellite, xform * glm::vec4(v.x, v.y, v.z, 1)); } int @@ -75,44 +106,40 @@ main() cameraInitPerspective( rs->cam, - glm::vec3(0, -100, 0), + glm::vec3(0, -75 / SCALING, 0), glm::vec3(0, 0, 0), glm::vec3(0,0,1) ); + g_orbit = initOrbit(); + ellipse_parameters ep = constructEllipse(g_orbit.ep.a, g_orbit.ep.b); + // TODO: (renderer) this needs to be more convenient shader_wrapper sw = { SIMPLE_SHADER, nullptr, rs->simple_shader }; rs->render_groups = renAllocateGroup(3, sw); rs->render_group_count = 1; entity& ellipse_entity = rs->render_groups[0].entities[0]; - ellipse_parameters ep = constructEllipse(50, 30); g_ellipse = constructEllipse3D(ep, 256); simple_mesh* sm = constructEllipseMesh(g_ellipse, glm::vec3(255, 0, 255)); - entInitMesh(ellipse_entity, sm, GL_LINES); + entInitMesh(ellipse_entity, sm, GL_LINE_LOOP); entRotate(ellipse_entity, (float) M_PI_2, glm::vec3(1, 0, 0)); entity& planet_entity = rs->render_groups[0].entities[1]; sm = constructCircleMesh(3, 24, glm::vec3(255, 255, 0)); entInitMesh(planet_entity, sm, GL_LINE_LOOP); entRotate(planet_entity, (float) M_PI_2, glm::vec3(1, 0, 0)); - entTranslate(planet_entity, glm::vec3(ep.f2.x, 0, 0)); entity& satellite_entity = rs->render_groups[0].entities[2]; - sm = meInitMesh(3); - glm::vec3 s_color(255, 0, 0); - sm->vertices[0] = glm::vec3(0, 1, 0); - sm->vert_colors[0] = s_color; - sm->vertices[1] = glm::vec3(-1, -1, 0); - sm->vert_colors[1] = s_color; - sm->vertices[2] = glm::vec3(1, -1, 0); - sm->vert_colors[2] = s_color; + sm = createSatelliteMesh(); entInitMesh(satellite_entity, sm, GL_TRIANGLES); entRotate(satellite_entity, (float) M_PI_2, glm::vec3(1, 0, 0)); entSetWorldPosition(satellite_entity, g_ellipse.vertices[0]); renDoRenderLoop(rs, 60, doFrameCallback); + // TODO: clean up mesh pointers? don't remember if renderer does that + // automatically renShutdown(rs); return 0; } diff --git a/src/orbits.cpp b/src/orbits.cpp index 43ca544..6b428c3 100644 --- a/src/orbits.cpp +++ b/src/orbits.cpp @@ -11,10 +11,10 @@ constructEllipse(double a, double b) ellipse_parameters ep = { a, b }; ep.c = sqrt(a * a - b * b); ep.e = ep.c / ep.a; + ep.p = ep.a * (1 - pow(ep.e, 2)); ep.f1.x = -1 * ep.c; ep.f2.x = ep.c; - return ep; -} + return ep; } ellipse_3d constructEllipse3D(ellipse_parameters ep, uint vert_count) @@ -29,10 +29,9 @@ constructEllipse3D(ellipse_parameters ep, uint vert_count) for (uint i = 0; i < vert_count; i++) { double a = angle * i; - // NOTE: solving for distance in polar coordinates - double r = ep.b / sqrt(1 - (ep.e * pow(cos(a), 2))); - // NOTE: converting from polar to rectangular - e3d.vertices[i] = glm::vec3(cos(a) * r, sin(a) * r, 0); + // NOTE: solving for distance in polar coordinates relative to focus + double r = ep.a * (1 - pow(ep.e, 2)) / (1 + ep.e * cos(a)); + e3d.vertices[i] = glm::vec3(polarToRect(a, r), 0); } return e3d; @@ -128,3 +127,15 @@ getPropagatedTrueAnomaly(orbital_elements orbit, return getTrueAnomFromEccAnom(orbit.ep.e, ecc_anom); } +double +getRadialPosition(ellipse_parameters ep, double true_anom) +{ + return ep.p / (1 + ep.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 70b9745..29c635b 100644 --- a/src/orbits.h +++ b/src/orbits.h @@ -14,6 +14,7 @@ struct ellipse_parameters double b; // NOTE: semi-minor axis double e; // NOTE: eccentricity double c; // NOTE: linear eccentricity + double p; // NOTE: semilatus rectum glm::vec2 f1; glm::vec2 f2; }; @@ -29,7 +30,7 @@ struct orbital_elements struct ellipse_3d { - ellipse_parameters params; + ellipse_parameters ep; glm::vec3* vertices; uint vert_count; }; @@ -82,3 +83,11 @@ getPropagatedTrueAnomaly(orbital_elements orbit, double initial_anom, unsigned int time_step); +// NOTE: aka) the trajectory equation (eq. 2.45) +// NOTE: returns radial distance in kilometers +double +getRadialPosition(ellipse_parameters ep, double true_anom); + +glm::vec2 +polarToRect(double true_anom, double r); + diff --git a/tests/orbit_test.cpp b/tests/orbit_test.cpp index ccfc4b6..3d8eacc 100644 --- a/tests/orbit_test.cpp +++ b/tests/orbit_test.cpp @@ -1,7 +1,7 @@ #include -#include "dumbLog.h" +#include #include "../src/orbits.cpp" @@ -11,21 +11,28 @@ using Catch::Matchers::WithinAbs; +TEST_CASE("orbit construction", "[orbits]") +{ + +} + // NOTE: example 4.6 in "Space Flight Dynamics" by Craig A. Kluever TEST_CASE("orbit propagation", "[orbits]") { orbital_elements orbit = {}; - ellipse_parameters ep = {}; - ep.a = 26564.5; // meters - ep.e = 0.7411; + // TODO: might want another contructor for ellipse_parameters + double a = 26564.5; // kilometers + double e = 0.7411; + double b = a * sqrt(1 - pow(e, 2.0)); + ellipse_parameters ep = constructEllipse(a, b); orbit.ep = ep; + double initial_anom = 260 * M_PI / 180; // radians + unsigned int time_step = 60 * 50; // seconds // NOTE: solve for gravitational parameter from mean motion double mean_motion = 0.00014582; // rad/s orbit.mu = pow(mean_motion, 2) * pow(ep.a, 3); - double initial_anom = 260 * M_PI / 180; // radians - unsigned int time_step = 60 * 50; // seconds double E1 = getEccAnomFromTrueAnom(ep.e, initial_anom); REQUIRE_THAT(E1, WithinAbs(-0.8615, 1e-4)); @@ -48,5 +55,12 @@ TEST_CASE("orbit propagation", "[orbits]") double true_anom = getPropagatedTrueAnomaly(orbit, initial_anom, time_step); REQUIRE_THAT(true_anom, WithinAbs(1.1339, 1e-4)); + + double r2 = getRadialPosition(ep, true_anom); + REQUIRE_THAT(r2, WithinAbs(9116.1, 0.1)); + + glm::vec2 pos = polarToRect(true_anom, r2); + REQUIRE_THAT(pos.x, WithinAbs(3856.9, 0.1)); + REQUIRE_THAT(pos.y, WithinAbs(8259.9, 0.1)); }