diff --git a/continue.md b/continue.md index a31cd49..6a47fd0 100644 --- a/continue.md +++ b/continue.md @@ -7,6 +7,8 @@ - **Use SCENARIO to share setup/teardown across multiple SECTIONs.** Catch2 re-initializes the fixture before each SECTION, so declare shared constants, structs, and variables in the SCENARIO body (between the opening `{` and the first `SECTION`). These persist across all SECTIONs within the SCENARIO. - Example: a `SimulationState* sim` created once in the SCENARIO body, then each SECTION mutates and tests it independently. - Embed expected values directly in `WithinAbs()` calls (see Section 4 for precalc script usage). No need to declare named constants unless the value is reused. +- **No decorative comments.** Do not add `// (Old: ...)` comments, `===` separators, `---` separators, or any other decorative annotations. The SECTION description string is the documentation. +- **Use `REQUIRE()` for integer comparisons**, `WithinAbs()` only for floating-point. E.g., `REQUIRE(sim->body_count == 2)` not `REQUIRE_THAT(sim->body_count, WithinAbs(2.0, 0.001))`. ### 2. Duplication Elimination - Use lambdas that capture the fixture for repeated setup→call→assert patterns @@ -47,6 +49,7 @@ All constants defined in `src/test_utilities.h` — use those, do not redefine l - Global distances are dominated by parent body positions (e.g., Earth-Sun distance swamps LEO orbit). - **Always output SI units** (meters, m/s, seconds) — C++ tests use SI internally. - Output C++-style comments with precalculated expected values for embedding in the test. Tolerances are chosen separately by the test writer using the Tolerance Reference table — the precalc script should not output tolerance values. +- **No decorative comments in precalc scripts.** Use simple blank lines between sections, no `# ====` or `# -----` separators. - Run with: `python3 scripts/precalc_.py` - If sim_engine.py lacks a feature, **stop to notify the user what feature is missing** diff --git a/scripts/precalc_hybrid_burns.py b/scripts/precalc_hybrid_burns.py new file mode 100644 index 0000000..9f69ce3 --- /dev/null +++ b/scripts/precalc_hybrid_burns.py @@ -0,0 +1,507 @@ +#!/usr/bin/env python3 +""" +Precalculate expected values for test_hybrid_burns.cpp refactoring. +Uses sim_engine.py for physics propagation. +""" + +import math +import sys +sys.path.insert(0, "/home/agent/dev/claudes_game") +from scripts.sim_engine import * + + +def simulate_continuous_burn(initial_orbit, parent_mass, total_dv, burn_duration, + num_steps, direction): + """Simulate continuous/low-thrust burn with sub-steps.""" + current_orbit = initial_orbit + dt_burn_step = burn_duration / num_steps + dv_per_step = total_dv / num_steps + + for _ in range(num_steps): + pos, vel = orbital_to_cartesian(current_orbit, parent_mass) + burn_dir = get_burn_direction(direction, pos, vel) + dv_vec = vscale(burn_dir, dv_per_step) + vel = vadd(vel, dv_vec) + current_orbit = cartesian_to_orbital_elements(pos, vel, parent_mass) + current_orbit = propagate(current_orbit, dt_burn_step, parent_mass) + + return current_orbit + + +def main(): + dt = 60.0 + earth_mass = 5.972e24 + mu = G * earth_mass + earth = None # filled below + + # Setup: load config and get Hohmann_Transfer craft + sim = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft = sim.spacecraft[0] # Hohmann_Transfer + earth = sim.bodies[1] + + # Initialize craft state from orbital elements + pos, vel = orbital_to_cartesian(craft.orbit, earth.mass) + craft.local_pos = pos + craft.local_vel = vel + + a0 = craft.orbit.a + e0 = craft.orbit.e + r0 = vmag(craft.local_pos) + v0 = vmag(craft.local_vel) + + print("// === Config loading ===") + print(f"// body_count = {len(sim.bodies)}") + print(f"// craft_count = {len(sim.spacecraft)}") + print(f"// maneuver_count = {len(sim.maneuvers)}") + print(f"// craft[0] = \"{craft.name}\", parent_index = {craft.parent_index}") + print() + + # Test: Hohmann transfer - first burn at perigee + sim_h1 = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_h1 = sim_h1.spacecraft[0] + earth_h1 = sim_h1.bodies[1] + pos_h1, vel_h1 = orbital_to_cartesian(craft_h1.orbit, earth_h1.mass) + craft_h1.local_pos = pos_h1 + craft_h1.local_vel = vel_h1 + + v_before = vmag(craft_h1.local_vel) + apply_impulsive_burn(craft_h1, BurnDirection.PROGRADE, 2440.0, earth_h1.mass) + v_after = vmag(craft_h1.local_vel) + r_after = vmag(craft_h1.local_pos) + + post_burn_els = cartesian_to_orbital_elements(craft_h1.local_pos, craft_h1.local_vel, earth_h1.mass) + a_after = post_burn_els.a + e_after = post_burn_els.e + + print("// === Hohmann transfer: first burn (2440 m/s prograde) ===") + print(f"// v_before = {v_before:.6f}") + print(f"// v_after = {v_after:.6f}") + print(f"// r_after = {r_after:.6f}") + print(f"// a_after = {a_after:.6f}") + print(f"// e_after = {e_after:.15f}") + print() + + # Test: Hohmann transfer - second burn at apogee + sim_h2 = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_h2 = sim_h2.spacecraft[0] + earth_h2 = sim_h2.bodies[1] + pos_h2, vel_h2 = orbital_to_cartesian(craft_h2.orbit, earth_h2.mass) + craft_h2.local_pos = pos_h2 + craft_h2.local_vel = vel_h2 + + # First burn + apply_impulsive_burn(craft_h2, BurnDirection.PROGRADE, 2440.0, earth_h2.mass) + els_after_1 = cartesian_to_orbital_elements(craft_h2.local_pos, craft_h2.local_vel, earth_h2.mass) + + # Propagate to apogee (true anomaly = pi) + els_apogee = els_after_1 + els_apogee.nu = math.pi + pos_apogee, vel_apogee = orbital_to_cartesian(els_apogee, earth_h2.mass) + craft_h2.local_pos = pos_apogee + craft_h2.local_vel = vel_apogee + + # Second burn + apply_impulsive_burn(craft_h2, BurnDirection.PROGRADE, 1500.0, earth_h2.mass) + final_els = cartesian_to_orbital_elements(craft_h2.local_pos, craft_h2.local_vel, earth_h2.mass) + a_final = final_els.a + e_final = final_els.e + + print("// === Hohmann transfer: second burn at apogee (1500 m/s prograde) ===") + print(f"// a_after_first = {els_after_1.a:.6f}") + print(f"// e_after_first = {els_after_1.e:.15f}") + print(f"// a_final = {a_final:.6f}") + print(f"// e_final = {e_final:.15f}") + print() + + # Test: Large burn -> hyperbolic orbit + sim_large = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_large = sim_large.spacecraft[5] # Large_Delta_v + earth_large = sim_large.bodies[1] + pos_l, vel_l = orbital_to_cartesian(craft_large.orbit, earth_large.mass) + craft_large.local_pos = pos_l + craft_large.local_vel = vel_l + + v_esc = math.sqrt(2.0 * G * earth_large.mass / vmag(craft_large.local_pos)) + v_before_l = vmag(craft_large.local_vel) + apply_impulsive_burn(craft_large, BurnDirection.PROGRADE, 12000.0, earth_large.mass) + v_after_l = vmag(craft_large.local_vel) + + hyper_els = cartesian_to_orbital_elements(craft_large.local_pos, craft_large.local_vel, earth_large.mass) + e_hyper = hyper_els.e + a_hyper = hyper_els.a + + # Vis-viva check + r_hyper = vmag(craft_large.local_pos) + vis_viva_expected = v_after_l ** 2 + vis_viva_calc = G * earth_large.mass * (2.0 / r_hyper - 1.0 / a_hyper) + vis_viva_err = abs(vis_viva_expected - vis_viva_calc) / vis_viva_expected + + print("// === Large burn (12000 m/s prograde) -> hyperbolic ===") + print(f"// v_before = {v_before_l:.6f}") + print(f"// v_escape = {v_esc:.6f}") + print(f"// v_after = {v_after_l:.6f}") + print(f"// e = {e_hyper:.15f}") + print(f"// a = {a_hyper:.6f}") + print(f"// vis_viva_error = {vis_viva_err:.15e}") + print() + + # Test: Energy conservation - prograde burn + sim_e1 = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_e1 = sim_e1.spacecraft[0] + earth_e1 = sim_e1.bodies[1] + pos_e1, vel_e1 = orbital_to_cartesian(craft_e1.orbit, earth_e1.mass) + craft_e1.local_pos = pos_e1 + craft_e1.local_vel = vel_e1 + + m_craft = craft_e1.mass + v_init = craft_e1.local_vel + ke_init = 0.5 * m_craft * vdot(v_init, v_init) + r_init = vmag(craft_e1.local_pos) + pe_init = -G * m_craft * earth_e1.mass / r_init + E_init = ke_init + pe_init + + v_before_e = vmag(v_init) + apply_impulsive_burn(craft_e1, BurnDirection.PROGRADE, 2440.0, earth_e1.mass) + v_final_e = craft_e1.local_vel + ke_final = 0.5 * m_craft * vdot(v_final_e, v_final_e) + pe_final = -G * m_craft * earth_e1.mass / vmag(craft_e1.local_pos) + E_final = ke_final + pe_final + + dE_actual = E_final - E_init + dv_vec = vsub(v_final_e, v_init) + dE_expected = vdot(v_init, dv_vec) * m_craft + 0.5 * m_craft * vdot(dv_vec, dv_vec) + dE_err = abs(dE_actual - dE_expected) / abs(dE_expected) + + print("// === Energy: prograde burn (2440 m/s) ===") + print(f"// E_init = {E_init:.6f}") + print(f"// E_final = {E_final:.6f}") + print(f"// dE_actual = {dE_actual:.6f}") + print(f"// dE_expected = {dE_expected:.6f}") + print(f"// dE_relative_error = {dE_err:.15e}") + print() + + # Test: Energy conservation - retrograde burn + sim_e2 = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_e2 = sim_e2.spacecraft[0] + earth_e2 = sim_e2.bodies[1] + pos_e2, vel_e2 = orbital_to_cartesian(craft_e2.orbit, earth_e2.mass) + craft_e2.local_pos = pos_e2 + craft_e2.local_vel = vel_e2 + + m_e2 = craft_e2.mass + v_init_e2 = craft_e2.local_vel + ke_init_e2 = 0.5 * m_e2 * vdot(v_init_e2, v_init_e2) + pe_init_e2 = -G * m_e2 * earth_e2.mass / vmag(craft_e2.local_pos) + E_init_e2 = ke_init_e2 + pe_init_e2 + + apply_custom_burn(craft_e2, vscale(vnorm(v_init_e2), -1000.0)) + v_final_e2 = craft_e2.local_vel + ke_final_e2 = 0.5 * m_e2 * vdot(v_final_e2, v_final_e2) + pe_final_e2 = -G * m_e2 * earth_e2.mass / vmag(craft_e2.local_pos) + E_final_e2 = ke_final_e2 + pe_final_e2 + + dE_actual_e2 = E_final_e2 - E_init_e2 + dv_vec_e2 = vsub(v_final_e2, v_init_e2) + dE_expected_e2 = vdot(v_init_e2, dv_vec_e2) * m_e2 + 0.5 * m_e2 * vdot(dv_vec_e2, dv_vec_e2) + dE_err_e2 = abs(dE_actual_e2 - dE_expected_e2) / abs(dE_expected_e2) + + print("// === Energy: retrograde burn (1000 m/s) ===") + print(f"// E_init = {E_init_e2:.6f}") + print(f"// E_final = {E_final_e2:.6f}") + print(f"// dE_actual = {dE_actual_e2:.6f}") + print(f"// dE_expected = {dE_expected_e2:.6f}") + print(f"// dE_relative_error = {dE_err_e2:.15e}") + print() + + # Test: Round-trip conversion stability + sim_rt = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_rt = sim_rt.spacecraft[0] + earth_rt = sim_rt.bodies[1] + pos_rt, vel_rt = orbital_to_cartesian(craft_rt.orbit, earth_rt.mass) + craft_rt.local_pos = pos_rt + craft_rt.local_vel = vel_rt + + orig_a = craft_rt.orbit.a + orig_e = craft_rt.orbit.e + + for _ in range(5): + els_rt = cartesian_to_orbital_elements(craft_rt.local_pos, craft_rt.local_vel, earth_rt.mass) + pos_rt, vel_rt = orbital_to_cartesian(els_rt, earth_rt.mass) + craft_rt.local_pos = pos_rt + craft_rt.local_vel = vel_rt + + final_els_rt = cartesian_to_orbital_elements(craft_rt.local_pos, craft_rt.local_vel, earth_rt.mass) + a_err_rt = abs(final_els_rt.a - orig_a) / orig_a + e_err_rt = abs(final_els_rt.e - orig_e) + + print("// === Round-trip conversion stability (5 iterations) ===") + print(f"// orig_a = {orig_a:.6f}") + print(f"// final_a = {final_els_rt.a:.6f}") + print(f"// a_relative_error = {a_err_rt:.15e}") + print(f"// orig_e = {orig_e:.15f}") + print(f"// final_e = {final_els_rt.e:.15f}") + print(f"// e_absolute_error = {e_err_rt:.15e}") + print() + + # Test: Burn direction orthogonality + sim_dir = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_dir = sim_dir.spacecraft[0] + earth_dir = sim_dir.bodies[1] + pos_dir, vel_dir = orbital_to_cartesian(craft_dir.orbit, earth_dir.mass) + craft_dir.local_pos = pos_dir + craft_dir.local_vel = vel_dir + + pro = get_burn_direction(BurnDirection.PROGRADE, pos_dir, vel_dir) + retro = get_burn_direction(BurnDirection.RETROGRADE, pos_dir, vel_dir) + norm_dir = get_burn_direction(BurnDirection.NORMAL, pos_dir, vel_dir) + anti = get_burn_direction(BurnDirection.ANTINORMAL, pos_dir, vel_dir) + rad_in = get_burn_direction(BurnDirection.RADIAL_IN, pos_dir, vel_dir) + rad_out = get_burn_direction(BurnDirection.RADIAL_OUT, pos_dir, vel_dir) + + dot_pro_retro = vdot(pro, retro) + dot_norm_anti = vdot(norm_dir, anti) + dot_rad_in_out = vdot(rad_in, rad_out) + + print("// === Burn direction orthogonality ===") + print(f"// prograde . retrograde = {dot_pro_retro:.15f}") + print(f"// normal . antinormal = {dot_norm_anti:.15f}") + print(f"// radial_in . radial_out = {dot_rad_in_out:.15f}") + print() + + # Test: Continuous burn (100 steps, 100 m/s total) + sim_cb = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_cb = sim_cb.spacecraft[6] # Low_Thrust_Ion + earth_cb = sim_cb.bodies[1] + + initial_a_cb = craft_cb.orbit.a + initial_e_cb = craft_cb.orbit.e + + final_cb = simulate_continuous_burn(craft_cb.orbit, earth_cb.mass, + 100.0, 5000.0, 100, BurnDirection.PROGRADE) + + a_cb = final_cb.a + e_cb = final_cb.e + v_circ_init = math.sqrt(mu / initial_a_cb) + v_circ_final = math.sqrt(mu / a_cb) + eps_init = -mu / (2.0 * initial_a_cb) + eps_final = -mu / (2.0 * a_cb) + delta_eps = eps_final - eps_init + expected_dv_from_energy = delta_eps / v_circ_init + rel_err_cb = abs(expected_dv_from_energy - 100.0) / 100.0 + + print("// === Continuous burn: 100 steps, 100 m/s total prograde ===") + print(f"// initial_a = {initial_a_cb:.6f}") + print(f"// final_a = {a_cb:.6f}") + print(f"// initial_e = {initial_e_cb:.15f}") + print(f"// final_e = {e_cb:.15f}") + print(f"// v_circ_initial = {v_circ_init:.6f}") + print(f"// v_circ_final = {v_circ_final:.6f}") + print(f"// delta_specific_energy = {delta_eps:.6f}") + print(f"// expected_dv_from_energy = {expected_dv_from_energy:.6f}") + print(f"// relative_error = {rel_err_cb:.15e}") + print() + + # Test: Multi-burn continuous sequence + sim_mb = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_mb = sim_mb.spacecraft[7] # Multi_Burn_Sequence + earth_mb = sim_mb.bodies[1] + + initial_a_mb = craft_mb.orbit.a + + orbit_after_1 = simulate_continuous_burn(craft_mb.orbit, earth_mb.mass, + 50.0, 2000.0, 20, BurnDirection.PROGRADE) + final_mb = simulate_continuous_burn(orbit_after_1, earth_mb.mass, + 75.0, 3000.0, 30, BurnDirection.PROGRADE) + + a_mb = final_mb.a + v_circ_init_mb = math.sqrt(mu / initial_a_mb) + eps_init_mb = -mu / (2.0 * initial_a_mb) + eps_final_mb = -mu / (2.0 * a_mb) + delta_eps_mb = eps_final_mb - eps_init_mb + expected_dv_mb = delta_eps_mb / v_circ_init_mb + rel_err_mb = abs(expected_dv_mb - 125.0) / 125.0 + + print("// === Multi-burn continuous: 50+75 m/s total prograde ===") + print(f"// initial_a = {initial_a_mb:.6f}") + print(f"// final_a = {a_mb:.6f}") + print(f"// total_dv = 125.0") + print(f"// expected_dv_from_energy = {expected_dv_mb:.6f}") + print(f"// relative_error = {rel_err_mb:.15e}") + print() + + # Test: Mode transition (elliptical orbit, continuous burn) + sim_mt = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_mt = sim_mt.spacecraft[8] # Mode_Transition + earth_mt = sim_mt.bodies[1] + + initial_a_mt = craft_mt.orbit.a + initial_e_mt = craft_mt.orbit.e + + final_mt = simulate_continuous_burn(craft_mt.orbit, earth_mt.mass, + 200.0, 4000.0, 80, BurnDirection.PROGRADE) + + a_mt = final_mt.a + e_mt = final_mt.e + mu_mt = G * earth_mt.mass + energy_before = -mu_mt / (2.0 * initial_a_mt) + energy_after = -mu_mt / (2.0 * a_mt) + energy_change = energy_after - energy_before + + v_init_mt = math.sqrt(mu_mt / initial_a_mt) + v_final_mt = math.sqrt(mu_mt / a_mt) + expected_energy_change = 0.5 * (v_init_mt + v_final_mt) * 200.0 + + print("// === Mode transition: 80 steps, 200 m/s total prograde (e=0.3) ===") + print(f"// initial_a = {initial_a_mt:.6f}") + print(f"// initial_e = {initial_e_mt:.15f}") + print(f"// final_a = {a_mt:.6f}") + print(f"// final_e = {e_mt:.15f}") + print(f"// energy_change = {energy_change:.6f}") + print(f"// expected_energy_change = {expected_energy_change:.6f}") + print() + + # Test: Continuous energy conservation + sim_ec = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_ec = sim_ec.spacecraft[9] # Energy_Conservation + earth_ec = sim_ec.bodies[1] + + ke_init_ec = 0.5 * craft_ec.mass * vdot(craft_ec.local_vel, craft_ec.local_vel) + pe_init_ec = -G * craft_ec.mass * earth_ec.mass / vmag(craft_ec.local_pos) + E_init_ec = ke_init_ec + pe_init_ec + + final_ec = simulate_continuous_burn(craft_ec.orbit, earth_ec.mass, + 150.0, 6000.0, 120, BurnDirection.PROGRADE) + + pos_ec, vel_ec = orbital_to_cartesian(final_ec, earth_ec.mass) + temp_craft = Spacecraft(name="temp", mass=craft_ec.mass, parent_index=craft_ec.parent_index, + orbit=final_ec, local_pos=pos_ec, local_vel=vel_ec, + global_pos=(0, 0, 0), global_vel=(0, 0, 0)) + ke_final_ec = 0.5 * craft_ec.mass * vdot(vel_ec, vel_ec) + pe_final_ec = -G * craft_ec.mass * earth_ec.mass / vmag(pos_ec) + E_final_ec = ke_final_ec + pe_final_ec + + total_dE_ec = E_final_ec - E_init_ec + expected_dE_approx = craft_ec.mass * math.sqrt(mu / craft_ec.orbit.a) * 150.0 + rel_err_ec = abs(total_dE_ec - expected_dE_approx) / expected_dE_approx + + print("// === Continuous energy conservation: 120 steps, 150 m/s ===") + print(f"// E_init = {E_init_ec:.6f}") + print(f"// E_final = {E_final_ec:.6f}") + print(f"// total_dE = {total_dE_ec:.6f}") + print(f"// expected_approx = {expected_dE_approx:.6f}") + print(f"// relative_error = {rel_err_ec:.15e}") + print() + + # Test: Continuous vs impulsive comparison + sim_cv = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_cv = sim_cv.spacecraft[6] # Low_Thrust_Ion + earth_cv = sim_cv.bodies[1] + + orbit_cont = simulate_continuous_burn(craft_cv.orbit, earth_cv.mass, + 100.0, 5000.0, 100, BurnDirection.PROGRADE) + orbit_imp = simulate_continuous_burn(craft_cv.orbit, earth_cv.mass, + 100.0, 5000.0, 1, BurnDirection.PROGRADE) + + diff_a = abs(orbit_cont.a - orbit_imp.a) + rel_diff_a = diff_a / orbit_cont.a * 100.0 + + v_cont = math.sqrt(mu / orbit_cont.a) + v_imp = math.sqrt(mu / orbit_imp.a) + v_diff = abs(v_cont - v_imp) + + print("// === Continuous vs impulsive (100 steps vs 1 step) ===") + print(f"// continuous_a = {orbit_cont.a:.6f}") + print(f"// impulsive_a = {orbit_imp.a:.6f}") + print(f"// a_difference = {diff_a:.6f}") + print(f"// a_relative_diff_pct = {rel_diff_a:.6f}%") + print(f"// v_continuous = {v_cont:.6f}") + print(f"// v_impulsive = {v_imp:.6f}") + print(f"// v_difference = {v_diff:.6f}") + print() + + # Test: Propagation during burn - path length + sim_prop = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_prop = sim_prop.spacecraft[6] # Low_Thrust_Ion + earth_prop = sim_prop.bodies[1] + + current_orbit = craft_prop.orbit + dt_burn = 5000.0 / 100 + dv_per = 100.0 / 100 + + positions = [] + for i in range(101): + pos, vel = orbital_to_cartesian(current_orbit, earth_prop.mass) + positions.append(pos) + if i < 100: + burn_dir = get_burn_direction(BurnDirection.PROGRADE, pos, vel) + vel = vadd(vel, vscale(burn_dir, dv_per)) + current_orbit = cartesian_to_orbital_elements(pos, vel, earth_prop.mass) + current_orbit = propagate(current_orbit, dt_burn, earth_prop.mass) + + total_path = sum(vmag(vsub(positions[i], positions[i-1])) for i in range(1, len(positions))) + straight = vmag(vsub(positions[100], positions[0])) + r_start = vmag(positions[0]) + r_end = vmag(positions[100]) + + # Expected semi-major axis from energy + v_init_prop = math.sqrt(mu / craft_prop.orbit.a) + eps_init_prop = -mu / (2.0 * craft_prop.orbit.a) + eps_final_prop = eps_init_prop + v_init_prop * 100.0 + a_expected_prop = -mu / (2.0 * eps_final_prop) + e_init_prop = craft_prop.orbit.e + r_peri = a_expected_prop * (1.0 - e_init_prop) + r_apo = a_expected_prop * (1.0 + e_init_prop) + + print("// === Propagation during burn: path length ===") + print(f"// total_path_length = {total_path:.6f}") + print(f"// straight_line = {straight:.6f}") + print(f"// r_start = {r_start:.6f}") + print(f"// r_end = {r_end:.6f}") + print(f"// a_expected = {a_expected_prop:.6f}") + print(f"// r_peri_expected = {r_peri:.6f}") + print(f"// r_apo_expected = {r_apo:.6f}") + print() + + # Test: Numerical stability - monotonicity + sim_stab = Simulator("tests/test_hybrid_burns.toml", dt=dt) + craft_stab = sim_stab.spacecraft[6] + earth_stab = sim_stab.bodies[1] + + current_orbit = craft_stab.orbit + dt_burn_s = 5000.0 / 100 + dv_per_s = 100.0 / 100 + + a_history = [] + e_history = [] + for i in range(100): + pos, vel = orbital_to_cartesian(current_orbit, earth_stab.mass) + burn_dir = get_burn_direction(BurnDirection.PROGRADE, pos, vel) + vel = vadd(vel, vscale(burn_dir, dv_per_s)) + current_orbit = cartesian_to_orbital_elements(pos, vel, earth_stab.mass) + current_orbit = propagate(current_orbit, dt_burn_s, earth_stab.mass) + a_history.append(current_orbit.a) + e_history.append(current_orbit.e) + + monotonic = all(a_history[i] >= a_history[i-1] for i in range(1, len(a_history))) + max_e = max(e_history) + min_e = min(e_history) + + initial_a_s = craft_stab.orbit.a + final_a_s = a_history[-1] + total_change_s = final_a_s - initial_a_s + avg_change_s = total_change_s / 100 + + max_dev_s = max(abs(a_history[i] - (initial_a_s + (i+1) * avg_change_s)) for i in range(100)) + + print("// === Numerical stability: monotonicity ===") + print(f"// monotonic_increase = {monotonic}") + print(f"// max_eccentricity = {max_e:.15f}") + print(f"// min_eccentricity = {min_e:.15f}") + print(f"// total_a_change = {total_change_s:.6f}") + print(f"// max_deviation_from_linear = {max_dev_s:.6f}") + print(f"// max_deviation_pct = {max_dev_s / total_change_s * 100:.6f}%") + print() + + +if __name__ == "__main__": + main() diff --git a/tests/test_hybrid_burns.cpp b/tests/test_hybrid_burns.cpp new file mode 100644 index 0000000..06c81b8 --- /dev/null +++ b/tests/test_hybrid_burns.cpp @@ -0,0 +1,643 @@ +#include +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/orbital_objects.h" +#include "../src/maneuver.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include +#include + +using Catch::Matchers::WithinAbs; + +// Simulate continuous/low-thrust burn with sub-steps (replaces numerical integrator) +static OrbitalElements simulate_continuous_burn(OrbitalElements initial_orbit, double parent_mass, + double total_dv, double burn_duration, + int num_steps, BurnDirection direction) { + OrbitalElements current_orbit = initial_orbit; + double dt_burn = burn_duration / num_steps; + double dv_per = total_dv / num_steps; + + for (int i = 0; i < num_steps; i++) { + Vec3 pos, vel; + orbital_elements_to_cartesian(current_orbit, parent_mass, &pos, &vel); + + Vec3 dir = get_burn_direction_vector(direction, pos, vel); + Vec3 dv_vec = vec3_scale(dir, dv_per); + vel = vec3_add(vel, dv_vec); + + current_orbit = cartesian_to_orbital_elements(pos, vel, parent_mass); + current_orbit = propagate_orbital_elements(current_orbit, dt_burn, parent_mass); + } + return current_orbit; +} + +SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") { + const double TIME_STEP = 60.0; + const double MU_EARTH = G * 5.972e24; + + SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_hybrid_burns.toml")); + + // Helper: initialize a spacecraft from its orbital elements + auto init_craft = [&](Spacecraft* craft, CelestialBody* parent) { + Vec3 pos, vel; + orbital_elements_to_cartesian(craft->orbit, parent->mass, &pos, &vel); + craft->local_position = pos; + craft->local_velocity = vel; + }; + + // Helper: find maneuver by name + auto find_maneuver = [&](const char* name) -> int { + for (int i = 0; i < sim->maneuver_count; i++) { + if (strcmp(sim->maneuvers[i].name, name) == 0) return i; + } + return -1; + }; + + // Helper: execute maneuver by name (sets time, calls execute_maneuver) + auto exec_by_name = [&](const char* name, Spacecraft* craft) { + int idx = find_maneuver(name); + REQUIRE(idx >= 0); + Maneuver* m = &sim->maneuvers[idx]; + REQUIRE(!m->executed); + if (m->trigger_type == TRIGGER_TIME) { + sim->time = m->trigger_value; + } + execute_maneuver(m, craft, sim, sim->time); + REQUIRE(m->executed); + REQUIRE_THAT(m->executed_time, WithinAbs(sim->time, 0.001)); + }; + + // Shared fixtures + CelestialBody* sun = &sim->bodies[0]; + CelestialBody* earth = &sim->bodies[1]; + Spacecraft* hohmann = &sim->spacecraft[0]; + Spacecraft* large_dv = &sim->spacecraft[5]; + Spacecraft* low_thrust = &sim->spacecraft[6]; + Spacecraft* multi_burn = &sim->spacecraft[7]; + Spacecraft* mode_trans = &sim->spacecraft[8]; + Spacecraft* energy_cons = &sim->spacecraft[9]; + + SECTION("config loads correctly: 2 bodies, 10 spacecraft, 7 maneuvers") { + REQUIRE(sim->body_count == 2); + REQUIRE(std::string(sun->name) == "Sun"); + REQUIRE(std::string(earth->name) == "Earth"); + REQUIRE(sim->craft_count == 10); + REQUIRE(sim->maneuver_count == 7); + + REQUIRE(std::string(hohmann->name) == "Hohmann_Transfer"); + REQUIRE(hohmann->parent_index == 1); + REQUIRE(std::string(large_dv->name) == "Large_Delta_v"); + } + + SECTION("first burn at perigee raises apogee") { + init_craft(hohmann, earth); + + const double v_before = vec3_magnitude(hohmann->local_velocity); + exec_by_name("hohmann_burn_1", hohmann); + const double v_after = vec3_magnitude(hohmann->local_velocity); + + REQUIRE(v_after > v_before); + + const auto post_els = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + INFO("v_before: " << v_before << " m/s"); + INFO("v_after: " << v_after << " m/s"); + INFO("a_after: " << post_els.semi_major_axis << " m"); + INFO("e_after: " << post_els.eccentricity); + + REQUIRE_THAT(post_els.semi_major_axis, WithinAbs(25762376.160113, 1.0)); + REQUIRE_THAT(post_els.eccentricity, WithinAbs(0.737174864697325, 1e-6)); + } + + SECTION("second burn at apogee circularizes orbit") { + init_craft(hohmann, earth); + + exec_by_name("hohmann_burn_1", hohmann); + const auto after_1 = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + // Propagate to apogee + auto apogee_els = after_1; + apogee_els.true_anomaly = M_PI; + Vec3 apogee_pos, apogee_vel; + orbital_elements_to_cartesian(apogee_els, earth->mass, &apogee_pos, &apogee_vel); + hohmann->local_position = apogee_pos; + hohmann->local_velocity = apogee_vel; + + exec_by_name("hohmann_burn_2", hohmann); + const auto final_els = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + INFO("a_after_first: " << after_1.semi_major_axis); + INFO("a_final: " << final_els.semi_major_axis); + INFO("e_final: " << final_els.eccentricity); + + REQUIRE(final_els.semi_major_axis > after_1.semi_major_axis); + REQUIRE(final_els.eccentricity < after_1.eccentricity); + REQUIRE(final_els.eccentricity < 0.1); + } + + SECTION("large prograde burn produces hyperbolic orbit") { + init_craft(large_dv, earth); + + const double v_before = vec3_magnitude(large_dv->local_velocity); + const double r = vec3_magnitude(large_dv->local_position); + const double v_escape = sqrt(2.0 * G * earth->mass / r); + + exec_by_name("large_burn", large_dv); + const double v_after = vec3_magnitude(large_dv->local_velocity); + + INFO("v_before: " << v_before << " m/s"); + INFO("v_escape: " << v_escape << " m/s"); + INFO("v_after: " << v_after << " m/s"); + + REQUIRE(v_after > v_escape); + + const auto hyper_els = cartesian_to_orbital_elements( + large_dv->local_position, large_dv->local_velocity, earth->mass); + + INFO("e: " << hyper_els.eccentricity); + INFO("a: " << hyper_els.semi_major_axis); + + REQUIRE_THAT(hyper_els.eccentricity, WithinAbs(5.709434906871548, 1e-6)); + REQUIRE_THAT(hyper_els.semi_major_axis, WithinAbs(-1486377.906994, 1.0)); + } + + SECTION("large burn satisfies vis-viva equation") { + init_craft(large_dv, earth); + + exec_by_name("large_burn", large_dv); + const auto hyper_els = cartesian_to_orbital_elements( + large_dv->local_position, large_dv->local_velocity, earth->mass); + + const double v_sq = vec3_magnitude(large_dv->local_velocity) + * vec3_magnitude(large_dv->local_velocity); + const double r = vec3_magnitude(large_dv->local_position); + const double vis_viva_calc = G * earth->mass * (2.0 / r - 1.0 / hyper_els.semi_major_axis); + + INFO("vis_viva_expected: " << v_sq); + INFO("vis_viva_calculated: " << vis_viva_calc); + + const double err = fabs(v_sq - vis_viva_calc) / v_sq; + REQUIRE_THAT(err, WithinAbs(0.0, 1e-12)); + } + + SECTION("prograde burn increases total energy") { + init_craft(hohmann, earth); + + const double m = hohmann->mass; + const Vec3 v_init = hohmann->local_velocity; + const double ke_init = 0.5 * m * vec3_dot(v_init, v_init); + const double r_init = vec3_magnitude(hohmann->local_position); + const double pe_init = -G * m * earth->mass / r_init; + const double E_init = ke_init + pe_init; + + exec_by_name("hohmann_burn_1", hohmann); + + const Vec3 v_final = hohmann->local_velocity; + const Vec3 dv = vec3_sub(v_final, v_init); + const double ke_final = 0.5 * m * vec3_dot(v_final, v_final); + const double pe_final = -G * m * earth->mass / vec3_magnitude(hohmann->local_position); + const double E_final = ke_final + pe_final; + + const double dE_actual = E_final - E_init; + const double dE_expected = vec3_dot(v_init, dv) * m + 0.5 * m * vec3_dot(dv, dv); + + INFO("E_init: " << E_init); + INFO("E_final: " << E_final); + INFO("dE_actual: " << dE_actual); + INFO("dE_expected: " << dE_expected); + + REQUIRE(E_final > E_init); + + const double dE_err = fabs(dE_actual - dE_expected) / fabs(dE_expected); + REQUIRE_THAT(dE_err, WithinAbs(0.0, 1e-12)); + } + + SECTION("retrograde burn decreases total energy") { + init_craft(hohmann, earth); + + const double m = hohmann->mass; + const Vec3 v_init = hohmann->local_velocity; + const double ke_init = 0.5 * m * vec3_dot(v_init, v_init); + const double pe_init = -G * m * earth->mass / vec3_magnitude(hohmann->local_position); + const double E_init = ke_init + pe_init; + + const Vec3 retro_dir = calculate_retrograde_dir(v_init); + apply_custom_burn(hohmann, vec3_scale(retro_dir, 1000.0)); + + const Vec3 v_final = hohmann->local_velocity; + const Vec3 dv = vec3_sub(v_final, v_init); + const double ke_final = 0.5 * m * vec3_dot(v_final, v_final); + const double pe_final = -G * m * earth->mass / vec3_magnitude(hohmann->local_position); + const double E_final = ke_final + pe_final; + + const double dE_actual = E_final - E_init; + const double dE_expected = vec3_dot(v_init, dv) * m + 0.5 * m * vec3_dot(dv, dv); + + INFO("E_init: " << E_init); + INFO("E_final: " << E_final); + INFO("dE_actual: " << dE_actual); + INFO("dE_expected: " << dE_expected); + + REQUIRE(E_final < E_init); + + const double dE_err = fabs(dE_actual - dE_expected) / fabs(dE_expected); + REQUIRE_THAT(dE_err, WithinAbs(0.0, 1e-12)); + } + + SECTION("orbital elements -> Cartesian -> burn -> orbital elements") { + init_craft(hohmann, earth); + + const auto orig_els = hohmann->orbit; + const auto recovered = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + INFO("orig_a: " << orig_els.semi_major_axis); + INFO("recovered_a: " << recovered.semi_major_axis); + INFO("orig_e: " << orig_els.eccentricity); + INFO("recovered_e: " << recovered.eccentricity); + + REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(orig_els.semi_major_axis, A_TOL)); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(orig_els.eccentricity, E_TOL)); + + exec_by_name("hohmann_burn_1", hohmann); + const auto post_burn = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + INFO("post_burn_a: " << post_burn.semi_major_axis); + INFO("post_burn_e: " << post_burn.eccentricity); + + REQUIRE(post_burn.semi_major_axis != recovered.semi_major_axis); + REQUIRE(post_burn.eccentricity != recovered.eccentricity); + } + + SECTION("multiple round-trip conversions maintain stability") { + init_craft(hohmann, earth); + + const auto orig_els = hohmann->orbit; + Vec3 pos = hohmann->local_position; + Vec3 vel = hohmann->local_velocity; + + for (int i = 0; i < 5; i++) { + const auto els = cartesian_to_orbital_elements(pos, vel, earth->mass); + orbital_elements_to_cartesian(els, earth->mass, &pos, &vel); + } + + const auto final_els = cartesian_to_orbital_elements(pos, vel, earth->mass); + const double a_err = fabs(final_els.semi_major_axis - orig_els.semi_major_axis) + / orig_els.semi_major_axis; + const double e_err = fabs(final_els.eccentricity - orig_els.eccentricity); + + INFO("orig_a: " << orig_els.semi_major_axis); + INFO("final_a: " << final_els.semi_major_axis); + INFO("orig_e: " << orig_els.eccentricity); + INFO("final_e: " << final_els.eccentricity); + + REQUIRE_THAT(a_err, WithinAbs(0.0, 1e-12)); + REQUIRE_THAT(e_err, WithinAbs(0.0, 1e-12)); + } + + SECTION("two-burn sequence raises orbit") { + init_craft(hohmann, earth); + + const auto init_els = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + exec_by_name("hohmann_burn_1", hohmann); + const auto after_1 = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + REQUIRE(after_1.semi_major_axis > init_els.semi_major_axis); + + // Propagate to apogee + auto apogee_els = after_1; + apogee_els.true_anomaly = M_PI; + Vec3 apogee_pos, apogee_vel; + orbital_elements_to_cartesian(apogee_els, earth->mass, &apogee_pos, &apogee_vel); + hohmann->local_position = apogee_pos; + hohmann->local_velocity = apogee_vel; + + exec_by_name("hohmann_burn_2", hohmann); + const auto after_2 = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + INFO("a_after_2: " << after_2.semi_major_axis); + INFO("e_after_2: " << after_2.eccentricity); + + REQUIRE(after_2.semi_major_axis > after_1.semi_major_axis); + REQUIRE(after_2.eccentricity < after_1.eccentricity); + } + + SECTION("three-burn sequence with plane change") { + init_craft(hohmann, earth); + + const auto init_els = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + // Burn 1: prograde 500 m/s + apply_custom_burn(hohmann, vec3_scale(calculate_prograde_dir(hohmann->local_velocity), 500.0)); + + // Burn 2: normal 300 m/s + apply_custom_burn(hohmann, vec3_scale(calculate_normal_dir(hohmann->local_position, + hohmann->local_velocity), 300.0)); + + // Burn 3: prograde 200 m/s + apply_custom_burn(hohmann, vec3_scale(calculate_prograde_dir(hohmann->local_velocity), 200.0)); + const auto after_3 = cartesian_to_orbital_elements( + hohmann->local_position, hohmann->local_velocity, earth->mass); + + INFO("init_a: " << init_els.semi_major_axis); + INFO("final_a: " << after_3.semi_major_axis); + INFO("init_inc: " << init_els.inclination); + INFO("final_inc: " << after_3.inclination); + + REQUIRE(after_3.semi_major_axis > init_els.semi_major_axis); + REQUIRE(after_3.inclination > init_els.inclination); + } + + SECTION("prograde and retrograde are opposite") { + init_craft(hohmann, earth); + + const Vec3 pro = calculate_prograde_dir(hohmann->local_velocity); + const Vec3 retro = calculate_retrograde_dir(hohmann->local_velocity); + const double dot = vec3_dot(pro, retro); + + INFO("prograde . retrograde: " << dot); + + REQUIRE_THAT(dot, WithinAbs(-1.0, V_TOL)); + } + + SECTION("normal and antinormal are opposite") { + init_craft(hohmann, earth); + + const Vec3 norm = calculate_normal_dir(hohmann->local_position, hohmann->local_velocity); + const Vec3 anti = calculate_antinormal_dir(hohmann->local_position, hohmann->local_velocity); + const double dot = vec3_dot(norm, anti); + + INFO("normal . antinormal: " << dot); + + REQUIRE_THAT(dot, WithinAbs(-1.0, V_TOL)); + } + + SECTION("radial in and radial out are opposite") { + init_craft(hohmann, earth); + + const Vec3 rad_in = calculate_radial_in_dir(hohmann->local_position); + const Vec3 rad_out = calculate_radial_out_dir(hohmann->local_position); + const double dot = vec3_dot(rad_in, rad_out); + + INFO("radial_in . radial_out: " << dot); + + REQUIRE_THAT(dot, WithinAbs(-1.0, V_TOL)); + } + + SECTION("continuous prograde burn raises semi-major axis") { + const auto final_els = simulate_continuous_burn( + low_thrust->orbit, earth->mass, 100.0, 5000.0, 100, BURN_PROGRADE); + + INFO("initial_a: " << low_thrust->orbit.semi_major_axis); + INFO("final_a: " << final_els.semi_major_axis); + INFO("final_e: " << final_els.eccentricity); + + REQUIRE(final_els.semi_major_axis > low_thrust->orbit.semi_major_axis); + + const double v_circ_init = sqrt(MU_EARTH / low_thrust->orbit.semi_major_axis); + // v_circ_final not needed for assertions + const double eps_init = -MU_EARTH / (2.0 * low_thrust->orbit.semi_major_axis); + const double eps_final = -MU_EARTH / (2.0 * final_els.semi_major_axis); + const double expected_dv = (eps_final - eps_init) / v_circ_init; + const double rel_err = fabs(expected_dv - 100.0) / 100.0; + + INFO("v_circ_init: " << v_circ_init); + INFO("expected_dv: " << expected_dv); + INFO("relative_error: " << rel_err); + + REQUIRE_THAT(rel_err, WithinAbs(0.0, 0.01)); + REQUIRE(final_els.eccentricity < 0.01); + } + + SECTION("continuous multi-burn sequence raises orbit") { + const auto after_1 = simulate_continuous_burn( + multi_burn->orbit, earth->mass, 50.0, 2000.0, 20, BURN_PROGRADE); + + REQUIRE(after_1.semi_major_axis > multi_burn->orbit.semi_major_axis); + + const auto final_els = simulate_continuous_burn( + after_1, earth->mass, 75.0, 3000.0, 30, BURN_PROGRADE); + + INFO("final_a: " << final_els.semi_major_axis); + + REQUIRE(final_els.semi_major_axis > after_1.semi_major_axis); + + const double v_circ_init = sqrt(MU_EARTH / multi_burn->orbit.semi_major_axis); + const double eps_init = -MU_EARTH / (2.0 * multi_burn->orbit.semi_major_axis); + const double eps_final = -MU_EARTH / (2.0 * final_els.semi_major_axis); + const double expected_dv = (eps_final - eps_init) / v_circ_init; + const double rel_err = fabs(expected_dv - 125.0) / 125.0; + + INFO("expected_dv: " << expected_dv); + INFO("relative_error: " << rel_err); + + REQUIRE_THAT(rel_err, WithinAbs(0.0, 0.01)); + } + + SECTION("continuous burn on elliptical orbit raises semi-major axis") { + const auto final_els = simulate_continuous_burn( + mode_trans->orbit, earth->mass, 200.0, 4000.0, 80, BURN_PROGRADE); + + INFO("initial_a: " << mode_trans->orbit.semi_major_axis); + INFO("final_a: " << final_els.semi_major_axis); + INFO("initial_e: " << mode_trans->orbit.eccentricity); + INFO("final_e: " << final_els.eccentricity); + + REQUIRE(final_els.semi_major_axis > mode_trans->orbit.semi_major_axis); + + const double energy_before = -MU_EARTH / (2.0 * mode_trans->orbit.semi_major_axis); + const double energy_after = -MU_EARTH / (2.0 * final_els.semi_major_axis); + const double energy_change = energy_after - energy_before; + + INFO("energy_change: " << energy_change); + + REQUIRE(fabs(energy_change) > 0.0); + } + + SECTION("continuous burn energy increases monotonically") { + const auto final_els = simulate_continuous_burn( + energy_cons->orbit, earth->mass, 150.0, 6000.0, 120, BURN_PROGRADE); + + const double m = energy_cons->mass; + const double ke_init = 0.5 * m * vec3_dot(energy_cons->local_velocity, + energy_cons->local_velocity); + const double pe_init = -G * m * earth->mass / vec3_magnitude(energy_cons->local_position); + const double E_init = ke_init + pe_init; + + Vec3 pos, vel; + orbital_elements_to_cartesian(final_els, earth->mass, &pos, &vel); + const double ke_final = 0.5 * m * vec3_dot(vel, vel); + const double pe_final = -G * m * earth->mass / vec3_magnitude(pos); + const double E_final = ke_final + pe_final; + + const double total_dE = E_final - E_init; + const double v_circ = sqrt(MU_EARTH / energy_cons->orbit.semi_major_axis); + const double expected_approx = m * v_circ * 150.0; + const double rel_err = fabs(total_dE - expected_approx) / expected_approx; + + INFO("E_init: " << E_init); + INFO("E_final: " << E_final); + INFO("total_dE: " << total_dE); + INFO("expected_approx: " << expected_approx); + INFO("relative_error: " << rel_err); + + REQUIRE(total_dE > 0.0); + REQUIRE_THAT(rel_err, WithinAbs(0.0, 0.01)); + } + + SECTION("continuous vs impulsive burn agree within 1%") { + const auto orbit_cont = simulate_continuous_burn( + low_thrust->orbit, earth->mass, 100.0, 5000.0, 100, BURN_PROGRADE); + const auto orbit_imp = simulate_continuous_burn( + low_thrust->orbit, earth->mass, 100.0, 5000.0, 1, BURN_PROGRADE); + + const double diff_a = fabs(orbit_cont.semi_major_axis - orbit_imp.semi_major_axis); + const double rel_diff = diff_a / orbit_cont.semi_major_axis * 100.0; + + const double v_cont = sqrt(MU_EARTH / orbit_cont.semi_major_axis); + const double v_imp = sqrt(MU_EARTH / orbit_imp.semi_major_axis); + const double v_diff = fabs(v_cont - v_imp); + + INFO("continuous_a: " << orbit_cont.semi_major_axis); + INFO("impulsive_a: " << orbit_imp.semi_major_axis); + INFO("rel_diff_a: " << rel_diff << "%"); + INFO("v_difference: " << v_diff); + + REQUIRE_THAT(rel_diff, WithinAbs(0.0, 1.0)); + REQUIRE(v_diff < 2.0); + } + + SECTION("propagation during burn: path length > straight line") { + OrbitalElements current = low_thrust->orbit; + const double dt_burn = 5000.0 / 100; + const double dv_per = 100.0 / 100; + + Vec3 pos_start, pos_end; + double total_path = 0.0; + Vec3 prev_pos; + bool first = true; + + for (int i = 0; i <= 100; i++) { + Vec3 pos, vel; + orbital_elements_to_cartesian(current, earth->mass, &pos, &vel); + + if (i == 0) pos_start = pos; + if (i == 100) pos_end = pos; + + if (!first) { + total_path += vec3_distance(prev_pos, pos); + } + first = false; + prev_pos = pos; + + if (i < 100) { + Vec3 dir = get_burn_direction_vector(BURN_PROGRADE, pos, vel); + vel = vec3_add(vel, vec3_scale(dir, dv_per)); + current = cartesian_to_orbital_elements(pos, vel, earth->mass); + current = propagate_orbital_elements(current, dt_burn, earth->mass); + } + } + + const double straight = vec3_distance(pos_start, pos_end); + const double r_start = vec3_magnitude(pos_start); + const double r_end = vec3_magnitude(pos_end); + + INFO("total_path: " << total_path); + INFO("straight_line: " << straight); + INFO("r_start: " << r_start); + INFO("r_end: " << r_end); + + REQUIRE(total_path > straight); + REQUIRE(r_end > r_start); + + // Check final radius within expected bounds + const double v_init = sqrt(MU_EARTH / low_thrust->orbit.semi_major_axis); + const double eps_init = -MU_EARTH / (2.0 * low_thrust->orbit.semi_major_axis); + const double eps_final = eps_init + v_init * 100.0; + const double a_expected = -MU_EARTH / (2.0 * eps_final); + const double e_init = low_thrust->orbit.eccentricity; + const double r_peri = a_expected * (1.0 - e_init); + const double r_apo = a_expected * (1.0 + e_init); + + INFO("a_expected: " << a_expected); + INFO("r_peri: " << r_peri); + INFO("r_apo: " << r_apo); + + REQUIRE(r_end >= r_peri - 1e5); + REQUIRE(r_end <= r_apo + 1e5); + } + + SECTION("continuous burn: semi-major axis increases monotonically") { + OrbitalElements current = low_thrust->orbit; + const double dt_burn = 5000.0 / 100; + const double dv_per = 100.0 / 100; + + bool monotonic = true; + double max_e = 0.0; + double min_e = 1.0; + double initial_a = current.semi_major_axis; + double a_prev = current.semi_major_axis; + + for (int i = 0; i < 100; i++) { + Vec3 pos, vel; + orbital_elements_to_cartesian(current, earth->mass, &pos, &vel); + + Vec3 dir = get_burn_direction_vector(BURN_PROGRADE, pos, vel); + vel = vec3_add(vel, vec3_scale(dir, dv_per)); + + current = cartesian_to_orbital_elements(pos, vel, earth->mass); + current = propagate_orbital_elements(current, dt_burn, earth->mass); + + if (current.semi_major_axis < a_prev) monotonic = false; + a_prev = current.semi_major_axis; + + if (current.eccentricity > max_e) max_e = current.eccentricity; + if (current.eccentricity < min_e) min_e = current.eccentricity; + } + + const double final_a = current.semi_major_axis; + const double total_change = final_a - initial_a; + + // Check deviation from linear trend + double max_dev = 0.0; + OrbitalElements cur = low_thrust->orbit; + a_prev = initial_a; + for (int i = 0; i < 100; i++) { + Vec3 pos, vel; + orbital_elements_to_cartesian(cur, earth->mass, &pos, &vel); + Vec3 dir = get_burn_direction_vector(BURN_PROGRADE, pos, vel); + vel = vec3_add(vel, vec3_scale(dir, dv_per)); + cur = cartesian_to_orbital_elements(pos, vel, earth->mass); + cur = propagate_orbital_elements(cur, dt_burn, earth->mass); + + const double expected = initial_a + (i + 1) * (total_change / 100.0); + const double dev = fabs(cur.semi_major_axis - expected); + if (dev > max_dev) max_dev = dev; + } + + INFO("monotonic: " << (monotonic ? "yes" : "no")); + INFO("max_ecc: " << max_e); + INFO("total_a_change: " << total_change); + INFO("max_deviation: " << max_dev); + INFO("max_deviation_pct: " << (max_dev / total_change * 100.0) << "%"); + + REQUIRE(monotonic); + REQUIRE(max_e < 0.1); + REQUIRE_THAT(max_dev, WithinAbs(0.0, total_change * 0.5)); + } + + destroy_simulation(sim); +} diff --git a/tests/test_hybrid_burns.toml b/tests/test_hybrid_burns.toml new file mode 100644 index 0000000..eb117df --- /dev/null +++ b/tests/test_hybrid_burns.toml @@ -0,0 +1,131 @@ +[[bodies]] +name = "Sun" +mass = 1.989e30 +radius = 6.96e8 +parent_index = -1 +color = { r = 1.0, g = 1.0, b = 0.0 } +orbit = { semi_major_axis = 0.0, eccentricity = 0.0, true_anomaly = 0.0 } + +[[bodies]] +name = "Earth" +mass = 5.972e24 +radius = 6.371e6 +parent_index = 0 +color = { r = 0.0, g = 0.5, b = 1.0 } +orbit = { semi_major_axis = 1.496e11, eccentricity = 0.0, true_anomaly = 0.0 } + +[[spacecraft]] +name = "Hohmann_Transfer" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 6.771e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[maneuvers]] +name = "hohmann_burn_1" +spacecraft_name = "Hohmann_Transfer" +trigger_type = "time" +trigger_value = 0.0 +direction = "prograde" +delta_v = 2440.0 + +[[maneuvers]] +name = "hohmann_burn_2" +spacecraft_name = "Hohmann_Transfer" +trigger_type = "time" +trigger_value = 5400.0 +direction = "prograde" +delta_v = 1500.0 + +[[spacecraft]] +name = "Plane_Change" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 7.0e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.2, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[maneuvers]] +name = "plane_change_burn" +spacecraft_name = "Plane_Change" +trigger_type = "time" +trigger_value = 0.0 +direction = "normal" +delta_v = 1400.0 + +[[spacecraft]] +name = "Periapsis_Burn" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 1.5e7, eccentricity = 0.5, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[maneuvers]] +name = "periapsis_burn" +spacecraft_name = "Periapsis_Burn" +trigger_type = "time" +trigger_value = 0.0 +direction = "prograde" +delta_v = 500.0 + +[[spacecraft]] +name = "Apoapsis_Burn" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 1.5e7, eccentricity = 0.5, true_anomaly = 3.141592653589793, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[maneuvers]] +name = "apoapsis_burn" +spacecraft_name = "Apoapsis_Burn" +trigger_type = "time" +trigger_value = 0.0 +direction = "prograde" +delta_v = 500.0 + +[[spacecraft]] +name = "Small_Delta_v" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 7.0e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[maneuvers]] +name = "small_burn" +spacecraft_name = "Small_Delta_v" +trigger_type = "time" +trigger_value = 0.0 +direction = "prograde" +delta_v = 0.5 + +[[spacecraft]] +name = "Large_Delta_v" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 7.0e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[maneuvers]] +name = "large_burn" +spacecraft_name = "Large_Delta_v" +trigger_type = "time" +trigger_value = 0.0 +direction = "prograde" +delta_v = 12000.0 + +[[spacecraft]] +name = "Low_Thrust_Ion" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 6.771e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[spacecraft]] +name = "Multi_Burn_Sequence" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 7.0e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[spacecraft]] +name = "Mode_Transition" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 1.2e7, eccentricity = 0.3, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[spacecraft]] +name = "Energy_Conservation" +mass = 1000.0 +parent_index = 1 +orbit = { semi_major_axis = 8.0e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }