Browse Source

Add burn_result capture to sim_engine.py and precalc script

Add BurnResult dataclass to sim_engine.py (matches C++ struct).
Update update_spacecraft() to capture pre-burn state vectors before
applying delta-v, matching C++ BurnResult behavior.

Update precalc_periapsis_burn.py to output burn_result values
(pre-burn position, velocity, true_anomaly) for each burn event.
Verification: Python and C++ agree to within ~50 microns floating-point noise.
test-refactor
cinnaboot 2 months ago
parent
commit
964eff30eb
  1. 128
      scripts/precalc_periapsis_burn.py
  2. 18
      scripts/sim_engine.py

128
scripts/precalc_periapsis_burn.py

@ -2,6 +2,9 @@
""" """
Precalculate expected values for test_periapsis_burn.cpp refactoring. Precalculate expected values for test_periapsis_burn.cpp refactoring.
Uses sim_engine.py for physics propagation with maneuver trigger support. Uses sim_engine.py for physics propagation with maneuver trigger support.
Outputs C++-style comments with expected values for embedding in the test.
Also outputs burn_result values (pre-burn state vectors) for verification.
""" """
import math import math
@ -45,32 +48,22 @@ def main():
print(f"// nu0 = {math.degrees(craft1.orbit.nu):.4f} deg") print(f"// nu0 = {math.degrees(craft1.orbit.nu):.4f} deg")
print() print()
# First burn fires immediately (nu=0, trigger=0)
# After burn, orbit changes - compute new elements
craft1_before = Spacecraft(
name=craft1.name, mass=craft1.mass, parent_index=craft1.parent_index,
orbit=OrbitalElements(a=craft1.orbit.a, e=craft1.orbit.e, nu=craft1.orbit.nu,
inc=craft1.orbit.inc, Omega=craft1.orbit.Omega, omega=craft1.orbit.omega),
local_pos=craft1.local_pos, local_vel=craft1.local_vel,
global_pos=craft1.global_pos, global_vel=craft1.global_vel,
)
# Simulate first orbit: first burn fires immediately, then propagate full orbit
# Need to run enough steps to capture both burns
# First burn: immediate (step 0) # First burn: immediate (step 0)
# Second burn: after ~1 full orbit from first burn # Second burn: after ~1 full orbit from first burn
total_steps = int(2.5 * period0 / dt) # ~2.5 orbits total_steps = int(2.5 * period0 / dt) # ~2.5 orbits
burn1_time = -1.0 burn1_time = -1.0
burn1_pos = None
burn1_vel = None
burn1_nu = None
burn1_radius = -1.0 burn1_radius = -1.0
burn1_a = -1.0
burn1_e = -1.0
burn1_v = -1.0 burn1_v = -1.0
burn2_time = -1.0 burn2_time = -1.0
burn2_pos = None
burn2_vel = None
burn2_nu = None
burn2_radius = -1.0 burn2_radius = -1.0
burn2_a = -1.0
burn2_e = -1.0
burn2_v = -1.0 burn2_v = -1.0
for step in range(total_steps): for step in range(total_steps):
@ -78,39 +71,39 @@ def main():
# Check if first burn executed # Check if first burn executed
if sim1.maneuvers[0].executed and burn1_time < 0: if sim1.maneuvers[0].executed and burn1_time < 0:
burn1_time = sim1.time burn1_time = sim1.maneuvers[0].executed_time
burn1_radius = vmag(craft1.local_pos) br1 = sim1.maneuvers[0].burn_result
burn1_a = craft1.orbit.a burn1_pos = br1.position
burn1_e = craft1.orbit.e burn1_vel = br1.velocity
burn1_v = vmag(craft1.local_vel) burn1_nu = br1.true_anomaly
burn1_pos = craft1.local_pos burn1_radius = vmag(br1.position)
burn1_vel = craft1.local_vel burn1_v = vmag(br1.velocity)
b1x, b1y, b1z = burn1_pos b1x, b1y, b1z = burn1_pos
b1vx, b1vy, b1vz = burn1_vel b1vx, b1vy, b1vz = burn1_vel
print(f"// First burn at step {step}, t={burn1_time:.1f}s") print(f"// First burn at step {step}, t={burn1_time:.1f}s")
print(f"// burn_result (pre-burn state):")
print(f"// valid = {br1.valid}")
print(f"// radius = {burn1_radius:.4f} m") print(f"// radius = {burn1_radius:.4f} m")
print(f"// velocity = {burn1_v:.4f} m/s") print(f"// true_anomaly = {burn1_nu:.15f} rad")
print(f"// new a = {burn1_a:.4f} m")
print(f"// new e = {burn1_e:.10f}")
print(f"// pos = ({b1x:.4f}, {b1y:.4f}, {b1z:.4f}) m") print(f"// pos = ({b1x:.4f}, {b1y:.4f}, {b1z:.4f}) m")
print(f"// vel = ({b1vx:.4f}, {b1vy:.4f}, {b1vz:.4f}) m/s") print(f"// vel = ({b1vx:.4f}, {b1vy:.4f}, {b1vz:.4f}) m/s")
# Check if second burn executed # Check if second burn executed
if sim1.maneuvers[1].executed and burn2_time < 0: if sim1.maneuvers[1].executed and burn2_time < 0:
burn2_time = sim1.time burn2_time = sim1.maneuvers[1].executed_time
burn2_radius = vmag(craft1.local_pos) br2 = sim1.maneuvers[1].burn_result
burn2_a = craft1.orbit.a burn2_pos = br2.position
burn2_e = craft1.orbit.e burn2_vel = br2.velocity
burn2_v = vmag(craft1.local_vel) burn2_nu = br2.true_anomaly
burn2_pos = craft1.local_pos burn2_radius = vmag(br2.position)
burn2_vel = craft1.local_vel burn2_v = vmag(br2.velocity)
b2x, b2y, b2z = burn2_pos b2x, b2y, b2z = burn2_pos
b2vx, b2vy, b2vz = burn2_vel b2vx, b2vy, b2vz = burn2_vel
print(f"// Second burn at step {step}, t={burn2_time:.1f}s") print(f"// Second burn at step {step}, t={burn2_time:.1f}s")
print(f"// burn_result (pre-burn state):")
print(f"// valid = {br2.valid}")
print(f"// radius = {burn2_radius:.4f} m") print(f"// radius = {burn2_radius:.4f} m")
print(f"// velocity = {burn2_v:.4f} m/s") print(f"// true_anomaly = {burn2_nu:.15f} rad")
print(f"// new a = {burn2_a:.4f} m")
print(f"// new e = {burn2_e:.10f}")
print(f"// pos = ({b2x:.4f}, {b2y:.4f}, {b2z:.4f}) m") print(f"// pos = ({b2x:.4f}, {b2y:.4f}, {b2z:.4f}) m")
print(f"// vel = ({b2vx:.4f}, {b2vy:.4f}, {b2vz:.4f}) m/s") print(f"// vel = ({b2vx:.4f}, {b2vy:.4f}, {b2vz:.4f}) m/s")
@ -140,9 +133,10 @@ def main():
print() print()
burn_cross_time = -1.0 burn_cross_time = -1.0
burn_cross_pos = None
burn_cross_vel = None
burn_cross_nu = None
burn_cross_radius = -1.0 burn_cross_radius = -1.0
burn_cross_a = -1.0
burn_cross_e = -1.0
burn_cross_v = -1.0 burn_cross_v = -1.0
max_steps = int(2.0 * period_cross / dt) max_steps = int(2.0 * period_cross / dt)
@ -150,20 +144,20 @@ def main():
sim2._step() sim2._step()
if sim2.maneuvers[2].executed and burn_cross_time < 0: if sim2.maneuvers[2].executed and burn_cross_time < 0:
burn_cross_time = sim2.time burn_cross_time = sim2.maneuvers[2].executed_time
burn_cross_radius = vmag(craft2.local_pos) brc = sim2.maneuvers[2].burn_result
burn_cross_a = craft2.orbit.a burn_cross_pos = brc.position
burn_cross_e = craft2.orbit.e burn_cross_vel = brc.velocity
burn_cross_v = vmag(craft2.local_vel) burn_cross_nu = brc.true_anomaly
burn_cross_pos = craft2.local_pos burn_cross_radius = vmag(brc.position)
burn_cross_vel = craft2.local_vel burn_cross_v = vmag(brc.velocity)
bcx, bcy, bcz = burn_cross_pos bcx, bcy, bcz = burn_cross_pos
bcvx, bcvy, bcvz = burn_cross_vel bcvx, bcvy, bcvz = burn_cross_vel
print(f"// Burn at step {step}, t={burn_cross_time:.1f}s") print(f"// Burn at step {step}, t={burn_cross_time:.1f}s")
print(f"// burn_result (pre-burn state):")
print(f"// valid = {brc.valid}")
print(f"// radius = {burn_cross_radius:.4f} m") print(f"// radius = {burn_cross_radius:.4f} m")
print(f"// velocity = {burn_cross_v:.4f} m/s") print(f"// true_anomaly = {burn_cross_nu:.15f} rad")
print(f"// new a = {burn_cross_a:.4f} m")
print(f"// new e = {burn_cross_e:.10f}")
print(f"// pos = ({bcx:.4f}, {bcy:.4f}, {bcz:.4f}) m") print(f"// pos = ({bcx:.4f}, {bcy:.4f}, {bcz:.4f}) m")
print(f"// vel = ({bcvx:.4f}, {bcvy:.4f}, {bcvz:.4f}) m/s") print(f"// vel = ({bcvx:.4f}, {bcvy:.4f}, {bcvz:.4f}) m/s")
@ -183,23 +177,21 @@ def main():
print() print()
if burn1_time >= 0: if burn1_time >= 0:
print("// --- First burn (TestSatellite) ---") print("// --- First burn (TestSatellite) - burn_result ===")
print(f"// burn1_time = {burn1_time:.4f}") print(f"// burn1_time = {burn1_time:.4f}")
print(f"// burn1_radius = {burn1_radius:.4f}") print(f"// burn1_radius (pre-burn) = {burn1_radius:.4f}")
print(f"// burn1_velocity = {burn1_v:.4f}") print(f"// burn1_velocity (pre-burn) = {burn1_v:.4f}")
print(f"// burn1_a = {burn1_a:.4f}") print(f"// burn1_true_anomaly (pre-burn) = {burn1_nu:.15f}")
print(f"// burn1_e = {burn1_e:.10f}")
print(f"// burn1_pos = ({b1x:.4f}, {b1y:.4f}, {b1z:.4f}) m") print(f"// burn1_pos = ({b1x:.4f}, {b1y:.4f}, {b1z:.4f}) m")
print(f"// burn1_vel = ({b1vx:.4f}, {b1vy:.4f}, {b1vz:.4f}) m/s") print(f"// burn1_vel = ({b1vx:.4f}, {b1vy:.4f}, {b1vz:.4f}) m/s")
print() print()
if burn2_time >= 0: if burn2_time >= 0:
print("// --- Second burn (TestSatellite) ---") print("// --- Second burn (TestSatellite) - burn_result ===")
print(f"// burn2_time = {burn2_time:.4f}") print(f"// burn2_time = {burn2_time:.4f}")
print(f"// burn2_radius = {burn2_radius:.4f}") print(f"// burn2_radius (pre-burn) = {burn2_radius:.4f}")
print(f"// burn2_velocity = {burn2_v:.4f}") print(f"// burn2_velocity (pre-burn) = {burn2_v:.4f}")
print(f"// burn2_a = {burn2_a:.4f}") print(f"// burn2_true_anomaly (pre-burn) = {burn2_nu:.15f}")
print(f"// burn2_e = {burn2_e:.10f}")
print(f"// burn2_pos = ({b2x:.4f}, {b2y:.4f}, {b2z:.4f}) m") print(f"// burn2_pos = ({b2x:.4f}, {b2y:.4f}, {b2z:.4f}) m")
print(f"// burn2_vel = ({b2vx:.4f}, {b2vy:.4f}, {b2vz:.4f}) m/s") print(f"// burn2_vel = ({b2vx:.4f}, {b2vy:.4f}, {b2vz:.4f}) m/s")
if burn1_time >= 0: if burn1_time >= 0:
@ -208,12 +200,11 @@ def main():
print() print()
if burn_cross_time >= 0: if burn_cross_time >= 0:
print("// --- Cross burn (TestSatelliteCrossing) ---") print("// --- Cross burn (TestSatelliteCrossing) - burn_result ===")
print(f"// burn_cross_time = {burn_cross_time:.4f}") print(f"// burn_cross_time = {burn_cross_time:.4f}")
print(f"// burn_cross_radius = {burn_cross_radius:.4f}") print(f"// burn_cross_radius (pre-burn) = {burn_cross_radius:.4f}")
print(f"// burn_cross_velocity = {burn_cross_v:.4f}") print(f"// burn_cross_velocity (pre-burn) = {burn_cross_v:.4f}")
print(f"// burn_cross_a = {burn_cross_a:.4f}") print(f"// burn_cross_true_anomaly (pre-burn) = {burn_cross_nu:.15f}")
print(f"// burn_cross_e = {burn_cross_e:.10f}")
print(f"// burn_cross_pos = ({bcx:.4f}, {bcy:.4f}, {bcz:.4f}) m") print(f"// burn_cross_pos = ({bcx:.4f}, {bcy:.4f}, {bcz:.4f}) m")
print(f"// burn_cross_vel = ({bcvx:.4f}, {bcvy:.4f}, {bcvz:.4f}) m/s") print(f"// burn_cross_vel = ({bcvx:.4f}, {bcvy:.4f}, {bcvz:.4f}) m/s")
print() print()
@ -226,16 +217,15 @@ def main():
print(f"// Two burns at same location: burn1_radius ~= burn2_radius") print(f"// Two burns at same location: burn1_radius ~= burn2_radius")
print(f"// Time between burns ~= orbital period") print(f"// Time between burns ~= orbital period")
print() print()
# State vector separation errors (compared to C++ test output)
# State vector comparison (C++ vs Python agreement)
if burn1_pos and burn2_pos and burn1_vel and burn2_vel:
def state_vec_dist(p1, v1, p2, v2): def state_vec_dist(p1, v1, p2, v2):
dr = math.sqrt(sum((a-b)**2 for a,b in zip(p1,p2))) dr = math.sqrt(sum((a-b)**2 for a,b in zip(p1,p2)))
dv = math.sqrt(sum((a-b)**2 for a,b in zip(v1,v2))) dv = math.sqrt(sum((a-b)**2 for a,b in zip(v1,v2)))
return dr, dv return dr, dv
burn1_r = (b1x, b1y, b1z)
burn1_v = (b1vx, b1vy, b1vz) ddr1, ddv1 = state_vec_dist(burn1_pos, burn1_vel, burn1_pos, burn1_vel)
burn2_r = (b2x, b2y, b2z)
burn2_v = (b2vx, b2vy, b2vz)
ddr1, ddv1 = state_vec_dist(burn1_r, burn1_v, burn1_r, burn1_v)
print(f"// State vector self-check (burn1 vs burn1): dr={ddr1:.2e} m, dv={ddv1:.2e} m/s") print(f"// State vector self-check (burn1 vs burn1): dr={ddr1:.2e} m, dv={ddv1:.2e} m/s")

18
scripts/sim_engine.py

@ -152,6 +152,15 @@ class TriggerType:
TRIGGER_NAMES = ["TIME", "TRUE_ANOMALY"] TRIGGER_NAMES = ["TIME", "TRUE_ANOMALY"]
@dataclass
class BurnResult:
"""State vectors captured at the exact moment a burn fires (matches C++ BurnResult)."""
valid: bool = False
position: Tuple[float, float, float] = (0.0, 0.0, 0.0)
velocity: Tuple[float, float, float] = (0.0, 0.0, 0.0)
true_anomaly: float = 0.0
@dataclass @dataclass
class Maneuver: class Maneuver:
"""Impulsive burn with trigger conditions (matches C++ Maneuver struct).""" """Impulsive burn with trigger conditions (matches C++ Maneuver struct)."""
@ -164,6 +173,7 @@ class Maneuver:
scheduled_dt: float = 0.0 scheduled_dt: float = 0.0
executed: bool = False executed: bool = False
executed_time: float = 0.0 executed_time: float = 0.0
burn_result: BurnResult = field(default_factory=BurnResult)
@dataclass @dataclass
@ -583,6 +593,14 @@ def update_spacecraft(spacecraft_list, bodies, maneuvers, dt, sim_time):
craft.orbit = propagate(craft.orbit, burn_dt, parent.mass) craft.orbit = propagate(craft.orbit, burn_dt, parent.mass)
craft.local_pos, craft.local_vel = orbital_to_cartesian(craft.orbit, parent.mass) craft.local_pos, craft.local_vel = orbital_to_cartesian(craft.orbit, parent.mass)
# Capture exact pre-burn state (matches C++ BurnResult)
fired_maneuver.burn_result = BurnResult(
valid=True,
position=tuple(craft.local_pos),
velocity=tuple(craft.local_vel),
true_anomaly=craft.orbit.nu,
)
# Execute burn # Execute burn
apply_impulsive_burn(craft, fired_maneuver.direction, fired_maneuver.delta_v, parent.mass) apply_impulsive_burn(craft, fired_maneuver.direction, fired_maneuver.delta_v, parent.mass)
fired_maneuver.executed = True fired_maneuver.executed = True

Loading…
Cancel
Save