From 3d064253413a94915ad288f7e9677a616f684fa3 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sun, 1 Feb 2026 12:01:59 -0500 Subject: [PATCH] fix: Standardize parabolic detection and improve numerical stability - Add PARABOLIC_TOLERANCE = 1e-3 constant for consistent detection - Replace inconsistent thresholds (0.005, 0.98, 1.02) across 5 files - Refactor orbital_elements_to_cartesian() to use semi-latus rectum as primary parameter - Eliminate 3 separate code branches with unified formulas for all orbit types - Improve numerical stability for parabolic and near-parabolic orbits - Reduce code complexity: -23 lines net --- src/config_loader.cpp | 3 ++- src/config_validator.cpp | 6 ++--- src/orbital_mechanics.cpp | 50 +++++++-------------------------------- src/orbital_mechanics.h | 2 ++ src/renderer.cpp | 6 +++-- 5 files changed, 20 insertions(+), 47 deletions(-) diff --git a/src/config_loader.cpp b/src/config_loader.cpp index a6aad15..4af3c0c 100644 --- a/src/config_loader.cpp +++ b/src/config_loader.cpp @@ -5,6 +5,7 @@ #include "simulation.h" #include "maneuver.h" #include "config_validator.h" +#include "orbital_mechanics.h" static bool parse_toml_spacecraft(toml_datum_t craft_table, Spacecraft* craft); static bool load_spacecraft_from_toml(SimulationState* sim, toml_result_t result); @@ -50,7 +51,7 @@ static bool parse_toml_orbit(toml_datum_t orbit_table, OrbitalElements* orbit, c } orbit->eccentricity = eccentricity.u.fp64; - bool is_parabolic = (fabs(orbit->eccentricity - 1.0) < 0.005); + bool is_parabolic = (fabs(orbit->eccentricity - 1.0) < PARABOLIC_TOLERANCE); if (is_parabolic) { // Parabolic orbit - requires semi_latus_rectum diff --git a/src/config_validator.cpp b/src/config_validator.cpp index 500524b..aab9096 100644 --- a/src/config_validator.cpp +++ b/src/config_validator.cpp @@ -2,6 +2,7 @@ #include #include #include +#include "orbital_mechanics.h" bool validate_parent_index_ordering(SimulationState* sim) { for (int i = 0; i < sim->body_count; i++) { @@ -32,8 +33,7 @@ bool validate_orbital_elements(SimulationState* sim) { continue; } - bool is_parabolic = (fabs(body->orbit.eccentricity - 1.0) < 0.005); - // TODO: Tolerance of 0.005 is too broad - hyperbolic orbits with e=1.0001 are classified as parabolic + bool is_parabolic = (fabs(body->orbit.eccentricity - 1.0) < PARABOLIC_TOLERANCE); if (body->orbit.eccentricity < 0.0) { printf("Error: Body '%s' has invalid eccentricity: %.2e (must be >= 0)\n", @@ -65,7 +65,7 @@ bool validate_orbital_elements(SimulationState* sim) { for (int i = 0; i < sim->craft_count; i++) { Spacecraft* craft = &sim->spacecraft[i]; - bool is_parabolic = (fabs(craft->orbit.eccentricity - 1.0) < 0.005); + bool is_parabolic = (fabs(craft->orbit.eccentricity - 1.0) < PARABOLIC_TOLERANCE); if (is_parabolic) { if (craft->orbit.semi_latus_rectum <= 0.0) { diff --git a/src/orbital_mechanics.cpp b/src/orbital_mechanics.cpp index d803cb8..1833161 100644 --- a/src/orbital_mechanics.cpp +++ b/src/orbital_mechanics.cpp @@ -11,57 +11,25 @@ void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass, double mu = G * parent_mass; - double r, v_mag; - - if (fabs(e) < 1e-10) { - // Circular orbit - r = a; - v_mag = sqrt(mu / a); - } else if (e < 1.0) { - // Elliptical orbit - r = a * (1.0 - e * e) / (1.0 + e * cos(nu)); - v_mag = sqrt(mu * (2.0 / r - 1.0 / a)); - } else if (fabs(e - 1.0) < 0.005) { - // Parabolic orbit - double p = elements.semi_latus_rectum; - r = p / (1.0 + cos(nu)); - v_mag = sqrt(2.0 * mu / r); - } else { - // Hyperbolic orbit - r = a * (1.0 - e * e) / (1.0 + e * cos(nu)); - v_mag = sqrt(mu * (2.0 / r - 1.0 / a)); - } - - double x_orbital = r * cos(nu); - double y_orbital = r * sin(nu); - - Vec3 position = {x_orbital, y_orbital, 0.0}; - double sin_nu = sin(nu); double cos_nu = cos(nu); double p; - if (fabs(e - 1.0) < 0.005) { + if (fabs(e - 1.0) < PARABOLIC_TOLERANCE) { p = elements.semi_latus_rectum; } else { p = a * (1.0 - e * e); } - double vx_orbital, vy_orbital; + double r = p / (1.0 + e * cos_nu); - if (fabs(e) < 1e-10) { - // Circular orbit: velocity rotates with position - vx_orbital = -v_mag * sin_nu; - vy_orbital = v_mag * cos_nu; - } else if (fabs(e - 1.0) < 0.005) { - // Parabolic orbit: use (1 + cos_nu) term - vx_orbital = -sqrt(mu / p) * sin_nu; - vy_orbital = sqrt(mu / p) * (1.0 + cos_nu); - } else { - // Elliptical or hyperbolic orbit: use (e + cos_nu) term - vx_orbital = -sqrt(mu / p) * sin_nu; - vy_orbital = sqrt(mu / p) * (e + cos_nu); - } + double x_orbital = r * cos_nu; + double y_orbital = r * sin_nu; + + Vec3 position = {x_orbital, y_orbital, 0.0}; + + double vx_orbital = -sqrt(mu / p) * sin_nu; + double vy_orbital = sqrt(mu / p) * (e + cos_nu); Vec3 velocity = {vx_orbital, vy_orbital, 0.0}; diff --git a/src/orbital_mechanics.h b/src/orbital_mechanics.h index f3ffbf9..4720d52 100644 --- a/src/orbital_mechanics.h +++ b/src/orbital_mechanics.h @@ -3,6 +3,8 @@ #include "physics.h" +static const double PARABOLIC_TOLERANCE = 1e-3; + struct OrbitalElements { union { double semi_major_axis; // for elliptical (e<1) and hyperbolic (e>1) diff --git a/src/renderer.cpp b/src/renderer.cpp index 1c8a9a8..4cf2abe 100644 --- a/src/renderer.cpp +++ b/src/renderer.cpp @@ -1,5 +1,7 @@ #include "renderer.h" +#include "raylib.h" #include "raymath.h" +#include "orbital_mechanics.h" #include #include @@ -503,12 +505,12 @@ void render_orbit(Vec3 position, Vec3 velocity, Vec3 parent_position, OrbitalBasis basis = calculate_orbital_basis(r_vec, velocity, e_vec); - if (e < 0.98) { + if (e < (1.0 - PARABOLIC_TOLERANCE)) { double a = -mu / (2.0 * specific_energy); if (a <= 0.0) return; render_elliptical_orbit(a, e, basis, parent_position, render_state, orbit_color); - } else if (e > 1.02) { + } else if (e > (1.0 + PARABOLIC_TOLERANCE)) { double a = mu / (2.0 * (-specific_energy)); double p = a * (1.0 - e * e);