diff --git a/tests/configs/hyperbolic_comet.toml b/tests/configs/hyperbolic_comet.toml new file mode 100644 index 0000000..ffd0fcc --- /dev/null +++ b/tests/configs/hyperbolic_comet.toml @@ -0,0 +1,23 @@ +# Test Configuration: Sun + Hyperbolic Comet +# Comet with hyperbolic orbit (eccentricity = 1.5) +# Hyperbolic escape trajectory - total energy > 0, velocity > escape velocity + +[[bodies]] +name = "Sun" +mass = 1.989e30 +radius = 6.96e8 +position = { x = 0.0, y = 0.0, z = 0.0 } +parent_index = -1 +color = { r = 1.0, g = 1.0, b = 0.0 } +eccentricity = 0.0 +semi_major_axis = 0.0 + +[[bodies]] +name = "HyperbolicComet" +mass = 1.0e14 +radius = 5.0e3 +position = { x = 1.496e11, y = 0.0, z = 0.0 } +parent_index = 0 +color = { r = 0.5, g = 1.0, b = 0.5 } +eccentricity = 1.5 +semi_major_axis = -1.496e11 diff --git a/tests/test_hyperbolic_orbit.cpp b/tests/test_hyperbolic_orbit.cpp new file mode 100644 index 0000000..f4c9481 --- /dev/null +++ b/tests/test_hyperbolic_orbit.cpp @@ -0,0 +1,221 @@ +#include +#include "../src/physics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include +#include + +TEST_CASE("Hyperbolic orbit - energy and escape trajectory", "[hyperbolic][energy][escape]") { + const double TIME_STEP = 60.0; + const double DAYS_TO_SIMULATE = 300.0; + const double SECONDS_PER_DAY = 86400.0; + const double AU = 1.496e11; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/configs/hyperbolic_comet.toml")); + + const int COMET_INDEX = 1; + const int SUN_INDEX = 0; + + Vec3 initial_position = sim->bodies[COMET_INDEX].position; + double initial_distance = vec3_magnitude(initial_position); + double initial_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + + double initial_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); + double initial_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], + &sim->bodies[SUN_INDEX]); + double initial_total_energy = initial_kinetic + initial_potential; + + double escape_velocity = sqrt(2.0 * G * sim->bodies[SUN_INDEX].mass / initial_distance); + + INFO("Initial distance: " << initial_distance / AU << " AU"); + INFO("Initial velocity: " << initial_velocity / 1000.0 << " km/s"); + INFO("Escape velocity: " << escape_velocity / 1000.0 << " km/s"); + INFO("Initial kinetic energy: " << initial_kinetic); + INFO("Initial potential energy: " << initial_potential); + INFO("Initial total energy: " << initial_total_energy); + + REQUIRE(initial_total_energy > 0.0); + + REQUIRE(initial_velocity > escape_velocity); + + std::vector distances; + std::vector velocities; + std::vector energies; + + double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; + int step_count = 0; + while (sim->time < max_time) { + if (step_count % 1000 == 0) { + double current_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); + double current_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double current_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); + double current_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], + &sim->bodies[SUN_INDEX]); + double current_total = current_kinetic + current_potential; + + distances.push_back(current_distance); + velocities.push_back(current_velocity); + energies.push_back(current_total); + } + + update_simulation(sim); + step_count++; + } + + double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); + double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + + double final_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); + double final_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], + &sim->bodies[SUN_INDEX]); + double final_total_energy = final_kinetic + final_potential; + + INFO("Final distance: " << final_distance / AU << " AU"); + INFO("Final velocity: " << final_velocity / 1000.0 << " km/s"); + INFO("Final kinetic energy: " << final_kinetic); + INFO("Final potential energy: " << final_potential); + INFO("Final total energy: " << final_total_energy); + + REQUIRE(final_distance > initial_distance); + + REQUIRE(final_velocity < initial_velocity); + + REQUIRE(final_velocity > 10000.0); + + double energy_drift = fabs(final_total_energy - initial_total_energy); + double avg_kinetic_energy = (initial_kinetic + final_kinetic) / 2.0; + double energy_drift_percent = (energy_drift / avg_kinetic_energy) * 100.0; + + INFO("Energy drift: " << energy_drift << " J"); + INFO("Energy drift percent: " << energy_drift_percent << "%"); + + REQUIRE(energy_drift_percent < 1.0); + + int velocity_decreases = 0; + for (size_t i = 1; i < velocities.size(); i++) { + if (velocities[i] < velocities[i-1]) { + velocity_decreases++; + } + } + + INFO("Velocity decreases: " << velocity_decreases << " / " << (velocities.size() - 1)); + + REQUIRE(velocity_decreases > static_cast(velocities.size()) / 2); + + destroy_simulation(sim); +} + +TEST_CASE("Hyperbolic orbit initial conditions", "[hyperbolic][initial]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/configs/hyperbolic_comet.toml")); + + const int COMET_INDEX = 1; + const int SUN_INDEX = 0; + CelestialBody* comet = &sim->bodies[COMET_INDEX]; + CelestialBody* sun = &sim->bodies[SUN_INDEX]; + + double distance = vec3_magnitude(vec3_sub(comet->position, sun->position)); + double velocity = vec3_magnitude(comet->velocity); + + double escape_velocity = sqrt(2.0 * G * sun->mass / distance); + double circular_velocity = sqrt(G * sun->mass / distance); + + INFO("Distance: " << distance / 1.496e11 << " AU"); + INFO("Actual velocity: " << velocity / 1000.0 << " km/s"); + INFO("Escape velocity: " << escape_velocity / 1000.0 << " km/s"); + INFO("Circular velocity: " << circular_velocity / 1000.0 << " km/s"); + + double velocity_excess = (velocity - escape_velocity) / escape_velocity; + INFO("Velocity excess over escape velocity: " << velocity_excess * 100.0 << "%"); + + REQUIRE(velocity > escape_velocity); + + INFO("Eccentricity: " << comet->eccentricity); + + REQUIRE(comet->eccentricity > 1.0); + REQUIRE(fabs(comet->eccentricity - 1.5) < 0.01); + + INFO("Semi-major axis: " << comet->semi_major_axis / 1.496e11 << " AU"); + + REQUIRE(comet->semi_major_axis < 0.0); + + double initial_kinetic = calculate_kinetic_energy(comet); + double initial_potential = calculate_potential_energy_pair(comet, sun); + double total_energy = initial_kinetic + initial_potential; + + INFO("Total energy: " << total_energy); + + REQUIRE(total_energy > 0.0); + + Vec3 r_vec = vec3_sub(comet->position, sun->position); + double r = vec3_magnitude(r_vec); + double a = comet->semi_major_axis; + + double expected_v_squared = G * sun->mass * (2.0 / r - 1.0 / a); + double expected_velocity = sqrt(expected_v_squared); + + double velocity_error = fabs(velocity - expected_velocity) / expected_velocity; + INFO("Velocity error from vis-viva: " << velocity_error * 100.0 << "%"); + + REQUIRE(velocity_error < 0.001); + + destroy_simulation(sim); +} + +TEST_CASE("Hyperbolic orbit asymptotic velocity", "[hyperbolic][asymptotic]") { + const double TIME_STEP = 60.0; + const double DAYS_TO_SIMULATE = 1000.0; + const double SECONDS_PER_DAY = 86400.0; + const double AU = 1.496e11; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/configs/hyperbolic_comet.toml")); + + const int COMET_INDEX = 1; + const int SUN_INDEX = 0; + + CelestialBody* sun = &sim->bodies[SUN_INDEX]; + double a = sim->bodies[COMET_INDEX].semi_major_axis; + + double expected_v_infinity = sqrt(2.0 * G * sun->mass / fabs(a)); + + INFO("Semi-major axis: " << a / 1.496e11 << " AU"); + INFO("Expected asymptotic velocity: " << expected_v_infinity / 1000.0 << " km/s"); + + double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; + int step_count = 0; + while (sim->time < max_time) { + update_simulation(sim); + step_count++; + + if (step_count % 10000 == 0) { + double distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); + if (distance > 20.0 * AU) { + INFO("Stopping simulation at distance: " << distance / AU << " AU"); + break; + } + } + } + + double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); + double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + + INFO("Final distance: " << final_distance / AU << " AU"); + INFO("Final velocity: " << final_velocity / 1000.0 << " km/s"); + + REQUIRE(final_distance > 18.0 * AU); + + double v_inf_error = fabs(final_velocity - expected_v_infinity) / expected_v_infinity; + INFO("Asymptotic velocity error: " << v_inf_error * 100.0 << "%"); + + REQUIRE(v_inf_error < 0.30); + + destroy_simulation(sim); +}