claudes_game/tests/informational/test_time_step_stability.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 <cstdio>

struct TestBody {
    const char* name;
    int body_index;
    int parent_index;
    double expected_period_days;
};

const double SECONDS_PER_DAY = 86400.0;
const double MIN_DT = 30.0;
const double MAX_DT = 600.0;
const double ENERGY_TOLERANCE = 1.0;
const int NUM_ORBITS = 100;

double calculate_orbital_period(CelestialBody* body, CelestialBody* parent) {
    Vec3 relative_pos = vec3_sub(body->global_position, parent->global_position);
    double r = vec3_magnitude(relative_pos);
    Vec3 relative_vel = vec3_sub(body->global_velocity, parent->global_velocity);
    double v = vec3_magnitude(relative_vel);
    
    double specific_energy = (v * v) / 2.0 - G * parent->mass / r;
    double semi_major_axis = -G * parent->mass / (2.0 * specific_energy);
    double period_seconds = 2.0 * M_PI * sqrt(pow(semi_major_axis, 3.0) / (G * parent->mass));
    
    return period_seconds;
}

bool is_dt_stable(SimulationState* sim, const TestBody& test_body, double dt, int num_orbits) {
    SimulationState* test_sim = create_simulation(sim->max_bodies, sim->max_craft, sim->max_maneuvers, dt);
    
    REQUIRE(load_system_config(test_sim, "tests/informational/test_time_step_stability.toml"));
    
    int body_index = test_body.body_index;
    int parent_index = test_body.parent_index;
    
    double initial_energy = calculate_system_total_energy(test_sim);
    
    Vec3 initial_pos_relative = vec3_sub(
        test_sim->bodies[body_index].global_position,
        test_sim->bodies[parent_index].global_position
    );
    double initial_distance = vec3_magnitude(initial_pos_relative);
    
    double period = calculate_orbital_period(&test_sim->bodies[body_index], &test_sim->bodies[parent_index]);
    double max_time = period * num_orbits;
    
    bool completed = true;
    while (test_sim->time < max_time) {
        update_simulation(test_sim);
        
        if (test_sim->bodies[body_index].parent_index != parent_index) {
            completed = false;
            break;
        }
    }
    
    double final_energy = calculate_system_total_energy(test_sim);
    double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0;
    
    Vec3 final_pos_relative = vec3_sub(
        test_sim->bodies[body_index].global_position,
        test_sim->bodies[parent_index].global_position
    );
    double final_distance = vec3_magnitude(final_pos_relative);
    double distance_drift_percent = fabs((final_distance - initial_distance) / initial_distance) * 100.0;
    
    bool stable = completed && (energy_drift_percent < ENERGY_TOLERANCE) && (distance_drift_percent < 5.0);
    
    destroy_simulation(test_sim);
    
    return stable;
}

double find_max_stable_dt(SimulationState* sim, const TestBody& test_body) {
    double low = MIN_DT;
    double high = MAX_DT;
    double max_stable = low;
    
    printf("Testing %s (period ~%.2f days):\n", test_body.name, test_body.expected_period_days);
    
    for (int iter = 0; iter < 10; iter++) {
        double mid = (low + high) / 2.0;
        bool stable = is_dt_stable(sim, test_body, mid, NUM_ORBITS);
        
        if (stable) {
            max_stable = mid;
            low = mid;
            printf("  dt=%.0fs: STABLE\n", mid);
        } else {
            high = mid;
            printf("  dt=%.0fs: UNSTABLE\n", mid);
        }
        
        if (high - low < 5.0) break;
    }
    
    printf("  Maximum stable dt: %.0f seconds\n\n", max_stable);
    return max_stable;
}

TEST_CASE("Time step stability - Mercury orbiter (MESSENGER-like)", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);
    
    TestBody mercury_orbiter = {"Mercury_Orbiter", 1, 0, 0.5};
    double max_dt = find_max_stable_dt(sim, mercury_orbiter);
    
    INFO("Mercury orbiter maximum stable dt: " << max_dt << " seconds");
    
    REQUIRE(max_dt >= MIN_DT);
    
    destroy_simulation(sim);
}

TEST_CASE("Time step stability - Io (Jupiter's moon)", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);
    
    TestBody io = {"Io", 3, 2, 1.77};
    double max_dt = find_max_stable_dt(sim, io);
    
    INFO("Io maximum stable dt: " << max_dt << " seconds");
    
    REQUIRE(max_dt >= MIN_DT);
    
    destroy_simulation(sim);
}

TEST_CASE("Time step stability - Moon (Earth's moon)", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);
    
    TestBody moon = {"Moon", 5, 4, 27.3};
    double max_dt = find_max_stable_dt(sim, moon);
    
    INFO("Moon maximum stable dt: " << max_dt << " seconds");
    
    REQUIRE(max_dt >= MIN_DT);
    
    destroy_simulation(sim);
}

TEST_CASE("Find minimum stable time step across all bodies", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);
    
    TestBody bodies[] = {
        {"Mercury_Orbiter", 1, 0, 0.5},
        {"Io", 3, 2, 1.77},
        {"Moon", 5, 4, 27.3}
    };
    
    double max_dt = MAX_DT;
    
    printf("\n=== Finding minimum stable dt across all bodies ===\n\n");
    
    for (int i = 0; i < 3; i++) {
        double body_max_dt = find_max_stable_dt(sim, bodies[i]);
        if (body_max_dt < max_dt) {
            max_dt = body_max_dt;
        }
    }
    
    printf("\n=== RESULTS ===\n");
    printf("Minimum stable time step: %.0f seconds\n", max_dt);
    printf("Recommended safe time step: %.0f seconds (%.0fx safety margin)\n", max_dt * 0.7, 1.0/0.7);
    
    INFO("Minimum stable dt: " << max_dt << " seconds");
    
    REQUIRE(max_dt >= MIN_DT);
    
    destroy_simulation(sim);
}

TEST_CASE("Verify current default dt (60s) stability", "[timestep][stability]") {
    const double DT = 60.0;
    const int NUM_ORBITS = 10;
    
    SimulationState* sim = create_simulation(10, 0, 0, DT);
    REQUIRE(load_system_config(sim, "tests/informational/test_time_step_stability.toml"));
    
    struct BodyTest {
        int body_index;
        int parent_index;
        const char* name;
    };
    
    BodyTest tests[] = {
        {1, 0, "Mercury_Orbiter"},
        {3, 2, "Io"},
        {5, 4, "Moon"}
    };
    
    for (int t = 0; t < 3; t++) {
        int body_index = tests[t].body_index;
        int parent_index = tests[t].parent_index;
        const char* name = tests[t].name;
        
        double period = calculate_orbital_period(&sim->bodies[body_index], &sim->bodies[parent_index]);
        double max_time = period * NUM_ORBITS;
        
        double initial_energy = calculate_system_total_energy(sim);
        
        INFO("Testing " << name << " with dt=" << DT << "s for " << NUM_ORBITS << " orbits");
        
        bool completed = true;
        while (sim->time < max_time) {
            update_simulation(sim);
            
            if (sim->bodies[body_index].parent_index != parent_index) {
                completed = false;
                break;
            }
        }
        
        double final_energy = calculate_system_total_energy(sim);
        double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0;
        
        INFO(name << " completed: " << (completed ? "yes" : "no"));
        INFO(name << " energy drift: " << energy_drift_percent << "%");
        
        REQUIRE(completed);
        REQUIRE(energy_drift_percent < ENERGY_TOLERANCE);
    }
    
    destroy_simulation(sim);
}