Browse Source

Add test_hybrid_impulse_burns for burn handling validation

- Tests Hohmann transfers (2 burns: apogee raise + circularization)
- Tests plane changes at nodes with normal burns
- Tests impulsive burns at periapsis and apoapsis
- Tests minimal burns (Δv < 1 m/s) for numerical precision
- Tests large burns (Δv > orbital velocity) for hyperbolic orbits
- Validates energy conservation during burns
- Validates round-trip conversion with burns
- Validates multiple burn sequences
- Critical for analytical propagation switch: validates burn handling workflow
main
cinnaboot 5 months ago
parent
commit
e132751c25
  1. 179
      tests/configs/test_hybrid_impulse_burns.toml
  2. 426
      tests/test_hybrid_impulse_burns.cpp

179
tests/configs/test_hybrid_impulse_burns.toml

@ -0,0 +1,179 @@
# Test Configuration: Hybrid Impulse Burns for Analytical Propagation
# Sun + Earth system with multiple spacecraft for impulsive maneuver testing
# Tests the critical workflow: orbital elements → Cartesian → burn → orbital elements
[[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
}
# 1. Hohmann Transfer Spacecraft
# Initial circular LEO orbit (altitude ~400 km)
# Two maneuvers: apogee raise, circularization
[[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
# 2. Plane Change Spacecraft
# Initial circular orbit with inclination 0.2 rad
# One maneuver: normal burn at ascending node to change inclination to 0.4 rad
[[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
# 3. Periapsis Burn Spacecraft
# Initial elliptical orbit (e = 0.5)
# One maneuver: prograde burn at periapsis to raise apoapsis
[[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
# 4. Apoapsis Burn Spacecraft
# Initial elliptical orbit (e = 0.5)
# One maneuver: prograde burn at apoapsis to raise periapsis
[[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
# 5. Small Delta-v Burn Spacecraft
# Initial circular orbit
# One maneuver: minimal prograde burn (Δv < 1 m/s)
[[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
# 6. Large Delta-v Burn Spacecraft
# Initial circular orbit
# One maneuver: large prograde burn (Δv > orbital velocity)
[[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

426
tests/test_hybrid_impulse_burns.cpp

@ -0,0 +1,426 @@
#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/spacecraft.h"
#include "../src/maneuver.h"
#include "../src/config_loader.h"
#include "../src/test_utilities.h"
#include <cmath>
const double POSITION_TOLERANCE = 1e-3;
const double VELOCITY_TOLERANCE = 1e-3;
const double ELEMENT_TOLERANCE = 1e-6;
const double ENERGY_TOLERANCE = 1e-6;
TEST_CASE("Config loading for hybrid impulse burns", "[hybrid][impulse][config]") {
const double TIME_STEP = 60.0;
SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/configs/test_hybrid_impulse_burns.toml"));
REQUIRE(sim->body_count == 2);
REQUIRE(std::string(sim->bodies[0].name) == "Sun");
REQUIRE(std::string(sim->bodies[1].name) == "Earth");
REQUIRE(sim->craft_count == 6);
REQUIRE(std::string(sim->spacecraft[0].name) == "Hohmann_Transfer");
REQUIRE(sim->spacecraft[0].parent_index == 1);
REQUIRE(std::string(sim->spacecraft[1].name) == "Plane_Change");
REQUIRE(sim->spacecraft[1].parent_index == 1);
REQUIRE(std::string(sim->spacecraft[2].name) == "Periapsis_Burn");
REQUIRE(sim->spacecraft[2].parent_index == 1);
REQUIRE(std::string(sim->spacecraft[3].name) == "Apoapsis_Burn");
REQUIRE(sim->spacecraft[3].parent_index == 1);
REQUIRE(std::string(sim->spacecraft[4].name) == "Small_Delta_v");
REQUIRE(sim->spacecraft[4].parent_index == 1);
REQUIRE(std::string(sim->spacecraft[5].name) == "Large_Delta_v");
REQUIRE(sim->spacecraft[5].parent_index == 1);
REQUIRE(sim->maneuver_count == 7);
destroy_simulation(sim);
}
SCENARIO("Hohmann transfer with two burns", "[hybrid][impulse][hohmann]") {
const double TIME_STEP = 60.0;
SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/configs/test_hybrid_impulse_burns.toml"));
Spacecraft* craft = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[1];
Vec3 initial_pos;
Vec3 initial_vel;
orbital_elements_to_cartesian(craft->orbit, earth->mass, &initial_pos, &initial_vel);
SECTION("First burn at perigee raises apogee") {
OrbitalElements initial_elements = craft->orbit;
double initial_velocity_mag = vec3_magnitude(initial_vel);
apply_impulsive_burn(craft, BURN_PROGRADE, 2440.0);
double new_velocity_mag = vec3_magnitude(craft->local_velocity);
REQUIRE(new_velocity_mag > initial_velocity_mag);
Vec3 new_pos = craft->local_position;
Vec3 new_vel = craft->local_velocity;
OrbitalElements new_elements = cartesian_to_orbital_elements(new_pos, new_vel, earth->mass);
INFO("Initial a: " << initial_elements.semi_major_axis);
INFO("New a: " << new_elements.semi_major_axis);
INFO("Initial e: " << initial_elements.eccentricity);
INFO("New e: " << new_elements.eccentricity);
REQUIRE(new_elements.semi_major_axis > initial_elements.semi_major_axis);
REQUIRE(new_elements.eccentricity > initial_elements.eccentricity);
}
destroy_simulation(sim);
}
SCENARIO("Large burns (Δv > orbital velocity)", "[hybrid][impulse][large_delta_v]") {
const double TIME_STEP = 60.0;
SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/configs/test_hybrid_impulse_burns.toml"));
Spacecraft* craft = &sim->spacecraft[5];
CelestialBody* earth = &sim->bodies[1];
Vec3 initial_pos;
Vec3 initial_vel;
orbital_elements_to_cartesian(craft->orbit, earth->mass, &initial_pos, &initial_vel);
OrbitalElements initial_elements = cartesian_to_orbital_elements(initial_pos, initial_vel, earth->mass);
double initial_velocity_mag = vec3_magnitude(initial_vel);
double escape_velocity = sqrt(2.0 * G * earth->mass / vec3_magnitude(initial_pos));
SECTION("Large prograde burn produces hyperbolic orbit") {
INFO("Initial velocity: " << initial_velocity_mag << " m/s");
INFO("Escape velocity: " << escape_velocity << " m/s");
double delta_v = 12000.0;
apply_impulsive_burn(craft, BURN_PROGRADE, delta_v);
double final_velocity_mag = vec3_magnitude(craft->local_velocity);
INFO("Final velocity: " << final_velocity_mag << " m/s");
REQUIRE(final_velocity_mag > escape_velocity);
OrbitalElements new_elements = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
INFO("Initial e: " << initial_elements.eccentricity);
INFO("New e: " << new_elements.eccentricity);
REQUIRE(new_elements.eccentricity > 1.0);
REQUIRE(new_elements.semi_major_axis < 0);
}
SECTION("Large burn produces correct hyperbolic trajectory") {
apply_impulsive_burn(craft, BURN_PROGRADE, 12000.0);
OrbitalElements new_elements = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
double final_velocity_mag = vec3_magnitude(craft->local_velocity);
double r = vec3_magnitude(craft->local_position);
double vis_viva_expected = final_velocity_mag * final_velocity_mag;
double vis_viva_calculated = G * earth->mass * (2.0 / r - 1.0 / new_elements.semi_major_axis);
INFO("Vis-viva expected: " << vis_viva_expected);
INFO("Vis-viva calculated: " << vis_viva_calculated);
double vis_viva_error = fabs(vis_viva_expected - vis_viva_calculated) / vis_viva_expected;
REQUIRE(vis_viva_error < 1e-6);
}
destroy_simulation(sim);
}
SCENARIO("Energy conservation during burns", "[hybrid][impulse][energy]") {
const double TIME_STEP = 60.0;
SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/configs/test_hybrid_impulse_burns.toml"));
Spacecraft* craft = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[1];
Vec3 initial_pos;
Vec3 initial_vel;
orbital_elements_to_cartesian(craft->orbit, earth->mass, &initial_pos, &initial_vel);
craft->local_position = initial_pos;
craft->local_velocity = initial_vel;
double initial_ke = 0.5 * craft->mass * vec3_dot(craft->local_velocity, craft->local_velocity);
double initial_pe = -G * craft->mass * earth->mass / vec3_magnitude(craft->local_position);
double initial_total_energy = initial_ke + initial_pe;
SECTION("Prograde burn increases total energy") {
double delta_v = 1000.0;
Vec3 v_initial = craft->local_velocity;
apply_impulsive_burn(craft, BURN_PROGRADE, delta_v);
Vec3 v_final = craft->local_velocity;
Vec3 dv = vec3_sub(v_final, v_initial);
double expected_energy_change = vec3_dot(v_initial, dv) * craft->mass + 0.5 * craft->mass * vec3_dot(dv, dv);
double final_ke = 0.5 * craft->mass * vec3_dot(craft->local_velocity, craft->local_velocity);
double final_pe = -G * craft->mass * earth->mass / vec3_magnitude(craft->local_position);
double final_total_energy = final_ke + final_pe;
double actual_energy_change = final_total_energy - initial_total_energy;
INFO("Initial energy: " << initial_total_energy);
INFO("Final energy: " << final_total_energy);
INFO("Expected ΔE: " << expected_energy_change);
INFO("Actual ΔE: " << actual_energy_change);
REQUIRE(final_total_energy > initial_total_energy);
double energy_error = fabs(actual_energy_change - expected_energy_change) / fabs(expected_energy_change);
REQUIRE(energy_error < 1e-6);
}
SECTION("Retrograde burn decreases total energy") {
double delta_v = 1000.0;
Vec3 v_initial = craft->local_velocity;
apply_impulsive_burn(craft, BURN_RETROGRADE, delta_v);
Vec3 v_final = craft->local_velocity;
Vec3 dv = vec3_sub(v_final, v_initial);
double expected_energy_change = vec3_dot(v_initial, dv) * craft->mass + 0.5 * craft->mass * vec3_dot(dv, dv);
double final_ke = 0.5 * craft->mass * vec3_dot(craft->local_velocity, craft->local_velocity);
double final_pe = -G * craft->mass * earth->mass / vec3_magnitude(craft->local_position);
double final_total_energy = final_ke + final_pe;
double actual_energy_change = final_total_energy - initial_total_energy;
INFO("Initial energy: " << initial_total_energy);
INFO("Final energy: " << final_total_energy);
INFO("Expected ΔE: " << expected_energy_change);
INFO("Actual ΔE: " << actual_energy_change);
REQUIRE(final_total_energy < initial_total_energy);
double energy_error = fabs(actual_energy_change - expected_energy_change) / fabs(expected_energy_change);
REQUIRE(energy_error < 1e-6);
}
destroy_simulation(sim);
}
SCENARIO("Round-trip conversion with burns", "[hybrid][impulse][roundtrip]") {
const double TIME_STEP = 60.0;
SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/configs/test_hybrid_impulse_burns.toml"));
Spacecraft* craft = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[1];
SECTION("Orbital elements → Cartesian → burn → orbital elements") {
OrbitalElements original_elements = craft->orbit;
Vec3 position_from_elements;
Vec3 velocity_from_elements;
orbital_elements_to_cartesian(original_elements, earth->mass, &position_from_elements, &velocity_from_elements);
INFO("Original semi_major_axis: " << original_elements.semi_major_axis);
INFO("Original eccentricity: " << original_elements.eccentricity);
OrbitalElements recovered_elements = cartesian_to_orbital_elements(position_from_elements, velocity_from_elements, earth->mass);
INFO("Recovered semi_major_axis: " << recovered_elements.semi_major_axis);
INFO("Recovered eccentricity: " << recovered_elements.eccentricity);
REQUIRE_THAT(recovered_elements.semi_major_axis, Catch::Matchers::WithinAbs(original_elements.semi_major_axis, ELEMENT_TOLERANCE));
REQUIRE_THAT(recovered_elements.eccentricity, Catch::Matchers::WithinAbs(original_elements.eccentricity, ELEMENT_TOLERANCE));
Vec3 burn_velocity = calculate_prograde_dir(velocity_from_elements);
Vec3 new_velocity = velocity_from_elements;
new_velocity.x += burn_velocity.x * 1000.0;
new_velocity.y += burn_velocity.y * 1000.0;
new_velocity.z += burn_velocity.z * 1000.0;
OrbitalElements post_burn_elements = cartesian_to_orbital_elements(position_from_elements, new_velocity, earth->mass);
INFO("Post-burn semi_major_axis: " << post_burn_elements.semi_major_axis);
INFO("Post-burn eccentricity: " << post_burn_elements.eccentricity);
REQUIRE(post_burn_elements.semi_major_axis != recovered_elements.semi_major_axis);
REQUIRE(post_burn_elements.eccentricity != recovered_elements.eccentricity);
}
SECTION("Multiple round-trip conversions with burns") {
OrbitalElements original_elements = craft->orbit;
Vec3 position;
Vec3 velocity;
orbital_elements_to_cartesian(original_elements, earth->mass, &position, &velocity);
for (int i = 0; i < 5; i++) {
OrbitalElements elements = cartesian_to_orbital_elements(position, velocity, earth->mass);
orbital_elements_to_cartesian(elements, earth->mass, &position, &velocity);
INFO("Iteration " << i << " complete");
}
OrbitalElements final_elements = cartesian_to_orbital_elements(position, velocity, earth->mass);
INFO("Original semi_major_axis: " << original_elements.semi_major_axis);
INFO("Final semi_major_axis: " << final_elements.semi_major_axis);
INFO("Original eccentricity: " << original_elements.eccentricity);
INFO("Final eccentricity: " << final_elements.eccentricity);
double a_error = fabs(final_elements.semi_major_axis - original_elements.semi_major_axis) / original_elements.semi_major_axis;
double e_error = fabs(final_elements.eccentricity - original_elements.eccentricity);
REQUIRE(a_error < 1e-9);
REQUIRE(e_error < 1e-9);
}
destroy_simulation(sim);
}
SCENARIO("Multiple burn sequences", "[hybrid][impulse][sequence]") {
const double TIME_STEP = 60.0;
SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/configs/test_hybrid_impulse_burns.toml"));
Spacecraft* craft = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[1];
Vec3 initial_pos;
Vec3 initial_vel;
orbital_elements_to_cartesian(craft->orbit, earth->mass, &initial_pos, &initial_vel);
craft->local_position = initial_pos;
craft->local_velocity = initial_vel;
SECTION("Two-burn sequence raises orbit") {
OrbitalElements initial_elements = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
INFO("Initial a: " << initial_elements.semi_major_axis);
INFO("Initial e: " << initial_elements.eccentricity);
apply_impulsive_burn(craft, BURN_PROGRADE, 500.0);
OrbitalElements after_first_burn = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
INFO("After first burn a: " << after_first_burn.semi_major_axis);
INFO("After first burn e: " << after_first_burn.eccentricity);
REQUIRE(after_first_burn.semi_major_axis > initial_elements.semi_major_axis);
OrbitalElements apogee_elements = after_first_burn;
apogee_elements.true_anomaly = M_PI;
Vec3 apogee_pos;
Vec3 apogee_vel;
orbital_elements_to_cartesian(apogee_elements, earth->mass, &apogee_pos, &apogee_vel);
craft->local_position = apogee_pos;
craft->local_velocity = apogee_vel;
apply_impulsive_burn(craft, BURN_PROGRADE, 300.0);
OrbitalElements after_second_burn = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
INFO("After second burn a: " << after_second_burn.semi_major_axis);
INFO("After second burn e: " << after_second_burn.eccentricity);
REQUIRE(after_second_burn.semi_major_axis > after_first_burn.semi_major_axis);
REQUIRE(after_second_burn.eccentricity < after_first_burn.eccentricity);
}
SECTION("Three-burn sequence with plane change") {
OrbitalElements initial_elements = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
apply_impulsive_burn(craft, BURN_PROGRADE, 500.0);
OrbitalElements after_burn1 = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
apply_impulsive_burn(craft, BURN_NORMAL, 300.0);
OrbitalElements after_burn2 = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
apply_impulsive_burn(craft, BURN_PROGRADE, 200.0);
OrbitalElements after_burn3 = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, earth->mass);
INFO("Initial a: " << initial_elements.semi_major_axis);
INFO("After burn 3 a: " << after_burn3.semi_major_axis);
INFO("Initial inclination: " << initial_elements.inclination);
INFO("After burn 3 inclination: " << after_burn3.inclination);
REQUIRE(after_burn3.semi_major_axis > initial_elements.semi_major_axis);
REQUIRE(after_burn3.inclination > initial_elements.inclination);
}
destroy_simulation(sim);
}
TEST_CASE("Burn direction vector calculation", "[hybrid][impulse][direction]") {
const double TIME_STEP = 60.0;
SimulationState* sim = create_simulation(10, 10, 100, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/configs/test_hybrid_impulse_burns.toml"));
Spacecraft* craft = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[1];
Vec3 position;
Vec3 velocity;
orbital_elements_to_cartesian(craft->orbit, earth->mass, &position, &velocity);
SECTION("Prograde and retrograde are opposite") {
Vec3 prograde = calculate_prograde_dir(velocity);
Vec3 retrograde = calculate_retrograde_dir(velocity);
double dot_product = vec3_dot(prograde, retrograde);
INFO("Prograde · Retrograde: " << dot_product);
REQUIRE_THAT(dot_product, Catch::Matchers::WithinAbs(-1.0, 1e-6));
}
SECTION("Normal and antinormal are opposite") {
Vec3 normal = calculate_normal_dir(position, velocity);
Vec3 antinormal = calculate_antinormal_dir(position, velocity);
double dot_product = vec3_dot(normal, antinormal);
INFO("Normal · Antinormal: " << dot_product);
REQUIRE_THAT(dot_product, Catch::Matchers::WithinAbs(-1.0, 1e-6));
}
SECTION("Radial in and radial out are opposite") {
Vec3 radial_in = calculate_radial_in_dir(position);
Vec3 radial_out = calculate_radial_out_dir(position);
double dot_product = vec3_dot(radial_in, radial_out);
INFO("Radial_in · Radial_out: " << dot_product);
REQUIRE_THAT(dot_product, Catch::Matchers::WithinAbs(-1.0, 1e-6));
}
destroy_simulation(sim);
}
Loading…
Cancel
Save