diff --git a/scripts/precalc_omega_debug.py b/scripts/precalc_omega_debug.py new file mode 100644 index 0000000..4d63312 --- /dev/null +++ b/scripts/precalc_omega_debug.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +""" +Precalculate expected values for test_omega_debug.cpp refactoring. +Computes expected orbital elements after a prograde burn at apoapsis. +""" + +import math +import sys +sys.path.insert(0, "/home/agent/dev/claudes_game") +from scripts.sim_engine import * + + +def main(): + earth_mass = 5.972e24 + mu = G * earth_mass + + # Initial orbit: zero inclination, omega = 0, start at apoapsis (nu = pi) + elements = OrbitalElements( + a=1.0e7, + e=0.3, + nu=math.pi, + inc=1e-12, + Omega=0.0, + omega=0.0, + ) + + pos, vel = orbital_to_cartesian(elements, earth_mass) + + r = vmag(pos) + v = vmag(vel) + + print("// === Initial state ===") + print(f"// pos = ({pos[0]:.15e}, {pos[1]:.15e}, {pos[2]:.15e}) m") + print(f"// vel = ({vel[0]:.15e}, {vel[1]:.15e}, {vel[2]:.15e}) m/s") + print(f"// r = {r:.15e} m") + print(f"// v = {v:.15e} m/s") + print() + + # Eccentricity vector + r_dot_v = vdot(pos, vel) + e_vec = ( + ((v * v - mu / r) * pos[0] - r_dot_v * vel[0]) / mu, + ((v * v - mu / r) * pos[1] - r_dot_v * vel[1]) / mu, + ((v * v - mu / r) * pos[2] - r_dot_v * vel[2]) / mu, + ) + e_mag = vmag(e_vec) + + print(f"// e_vec_initial = ({e_vec[0]:.15e}, {e_vec[1]:.15e}, {e_vec[2]:.15e})") + print(f"// e_initial = {e_mag:.15e}") + print() + + # Apply prograde burn (1000 m/s) + burn_dir = get_burn_direction(BurnDirection.PROGRADE, pos, vel) + dv = 1000.0 + dv_vec = vscale(burn_dir, dv) + vel_new = vadd(vel, dv_vec) + + v_new = vmag(vel_new) + + print("// === After prograde burn ===") + print(f"// burn_dir = ({burn_dir[0]:.15e}, {burn_dir[1]:.15e}, {burn_dir[2]:.15e})") + print(f"// vel_new = ({vel_new[0]:.15e}, {vel_new[1]:.15e}, {vel_new[2]:.15e}) m/s") + print(f"// v_new = {v_new:.15e} m/s") + print() + + # Reconstruct orbital elements + new_elements = cartesian_to_orbital_elements(pos, vel_new, earth_mass) + + print(f"// new elements:") + print(f"// a = {new_elements.a:.15e} m") + print(f"// e = {new_elements.e:.15e}") + print(f"// nu = {new_elements.nu:.15e} rad ({math.degrees(new_elements.nu):.6f} deg)") + print(f"// inc = {new_elements.inc:.15e} rad") + print(f"// Omega = {new_elements.Omega:.15e} rad") + print(f"// omega = {new_elements.omega:.15e} rad ({math.degrees(new_elements.omega):.6f} deg)") + print() + + # New eccentricity vector + r_dot_v_new = vdot(pos, vel_new) + e_vec_new = ( + ((v_new * v_new - mu / r) * pos[0] - r_dot_v_new * vel_new[0]) / mu, + ((v_new * v_new - mu / r) * pos[1] - r_dot_v_new * vel_new[1]) / mu, + ((v_new * v_new - mu / r) * pos[2] - r_dot_v_new * vel_new[2]) / mu, + ) + + print(f"// e_vec_new = ({e_vec_new[0]:.15e}, {e_vec_new[1]:.15e}, {e_vec_new[2]:.15e})") + print(f"// e_new = {vmag(e_vec_new):.15e}") + print() + + # Verify omega is in [0, 2*pi) + print("// === Omega range check ===") + print(f"// omega = {new_elements.omega:.15e} rad") + print(f"// omega in [0, 2*pi)? {0.0 <= new_elements.omega < 2.0 * math.pi}") + print() + + # After a prograde burn at apoapsis (nu=pi), the eccentricity vector flips + # direction because the increased velocity raises the opposite side of the orbit. + # This means omega should change from 0 to approximately pi. + # + # Rationale: at apoapsis, position and velocity are perpendicular. + # A prograde burn adds velocity along the velocity direction, increasing energy. + # The eccentricity vector formula: e_vec = (v^2 - mu/r)*r/μ - (r·v)*v/μ + # At apoapsis: r·v = 0, so e_vec = (v^2 - mu/r) * r / mu + # After prograde burn, v increases, so (v^2 - mu/r) becomes more positive, + # making e_vec more aligned with r direction. + # Since at apoapsis, r points opposite to periapsis direction (for ω=0), + # the eccentricity vector flips, meaning periapsis moves to the opposite side, + # so ω → π. + + print("// === Expected test values ===") + print(f"// a_expected = {new_elements.a:.15e}") + print(f"// e_expected = {new_elements.e:.15e}") + print(f"// omega_expected = {new_elements.omega:.15e} rad ({math.degrees(new_elements.omega):.6f} deg)") + print() + + # Also compute expected values using the same check as the original test + print("// For WithinAbs assertions:") + print(f"// a := {new_elements.a:.15e}") + print(f"// e := {new_elements.e:.15e}") + print(f"// omega := {new_elements.omega:.15e}") + print(f"// inc := {new_elements.inc:.15e}") + print(f"// Omega := {new_elements.Omega:.15e}") + print(f"// nu := {new_elements.nu:.15e}") + print(f"// r := {r:.15e}") + print(f"// v_new := {v_new:.15e}") + + +if __name__ == "__main__": + main() diff --git a/tests/test_omega_debug.cpp b/tests/test_omega_debug.cpp new file mode 100644 index 0000000..f61a252 --- /dev/null +++ b/tests/test_omega_debug.cpp @@ -0,0 +1,63 @@ +#include +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/test_utilities.h" +#include + +using Catch::Matchers::WithinAbs; + +SCENARIO("Omega reconstruction after prograde burn at apoapsis", "[omega][debug]") { + const double parent_mass = 5.972e24; + const double mu = G * parent_mass; + + OrbitalElements elements = {}; + elements.semi_major_axis = 1.0e7; + elements.eccentricity = 0.3; + elements.true_anomaly = M_PI; + elements.inclination = 1e-12; + elements.longitude_of_ascending_node = 0.0; + elements.argument_of_periapsis = 0.0; + + Vec3 pos = {}, vel = {}; + orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); + + const double r = vec3_magnitude(pos); + const double v = vec3_magnitude(vel); + + const double r_dot_v = vec3_dot(pos, vel); + const Vec3 e_vec = { + ((v * v - mu / r) * pos.x - r_dot_v * vel.x) / mu, + ((v * v - mu / r) * pos.y - r_dot_v * vel.y) / mu, + ((v * v - mu / r) * pos.z - r_dot_v * vel.z) / mu, + }; + const double e_initial = vec3_magnitude(e_vec); + + SECTION("initial apoapsis state is correct") { + REQUIRE_THAT(r, WithinAbs(1.3e7, R_TOL)); + REQUIRE_THAT(v, WithinAbs(4632.763232589246, V_TOL)); + REQUIRE_THAT(e_initial, WithinAbs(0.3, E_TOL)); + REQUIRE_THAT(e_vec.x / e_initial, WithinAbs(1.0, E_TOL)); + REQUIRE_THAT(e_vec.y, WithinAbs(0.0, E_TOL)); + REQUIRE_THAT(e_vec.z, WithinAbs(0.0, E_TOL)); + } + + SECTION("prograde burn at apoapsis flips periapsis") { + const Vec3 burn_dir = vec3_normalize(vel); + const Vec3 delta_v = vec3_scale(burn_dir, 1000.0); + const Vec3 vel_new = vec3_add(vel, delta_v); + const double v_new = vec3_magnitude(vel_new); + + const OrbitalElements new_elements = cartesian_to_orbital_elements(pos, vel_new, parent_mass); + + REQUIRE_THAT(new_elements.semi_major_axis, WithinAbs(1.346885753127762e7, A_TOL)); + REQUIRE(new_elements.semi_major_axis > elements.semi_major_axis); + REQUIRE_THAT(new_elements.eccentricity, WithinAbs(3.481049006486453e-2, E_TOL)); + REQUIRE(new_elements.eccentricity < elements.eccentricity); + REQUIRE_THAT(new_elements.argument_of_periapsis, WithinAbs(M_PI, ANG_TOL)); + REQUIRE_THAT(new_elements.true_anomaly, WithinAbs(0.0, ANG_TOL)); + REQUIRE_THAT(new_elements.inclination, WithinAbs(0.0, ANG_TOL)); + REQUIRE_THAT(v_new, WithinAbs(5632.763232589246, V_TOL)); + REQUIRE(v_new > v); + } +}