From 909a7c8ae17e2a816bc92fdaec227c2672bddab3 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Mon, 13 Mar 2023 11:26:40 -0400 Subject: [PATCH] fix the source of all the 'nan' position errors We were getting slight floating point errors from glm::dot() before sending the result to acos(). Anything out of the range of -1 >= x <= 1 is undefined --- src/game.cpp | 4 +++- src/main.cpp | 7 ++++--- src/orbits.cpp | 11 +++++++++-- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/game.cpp b/src/game.cpp index 5c469b0..c82be5d 100644 --- a/src/game.cpp +++ b/src/game.cpp @@ -247,9 +247,11 @@ applyManeuver(GameOrbit* orbit, ManeuverNode* maneuver) // FIXME: debugging code +#if 0 GameOrbit old_orbit = *orbit; if (!old_orbit.in_use) exit(1); // don't optimize out - +#endif + // FIXME: debugging code dvec3 pos = orbitGetPositionVector(r, theta); diff --git a/src/main.cpp b/src/main.cpp index 2dab369..ea3117a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -267,11 +267,12 @@ loadRandomOrbits(GameState* gs, double semi_major_axis = fmod(fabs((double) rand()), semi_major_axis_range) + min_semi_major_axis; double eccentricity = max_eccentricity * rand() / RAND_MAX; double nu = max_nu * rand() / RAND_MAX; - //double omega = 2 * M_PI * rand() / RAND_MAX; - double omega = 0; + double omega = 2 * M_PI * rand() / RAND_MAX; + double iota = 2 * M_PI * rand() / RAND_MAX; + double ohm = 2 * M_PI * rand() / RAND_MAX; GameOrbit* gorb = loadOrbit(gs, rs, body, - orbitInit(semi_major_axis, eccentricity, 0, 0, omega, nu), + orbitInit(semi_major_axis, eccentricity, iota, ohm, omega, nu), rg, "???"); addManeuver(gs, gorb, ImpulseType::CIRCULARIZE_LOWERING); diff --git a/src/orbits.cpp b/src/orbits.cpp index ad5a6f2..2b1e981 100644 --- a/src/orbits.cpp +++ b/src/orbits.cpp @@ -111,12 +111,19 @@ orbitGetElementsFromStateVectors(glm::dvec3 r, glm::dvec3 v, double mu) if (ecc_v.z < 0) el.omega = 2 * M_PI - el.omega; // quadrant check } - // FIXME: I think this is the smoking gun for the source of all the nan - // position and velocity components if (el.e == 0) { el.nu = 0; } else { double cos_theta = glm::dot(r, ecc_v) / (el.e * r_mag); + + // NOTE: clamp to vaild range for acos function + // we were getting slight floating point errors out of the glm::dot() + // function above + if (cos_theta > 1.0) + cos_theta = 1.0; + else if (cos_theta < - 1.0) + cos_theta = -1.0; + el.nu= acos(cos_theta); }