cinnaboot
/
claudes_game
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

refactor: test_hybrid_burns - impulse + continuous burn behavior 

- 14 TEST_CASEs consolidated into 1 SCENARIO with 22 SECTIONs
- Shared fixture: sim, bodies, all spacecraft, maneuvers
- Helper lambdas: init_craft, find_maneuver, exec_by_name
- All tolerances precalculated via sim_engine.py
- Named tolerance constants (A_TOL, E_TOL, V_TOL)
- 112 assertions, all passing
- Updated continue.md: no decorative comments, REQUIRE for integers
test-refactor
cinnaboot 2 months ago
parent
b69c55e793
commit
20327a75d1
4 changed files with 1284 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 3
      continue.md
  2. 507
      scripts/precalc_hybrid_burns.py
  3. 643
      tests/test_hybrid_burns.cpp
  4. 131
      tests/test_hybrid_burns.toml

3
continue.md
Unescape View File

@ -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_<test_name>.py`
- If sim_engine.py lacks a feature, **stop to notify the user what feature is missing**

507
scripts/precalc_hybrid_burns.py
Unescape View File

@ -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()

643
tests/test_hybrid_burns.cpp
Unescape View File

@ -0,0 +1,643 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#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 <cmath>
#include <cstring>
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);
}

131
tests/test_hybrid_burns.toml
Unescape View File

@ -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 }
Write Preview
Loading…
Cancel
Save