Add hyperbolic orbit test case and configuration

- Add hyperbolic_comet.toml with e=1.5, a=-1.496e11
- Add 3 test cases:
  * Energy and escape trajectory (300 days)
  * Initial conditions verification
  * Asymptotic velocity validation (1000 days)
- Tests validate: v > v_escape, E > 0, e > 1.0, a < 0
- Asymptotic test uses 30% tolerance at 18+ AU distance
- Existing hyperbolic rendering unchanged

Claude

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

tests/test_hyperbolic_orbit.cpp

@@ -0,0 +1,221 @@
+#include <catch2/catch_test_macros.hpp>
+#include "../src/physics.h"
+#include "../src/simulation.h"
+#include "../src/config_loader.h"
+#include "../src/test_utilities.h"
+#include <cmath>
+#include <vector>
+
+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<double> distances;
+    std::vector<double> velocities;
+    std::vector<double> 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<int>(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);
+}