claudes_game/tests/test_hyperbolic_orbit.cpp

#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, 0, 0, TIME_STEP);

    REQUIRE(load_system_config(sim, "tests/test_hyperbolic_orbit.toml"));

    const int COMET_INDEX = 1;
    const int SUN_INDEX = 0;

    Vec3 initial_position = sim->bodies[COMET_INDEX].global_position;
    double initial_distance = vec3_magnitude(initial_position);
    double initial_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_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].global_position);
            double current_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_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].global_position);
    double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_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, 0, 0, TIME_STEP);

    REQUIRE(load_system_config(sim, "tests/test_hyperbolic_orbit.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->global_position, sun->global_position));
    double velocity = vec3_magnitude(comet->global_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->orbit.eccentricity);

    REQUIRE(comet->orbit.eccentricity > 1.0);
    REQUIRE(fabs(comet->orbit.eccentricity - 1.5) < 0.01);

    INFO("Semi-major axis: " << comet->orbit.semi_major_axis / 1.496e11 << " AU");

    REQUIRE(comet->orbit.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->global_position, sun->global_position);
    double r = vec3_magnitude(r_vec);
    double a = comet->orbit.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, 0, 0, TIME_STEP);

    REQUIRE(load_system_config(sim, "tests/test_hyperbolic_orbit.toml"));

    const int COMET_INDEX = 1;
    const int SUN_INDEX = 0;

    CelestialBody* sun = &sim->bodies[SUN_INDEX];
    double a = sim->bodies[COMET_INDEX].orbit.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].global_position);
            if (distance > 20.0 * AU) {
                INFO("Stopping simulation at distance: " << distance / AU << " AU");
                break;
            }
        }
    }

    double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].global_position);
    double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_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);
}