cinnaboot
/
claudes_game
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
vibe coding an orbital mechanics simulation to try out claude code
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
515 Commits
2 Branches
0 Tags
1.9 MiB
 
 
 
 
 
Branch: test-refactor
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'test-refactor'
${ noResults }
claudes_game/scripts/precalc_hybrid_burns.py

545 lines
22 KiB
Raw Permalink Blame History

#!/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)
a_after_1_mb = orbit_after_1.a
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"// after_1_a = {a_after_1_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()
# Test: Three-burn sequence with plane change
sim_3b = Simulator("tests/test_hybrid_burns.toml", dt=dt)
craft_3b = sim_3b.spacecraft[0] # Hohmann_Transfer
earth_3b = sim_3b.bodies[1]
pos_3b, vel_3b = orbital_to_cartesian(craft_3b.orbit, earth_3b.mass)
craft_3b.local_pos = pos_3b
craft_3b.local_vel = vel_3b
init_a_3b = craft_3b.orbit.a
init_inc_3b = craft_3b.orbit.inc
# Burn 1: prograde 500 m/s
burn1_dir = get_burn_direction(BurnDirection.PROGRADE, craft_3b.local_pos, craft_3b.local_vel)
craft_3b.local_vel = vadd(craft_3b.local_vel, vscale(burn1_dir, 500.0))
craft_3b.orbit = cartesian_to_orbital_elements(craft_3b.local_pos, craft_3b.local_vel, earth_3b.mass)
# Burn 2: normal 300 m/s
burn2_dir = get_burn_direction(BurnDirection.NORMAL, craft_3b.local_pos, craft_3b.local_vel)
craft_3b.local_vel = vadd(craft_3b.local_vel, vscale(burn2_dir, 300.0))
craft_3b.orbit = cartesian_to_orbital_elements(craft_3b.local_pos, craft_3b.local_vel, earth_3b.mass)
# Burn 3: prograde 200 m/s
burn3_dir = get_burn_direction(BurnDirection.PROGRADE, craft_3b.local_pos, craft_3b.local_vel)
craft_3b.local_vel = vadd(craft_3b.local_vel, vscale(burn3_dir, 200.0))
craft_3b.orbit = cartesian_to_orbital_elements(craft_3b.local_pos, craft_3b.local_vel, earth_3b.mass)
final_a_3b = craft_3b.orbit.a
final_inc_3b = craft_3b.orbit.inc
print("// === Three-burn sequence with plane change ===")
print(f"// init_a = {init_a_3b:.6f}")
print(f"// init_inc = {init_inc_3b:.15f}")
print(f"// final_a = {final_a_3b:.6f}")
print(f"// final_inc = {final_inc_3b:.15f}")
print()
if __name__ == "__main__":
main()