diff --git a/src/simulation.cpp b/src/simulation.cpp index 737da51..a9da247 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -119,18 +119,21 @@ static void compute_orbital_velocity_from_vis_viva(CelestialBody* body, CelestialBody* parent) { Vec3 r = vec3_sub(body->position, parent->position); double distance = vec3_magnitude(r); - // double e = body->eccentricity; // FIXME: unused variable + double e = body->eccentricity; double a = body->semi_major_axis; - double v_squared = G * parent->mass * (2.0 / distance - 1.0 / a); - // FIXME: need error handling here - // if v_squared is negative, the inital body parameters are bad + double v_squared; + + if (fabs(e - 1.0) < 0.0001) { + v_squared = 2.0 * G * parent->mass / distance; + } else { + v_squared = G * parent->mass * (2.0 / distance - 1.0 / a); + } + assert(v_squared >= 0); double speed = (double) sqrt(v_squared); Vec3 z_axis = {0.0, 0.0, 1.0}; Vec3 vel_dir = vec3_cross(r, z_axis); - // NOTE: I suppose this is for the case of a high inclination orbit? - // FIXME: make a test to see if this is necessary if (vec3_magnitude(vel_dir) < 0.01) { Vec3 x_axis = {1.0, 0.0, 0.0}; vel_dir = vec3_cross(r, x_axis); @@ -138,7 +141,6 @@ static void compute_orbital_velocity_from_vis_viva(CelestialBody* body, vel_dir = vec3_normalize(vel_dir); body->velocity = vec3_scale(vel_dir, speed); - // FIXME: this is wrong, we should be simulating in parent coord frame body->velocity = vec3_add(body->velocity, parent->velocity); } diff --git a/tests/configs/parabolic_comet.toml b/tests/configs/parabolic_comet.toml new file mode 100644 index 0000000..7eb08e1 --- /dev/null +++ b/tests/configs/parabolic_comet.toml @@ -0,0 +1,23 @@ +# Test Configuration: Sun + Parabolic Comet +# Comet with parabolic orbit (eccentricity = 1.0) +# Escape trajectory - total energy = 0 + +[[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 = "ParabolicComet" +mass = 1.0e14 +radius = 5.0e3 +position = { x = 1.496e11, y = 0.0, z = 0.0 } +parent_index = 0 +color = { r = 0.7, g = 0.8, b = 0.9 } +eccentricity = 1.0 +semi_major_axis = 1.0e30 diff --git a/tests/test_parabolic_orbit.cpp b/tests/test_parabolic_orbit.cpp new file mode 100644 index 0000000..14e9819 --- /dev/null +++ b/tests/test_parabolic_orbit.cpp @@ -0,0 +1,137 @@ +#include +#include "../src/physics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include +#include + +TEST_CASE("Parabolic orbit - energy and escape trajectory", "[parabolic][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/parabolic_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; + + INFO("Initial distance: " << initial_distance / AU << " AU"); + INFO("Initial velocity: " << vec3_magnitude(sim->bodies[COMET_INDEX].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 >= -1e25); + + 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); + + 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("Parabolic orbit initial conditions", "[parabolic][initial]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/configs/parabolic_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_error = fabs(velocity - escape_velocity) / escape_velocity; + INFO("Velocity error from escape velocity: " << velocity_error * 100.0 << "%"); + + REQUIRE(velocity_error < 0.001); + + INFO("Eccentricity: " << comet->eccentricity); + + REQUIRE(fabs(comet->eccentricity - 1.0) < 0.0001); + + destroy_simulation(sim); +}