diff --git a/scripts/precalc_hybrid_burns.py b/scripts/precalc_hybrid_burns.py index 9f69ce3..e759110 100644 --- a/scripts/precalc_hybrid_burns.py +++ b/scripts/precalc_hybrid_burns.py @@ -310,6 +310,7 @@ def main(): 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) @@ -323,6 +324,7 @@ def main(): 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}") @@ -502,6 +504,42 @@ def main(): 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() diff --git a/tests/test_hybrid_burns.cpp b/tests/test_hybrid_burns.cpp index 06c81b8..cc9df33 100644 --- a/tests/test_hybrid_burns.cpp +++ b/tests/test_hybrid_burns.cpp @@ -68,7 +68,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") } execute_maneuver(m, craft, sim, sim->time); REQUIRE(m->executed); - REQUIRE_THAT(m->executed_time, WithinAbs(sim->time, 0.001)); + REQUIRE_THAT(m->executed_time, WithinAbs(sim->time, M_TOL)); }; // Shared fixtures @@ -100,7 +100,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") exec_by_name("hohmann_burn_1", hohmann); const double v_after = vec3_magnitude(hohmann->local_velocity); - REQUIRE(v_after > v_before); + REQUIRE_THAT(v_after, WithinAbs(10112.490413, V_TOL)); const auto post_els = cartesian_to_orbital_elements( hohmann->local_position, hohmann->local_velocity, earth->mass); @@ -110,8 +110,8 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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)); + REQUIRE_THAT(post_els.semi_major_axis, WithinAbs(25762376.160113, A_TOL)); + REQUIRE_THAT(post_els.eccentricity, WithinAbs(0.737174864697325, E_TOL)); } SECTION("second burn at apogee circularizes orbit") { @@ -137,9 +137,8 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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); + REQUIRE_THAT(final_els.semi_major_axis, WithinAbs(46176507.362571, A_TOL)); + REQUIRE_THAT(final_els.eccentricity, WithinAbs(0.030811231156453, E_TOL)); } SECTION("large prograde burn produces hyperbolic orbit") { @@ -156,7 +155,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("v_escape: " << v_escape << " m/s"); INFO("v_after: " << v_after << " m/s"); - REQUIRE(v_after > v_escape); + REQUIRE_THAT(v_after, WithinAbs(19545.946840, V_TOL)); const auto hyper_els = cartesian_to_orbital_elements( large_dv->local_position, large_dv->local_velocity, earth->mass); @@ -164,8 +163,8 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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)); + REQUIRE_THAT(hyper_els.eccentricity, WithinAbs(5.709434906871548, E_TOL)); + REQUIRE_THAT(hyper_els.semi_major_axis, WithinAbs(-1486377.906994, A_TOL)); } SECTION("large burn satisfies vis-viva equation") { @@ -184,7 +183,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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)); + REQUIRE_THAT(err, WithinAbs(0.0, D_TOL)); } SECTION("prograde burn increases total energy") { @@ -213,10 +212,10 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("dE_actual: " << dE_actual); INFO("dE_expected: " << dE_expected); - REQUIRE(E_final > E_init); + REQUIRE_THAT(E_final, WithinAbs(-7735877962.552383, A_TOL)); const double dE_err = fabs(dE_actual - dE_expected) / fabs(dE_expected); - REQUIRE_THAT(dE_err, WithinAbs(0.0, 1e-12)); + REQUIRE_THAT(dE_err, WithinAbs(0.0, D_TOL)); } SECTION("retrograde burn decreases total energy") { @@ -245,10 +244,10 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("dE_actual: " << dE_actual); INFO("dE_expected: " << dE_expected); - REQUIRE(E_final < E_init); + REQUIRE_THAT(E_final, WithinAbs(-36606044984.248001, A_TOL)); const double dE_err = fabs(dE_actual - dE_expected) / fabs(dE_expected); - REQUIRE_THAT(dE_err, WithinAbs(0.0, 1e-12)); + REQUIRE_THAT(dE_err, WithinAbs(0.0, D_TOL)); } SECTION("orbital elements -> Cartesian -> burn -> orbital elements") { @@ -306,14 +305,11 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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); + REQUIRE_THAT(after_1.semi_major_axis, WithinAbs(25762376.160113, A_TOL)); // Propagate to apogee auto apogee_els = after_1; @@ -330,8 +326,8 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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); + REQUIRE_THAT(after_2.semi_major_axis, WithinAbs(46176507.362571, A_TOL)); + REQUIRE_THAT(after_2.eccentricity, WithinAbs(0.030811231156453, E_TOL)); } SECTION("three-burn sequence with plane change") { @@ -357,8 +353,8 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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); + REQUIRE_THAT(after_3.semi_major_axis, WithinAbs(8383687.781504, A_TOL)); + REQUIRE_THAT(after_3.inclination, WithinAbs(0.036692041490386, ANG_TOL)); } SECTION("prograde and retrograde are opposite") { @@ -405,10 +401,9 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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); + REQUIRE_THAT(final_els.semi_major_axis, WithinAbs(6951054.544051, A_TOL)); 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; @@ -419,21 +414,21 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("relative_error: " << rel_err); REQUIRE_THAT(rel_err, WithinAbs(0.0, 0.01)); - REQUIRE(final_els.eccentricity < 0.01); + REQUIRE_THAT(final_els.eccentricity, WithinAbs(0.003347573985440, E_TOL)); } 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); + REQUIRE_THAT(after_1.semi_major_axis, WithinAbs(7094118.510013, A_TOL)); 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); + REQUIRE_THAT(final_els.semi_major_axis, WithinAbs(7237952.003198, A_TOL)); 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); @@ -456,7 +451,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") 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); + REQUIRE_THAT(final_els.semi_major_axis, WithinAbs(13012778.714495, A_TOL)); 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); @@ -464,7 +459,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("energy_change: " << energy_change); - REQUIRE(fabs(energy_change) > 0.0); + REQUIRE_THAT(energy_change, WithinAbs(1292584.077011, A_TOL)); } SECTION("continuous burn energy increases monotonically") { @@ -494,7 +489,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("expected_approx: " << expected_approx); INFO("relative_error: " << rel_err); - REQUIRE(total_dE > 0.0); + REQUIRE_THAT(total_dE, WithinAbs(1048578803.759296, A_TOL)); REQUIRE_THAT(rel_err, WithinAbs(0.0, 0.01)); } @@ -517,7 +512,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("v_difference: " << v_diff); REQUIRE_THAT(rel_diff, WithinAbs(0.0, 1.0)); - REQUIRE(v_diff < 2.0); + REQUIRE_THAT(v_diff, WithinAbs(1.297686, 0.5)); } SECTION("propagation during burn: path length > straight line") { @@ -560,8 +555,8 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("r_start: " << r_start); INFO("r_end: " << r_end); - REQUIRE(total_path > straight); - REQUIRE(r_end > r_start); + REQUIRE_THAT(total_path, WithinAbs(38113183.100583, A_TOL)); + REQUIRE_THAT(r_end, WithinAbs(6972172.241655, R_TOL)); // Check final radius within expected bounds const double v_init = sqrt(MU_EARTH / low_thrust->orbit.semi_major_axis); @@ -576,8 +571,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("r_peri: " << r_peri); INFO("r_apo: " << r_apo); - REQUIRE(r_end >= r_peri - 1e5); - REQUIRE(r_end <= r_apo + 1e5); + REQUIRE_THAT(r_end, WithinAbs(6972172.241655, 1e5)); } SECTION("continuous burn: semi-major axis increases monotonically") { @@ -635,7 +629,7 @@ SCENARIO("Hybrid burns: impulse + continuous burn behavior", "[hybrid][burns]") INFO("max_deviation_pct: " << (max_dev / total_change * 100.0) << "%"); REQUIRE(monotonic); - REQUIRE(max_e < 0.1); + REQUIRE_THAT(max_e, WithinAbs(0.009304764034330, 0.01)); REQUIRE_THAT(max_dev, WithinAbs(0.0, total_change * 0.5)); }