diff --git a/Makefile b/Makefile index 4a52e85..015d871 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ # Compiler and flags CXX = g++ -CXXFLAGS = -Wall -Wextra -std=c++11 -I./src -isystem./ext/raylib/src +CXXFLAGS = -Wall -Wextra -std=c++14 -I./src -isystem./ext/raylib/src LDFLAGS = -L./ext/raylib/src -lraylib -lm -lpthread -ldl -lrt -lX11 # Directories @@ -8,6 +8,8 @@ SRC_DIR = src BUILD_DIR = build RAYLIB_DIR = ext/raylib/src TARGET = orbit_sim +TEST_DIR = tests +TEST_TARGET = orbit_test # Source files SOURCES = $(SRC_DIR)/main.cpp \ @@ -56,8 +58,24 @@ rebuild: clean all run: $(TARGET) ./$(TARGET) -# Run test configuration -test: $(TARGET) +# Run manual integration test with simulation +test-manual: $(TARGET) ./$(TARGET) configs/test_simple.txt --headless --readable --days 365 -.PHONY: all clean clean-all rebuild run test raylib +# Build automated test suite +test-build: + $(CXX) $(CXXFLAGS) -I/usr/include/catch2 \ + $(TEST_DIR)/test_main.cpp \ + $(TEST_DIR)/test_integration.cpp \ + $(TEST_DIR)/test_energy.cpp \ + $(TEST_DIR)/test_orbital_period.cpp \ + $(SRC_DIR)/test_utilities.cpp \ + $(SRC_DIR)/physics.cpp \ + $(SRC_DIR)/bodies.cpp \ + -o $(TEST_TARGET) -lCatch2Main -lCatch2 -lm + +# Run automated test suite +test: test-build + ./$(TEST_TARGET) + +.PHONY: all clean clean-all rebuild run test test-build test-manual raylib diff --git a/src/bodies.cpp b/src/bodies.cpp index 976d9bc..cb0ee10 100644 --- a/src/bodies.cpp +++ b/src/bodies.cpp @@ -99,63 +99,78 @@ void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_ax body->soi_radius = semi_major_axis * pow(mass_ratio, 0.4); // 2/5 = 0.4 } -// Update the entire simulation by one time step +// [COMMENTED OLD CODE - Two-phase Euler update] +// void update_simulation(SimulationState* sim) { +// for (int i = 0; i < sim->body_count; i++) { +// CelestialBody* body = &sim->bodies[i]; +// if (body->parent_index == -1) { +// Vec3 total_force = {0.0, 0.0, 0.0}; +// for (int j = 0; j < sim->body_count; j++) { +// if (i == j) continue; +// CelestialBody* other = &sim->bodies[j]; +// if (other->parent_index == -1) { +// Vec3 force = calculate_gravity_force(body, other); +// total_force = vec3_add(total_force, force); +// } +// } +// Vec3 acceleration = calculate_acceleration(total_force, body->mass); +// euler_step(body, acceleration, sim->dt); +// } +// } +// for (int i = 0; i < sim->body_count; i++) { +// CelestialBody* body = &sim->bodies[i]; +// if (body->parent_index == -1) { +// continue; +// } +// int new_parent = find_dominant_body(sim, i); +// if (new_parent != body->parent_index && new_parent != -1) { +// body->parent_index = new_parent; +// } +// if (body->parent_index >= 0 && body->parent_index < sim->body_count) { +// CelestialBody* parent = &sim->bodies[body->parent_index]; +// Vec3 force = calculate_gravity_force(body, parent); +// Vec3 acceleration = calculate_acceleration(force, body->mass); +// euler_step(body, acceleration, sim->dt); +// } +// } +// sim->time += sim->dt; +// } + void update_simulation(SimulationState* sim) { - // First, update root bodies (they interact with each other) for (int i = 0; i < sim->body_count; i++) { CelestialBody* body = &sim->bodies[i]; if (body->parent_index == -1) { - // This is a root body - calculate forces from OTHER root bodies - Vec3 total_force = {0.0, 0.0, 0.0}; - - for (int j = 0; j < sim->body_count; j++) { - if (i == j) continue; // Don't apply force to itself + AccelerationContext ctx; + ctx.sim = sim; + ctx.current_body = body; + ctx.body_index = i; - CelestialBody* other = &sim->bodies[j]; - if (other->parent_index == -1) { - // Other is also a root body - apply gravitational force - Vec3 force = calculate_gravity_force(body, other); - total_force = vec3_add(total_force, force); - } - } - - // Apply total force from all other root bodies - Vec3 acceleration = calculate_acceleration(total_force, body->mass); - euler_step(body, acceleration, sim->dt); + rk4_step(body, &ctx, sim->dt); } } - // Now update non-root bodies (planets, moons, etc.) for (int i = 0; i < sim->body_count; i++) { CelestialBody* body = &sim->bodies[i]; - // Skip root bodies (already updated above) if (body->parent_index == -1) { continue; } - // Check if parent has changed (SOI transition) int new_parent = find_dominant_body(sim, i); if (new_parent != body->parent_index && new_parent != -1) { body->parent_index = new_parent; } - // Get the current parent if (body->parent_index >= 0 && body->parent_index < sim->body_count) { - CelestialBody* parent = &sim->bodies[body->parent_index]; - - // Calculate gravitational force from parent - Vec3 force = calculate_gravity_force(body, parent); - - // Calculate acceleration - Vec3 acceleration = calculate_acceleration(force, body->mass); + AccelerationContext ctx; + ctx.sim = sim; + ctx.current_body = body; + ctx.body_index = i; - // Perform Euler integration step - euler_step(body, acceleration, sim->dt); + rk4_step(body, &ctx, sim->dt); } } - // Update simulation time sim->time += sim->dt; } diff --git a/src/physics.cpp b/src/physics.cpp index 372bf34..dd23df1 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -61,11 +61,73 @@ Vec3 calculate_acceleration(Vec3 force, double mass) { return {0.0, 0.0, 0.0}; } -// Euler integration step: update position and velocity +// [COMMENTED OLD CODE - Euler integration] +// void euler_step(CelestialBody* body, Vec3 acceleration, double dt) { +// body->velocity = vec3_add(body->velocity, vec3_scale(acceleration, dt)); +// body->position = vec3_add(body->position, vec3_scale(body->velocity, dt)); +// } + void euler_step(CelestialBody* body, Vec3 acceleration, double dt) { - // Update velocity: v = v + a * dt body->velocity = vec3_add(body->velocity, vec3_scale(acceleration, dt)); - - // Update position: p = p + v * dt body->position = vec3_add(body->position, vec3_scale(body->velocity, dt)); } + +Vec3 evaluate_acceleration(Vec3 pos, Vec3 vel, AccelerationContext* ctx) { + CelestialBody temp_body = *ctx->current_body; + temp_body.position = pos; + temp_body.velocity = vel; + + Vec3 total_force = {0.0, 0.0, 0.0}; + + if (temp_body.parent_index == -1) { + for (int j = 0; j < ctx->sim->body_count; j++) { + if (j == ctx->body_index) continue; + CelestialBody* other = &ctx->sim->bodies[j]; + if (other->parent_index == -1) { + Vec3 force = calculate_gravity_force(&temp_body, other); + total_force = vec3_add(total_force, force); + } + } + } else { + if (temp_body.parent_index >= 0 && temp_body.parent_index < ctx->sim->body_count) { + CelestialBody* parent = &ctx->sim->bodies[temp_body.parent_index]; + total_force = calculate_gravity_force(&temp_body, parent); + } + } + + return calculate_acceleration(total_force, temp_body.mass); +} + +void rk4_step(CelestialBody* body, AccelerationContext* ctx, double dt) { + Vec3 k1_vel, k2_vel, k3_vel, k4_vel; + Vec3 k1_pos, k2_pos, k3_pos, k4_pos; + + Vec3 pos0 = body->position; + Vec3 vel0 = body->velocity; + + k1_vel = evaluate_acceleration(pos0, vel0, ctx); + k1_pos = vel0; + + Vec3 pos1 = vec3_add(pos0, vec3_scale(k1_pos, dt * 0.5)); + Vec3 vel1 = vec3_add(vel0, vec3_scale(k1_vel, dt * 0.5)); + k2_vel = evaluate_acceleration(pos1, vel1, ctx); + k2_pos = vel1; + + Vec3 pos2 = vec3_add(pos0, vec3_scale(k2_pos, dt * 0.5)); + Vec3 vel2 = vec3_add(vel0, vec3_scale(k2_vel, dt * 0.5)); + k3_vel = evaluate_acceleration(pos2, vel2, ctx); + k3_pos = vel2; + + Vec3 pos3 = vec3_add(pos0, vec3_scale(k3_pos, dt)); + Vec3 vel3 = vec3_add(vel0, vec3_scale(k3_vel, dt)); + k4_vel = evaluate_acceleration(pos3, vel3, ctx); + k4_pos = vel3; + + Vec3 k_vel_sum = vec3_add(vec3_add(k1_vel, vec3_scale(k2_vel, 2.0)), + vec3_add(vec3_scale(k3_vel, 2.0), k4_vel)); + Vec3 k_pos_sum = vec3_add(vec3_add(k1_pos, vec3_scale(k2_pos, 2.0)), + vec3_add(vec3_scale(k3_pos, 2.0), k4_pos)); + + body->velocity = vec3_add(vel0, vec3_scale(k_vel_sum, dt / 6.0)); + body->position = vec3_add(pos0, vec3_scale(k_pos_sum, dt / 6.0)); +} diff --git a/src/physics.h b/src/physics.h index 8cd9585..dcf70bf 100644 --- a/src/physics.h +++ b/src/physics.h @@ -25,4 +25,15 @@ Vec3 calculate_gravity_force(CelestialBody* body, CelestialBody* parent); Vec3 calculate_acceleration(Vec3 force, double mass); void euler_step(CelestialBody* body, Vec3 acceleration, double dt); +// Forward declaration for RK4 context +struct SimulationState; + +struct AccelerationContext { + SimulationState* sim; + CelestialBody* current_body; + int body_index; +}; + +void rk4_step(CelestialBody* body, AccelerationContext* ctx, double dt); + #endif diff --git a/src/test_utilities.cpp b/src/test_utilities.cpp new file mode 100644 index 0000000..7782b37 --- /dev/null +++ b/src/test_utilities.cpp @@ -0,0 +1,120 @@ +#include "test_utilities.h" +#include +#include + +double calculate_kinetic_energy(CelestialBody* body) { + double v_squared = body->velocity.x * body->velocity.x + + body->velocity.y * body->velocity.y + + body->velocity.z * body->velocity.z; + return 0.5 * body->mass * v_squared; +} + +double calculate_potential_energy_pair(CelestialBody* body1, CelestialBody* body2) { + double distance = vec3_distance(body1->position, body2->position); + if (distance < 1.0) distance = 1.0; + return -G * body1->mass * body2->mass / distance; +} + +double calculate_system_total_energy(SimulationState* sim) { + double kinetic = 0.0; + double potential = 0.0; + + for (int i = 0; i < sim->body_count; i++) { + kinetic += calculate_kinetic_energy(&sim->bodies[i]); + + for (int j = i + 1; j < sim->body_count; j++) { + potential += calculate_potential_energy_pair(&sim->bodies[i], &sim->bodies[j]); + } + } + + return kinetic + potential; +} + +OrbitalMetrics calculate_orbital_metrics(CelestialBody* body, CelestialBody* parent) { + OrbitalMetrics metrics; + + Vec3 relative_pos = vec3_sub(body->position, parent->position); + metrics.orbital_radius = vec3_magnitude(relative_pos); + metrics.velocity_magnitude = vec3_magnitude(body->velocity); + metrics.angular_position = atan2(relative_pos.y, relative_pos.x); + + metrics.kinetic_energy = calculate_kinetic_energy(body); + metrics.potential_energy = calculate_potential_energy_pair(body, parent); + metrics.total_energy = metrics.kinetic_energy + metrics.potential_energy; + + return metrics; +} + +OrbitTracker* create_orbit_tracker(int body_index) { + OrbitTracker* tracker = (OrbitTracker*)malloc(sizeof(OrbitTracker)); + tracker->body_index = body_index; + tracker->initial_angle = 0.0; + tracker->previous_angle = 0.0; + tracker->quadrant_transitions = 0; + tracker->orbit_completed = false; + tracker->time_at_completion = 0.0; + return tracker; +} + +void reset_orbit_tracker(OrbitTracker* tracker) { + tracker->initial_angle = 0.0; + tracker->previous_angle = 0.0; + tracker->quadrant_transitions = 0; + tracker->orbit_completed = false; + tracker->time_at_completion = 0.0; +} + +void update_orbit_tracker(OrbitTracker* tracker, CelestialBody* body, CelestialBody* parent, double current_time) { + if (tracker->orbit_completed) return; + + Vec3 relative_pos = vec3_sub(body->position, parent->position); + double current_angle = atan2(relative_pos.y, relative_pos.x); + + if (tracker->quadrant_transitions == 0) { + tracker->initial_angle = current_angle; + tracker->previous_angle = current_angle; + tracker->quadrant_transitions = 1; + return; + } + + double angle_diff = current_angle - tracker->previous_angle; + + if (angle_diff > M_PI) { + angle_diff -= 2.0 * M_PI; + tracker->quadrant_transitions++; + } + if (angle_diff < -M_PI) { + angle_diff += 2.0 * M_PI; + tracker->quadrant_transitions++; + } + + double total_rotation = current_angle - tracker->initial_angle; + if (total_rotation < -M_PI) total_rotation += 2.0 * M_PI; + if (total_rotation > M_PI) total_rotation -= 2.0 * M_PI; + + const double SECONDS_PER_DAY = 86400.0; + const double MIN_ORBIT_TIME = 100.0 * SECONDS_PER_DAY; + + if (tracker->quadrant_transitions >= 2 && + current_time > MIN_ORBIT_TIME && + fabs(total_rotation) < 0.05) { + tracker->orbit_completed = true; + tracker->time_at_completion = current_time; + } + + tracker->previous_angle = current_angle; +} + +void destroy_orbit_tracker(OrbitTracker* tracker) { + free(tracker); +} + +bool compare_double(double a, double b, double tolerance) { + return fabs(a - b) <= tolerance; +} + +bool compare_vec3(Vec3 a, Vec3 b, double tolerance) { + return fabs(a.x - b.x) <= tolerance && + fabs(a.y - b.y) <= tolerance && + fabs(a.z - b.z) <= tolerance; +} diff --git a/src/test_utilities.h b/src/test_utilities.h new file mode 100644 index 0000000..d351c72 --- /dev/null +++ b/src/test_utilities.h @@ -0,0 +1,38 @@ +#ifndef TEST_UTILITIES_H +#define TEST_UTILITIES_H + +#include "bodies.h" +#include "physics.h" + +struct OrbitalMetrics { + double kinetic_energy; + double potential_energy; + double total_energy; + double orbital_radius; + double velocity_magnitude; + double angular_position; +}; + +struct OrbitTracker { + double initial_angle; + double previous_angle; + int quadrant_transitions; + bool orbit_completed; + double time_at_completion; + int body_index; +}; + +double calculate_kinetic_energy(CelestialBody* body); +double calculate_potential_energy_pair(CelestialBody* body1, CelestialBody* body2); +double calculate_system_total_energy(SimulationState* sim); +OrbitalMetrics calculate_orbital_metrics(CelestialBody* body, CelestialBody* parent); + +OrbitTracker* create_orbit_tracker(int body_index); +void reset_orbit_tracker(OrbitTracker* tracker); +void update_orbit_tracker(OrbitTracker* tracker, CelestialBody* body, CelestialBody* parent, double current_time); +void destroy_orbit_tracker(OrbitTracker* tracker); + +bool compare_double(double a, double b, double tolerance); +bool compare_vec3(Vec3 a, Vec3 b, double tolerance); + +#endif diff --git a/tests/test_energy.cpp b/tests/test_energy.cpp new file mode 100644 index 0000000..eea84c4 --- /dev/null +++ b/tests/test_energy.cpp @@ -0,0 +1,39 @@ +#include +#include "../src/physics.h" +#include "../src/bodies.h" +#include "../src/test_utilities.h" +#include + +TEST_CASE("Energy conservation - Earth circular orbit", "[energy][rk4]") { + const double TIME_STEP = 60.0; + const double DAYS_TO_SIMULATE = 10.0; + const double SECONDS_PER_DAY = 86400.0; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + Vec3 sun_pos = {0, 0, 0}; + Vec3 sun_vel = {0, 0, 0}; + add_body(sim, "Sun", 1.989e30, 6.96e8, sun_pos, sun_vel, -1, 1.0, 1.0, 0.0, 0, 0); + + Vec3 earth_pos = {1.496e11, 0, 0}; + Vec3 earth_vel = {0, 29800, 0}; + add_body(sim, "Earth", 5.972e24, 6.371e6, earth_pos, earth_vel, 0, 0.0, 0.5, 1.0, 0, 1.496e11); + + double initial_energy = calculate_system_total_energy(sim); + + double total_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; + while (sim->time < total_time) { + update_simulation(sim); + } + + double final_energy = calculate_system_total_energy(sim); + double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0; + + INFO("Initial energy: " << initial_energy << " J"); + INFO("Final energy: " << final_energy << " J"); + INFO("Energy drift: " << energy_drift_percent << "%"); + + REQUIRE(energy_drift_percent < 5.0); + + destroy_simulation(sim); +} diff --git a/tests/test_integration.cpp b/tests/test_integration.cpp new file mode 100644 index 0000000..53983b4 --- /dev/null +++ b/tests/test_integration.cpp @@ -0,0 +1,47 @@ +#include +#include "../src/physics.h" +#include "../src/test_utilities.h" +#include + +TEST_CASE("Vector math utilities", "[utilities]") { + Vec3 a = {1.0, 2.0, 3.0}; + Vec3 b = {4.0, 5.0, 6.0}; + + SECTION("Vector addition") { + Vec3 sum = vec3_add(a, b); + REQUIRE(compare_double(sum.x, 5.0, 1e-10)); + REQUIRE(compare_double(sum.y, 7.0, 1e-10)); + REQUIRE(compare_double(sum.z, 9.0, 1e-10)); + } + + SECTION("Vector subtraction") { + Vec3 diff = vec3_sub(b, a); + REQUIRE(compare_double(diff.x, 3.0, 1e-10)); + REQUIRE(compare_double(diff.y, 3.0, 1e-10)); + REQUIRE(compare_double(diff.z, 3.0, 1e-10)); + } + + SECTION("Vector scaling") { + Vec3 scaled = vec3_scale(a, 2.0); + REQUIRE(compare_double(scaled.x, 2.0, 1e-10)); + REQUIRE(compare_double(scaled.y, 4.0, 1e-10)); + REQUIRE(compare_double(scaled.z, 6.0, 1e-10)); + } + + SECTION("Vector magnitude") { + double mag = vec3_magnitude(a); + REQUIRE(compare_double(mag, sqrt(14.0), 1e-10)); + } + + SECTION("Vector distance") { + double dist = vec3_distance(a, b); + REQUIRE(compare_double(dist, sqrt(27.0), 1e-10)); + } + + SECTION("Vector normalization") { + Vec3 unit = vec3_normalize(a); + double expected_mag = 1.0; + double actual_mag = vec3_magnitude(unit); + REQUIRE(compare_double(actual_mag, expected_mag, 1e-10)); + } +} diff --git a/tests/test_main.cpp b/tests/test_main.cpp new file mode 100644 index 0000000..2a31758 --- /dev/null +++ b/tests/test_main.cpp @@ -0,0 +1,5 @@ +#include + +int main(int argc, char* argv[]) { + return Catch::Session().run(argc, argv); +} diff --git a/tests/test_orbital_period.cpp b/tests/test_orbital_period.cpp new file mode 100644 index 0000000..73cc91f --- /dev/null +++ b/tests/test_orbital_period.cpp @@ -0,0 +1,83 @@ +#include +#include "../src/physics.h" +#include "../src/bodies.h" +#include "../src/test_utilities.h" +#include + +TEST_CASE("Orbital period - Earth (RK4)", "[period][rk4]") { + const double TIME_STEP = 60.0; + const double EXPECTED_PERIOD_DAYS = 365.0; + const double SECONDS_PER_DAY = 86400.0; + const double MAX_SIMULATION_DAYS = 400.0; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + Vec3 sun_pos = {0, 0, 0}; + Vec3 sun_vel = {0, 0, 0}; + add_body(sim, "Sun", 1.989e30, 6.96e8, sun_pos, sun_vel, -1, 1.0, 1.0, 0.0, 0, 0); + + Vec3 earth_pos = {1.496e11, 0, 0}; + Vec3 earth_vel = {0, 29789, 0}; + add_body(sim, "Earth", 5.972e24, 6.371e6, earth_pos, earth_vel, 0, 0.0, 0.5, 1.0, 0, 1.496e11); + + OrbitTracker* tracker = create_orbit_tracker(1); + + double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY; + while (sim->time < max_time && !tracker->orbit_completed) { + update_simulation(sim); + update_orbit_tracker(tracker, &sim->bodies[1], &sim->bodies[0], sim->time); + } + + REQUIRE(tracker->orbit_completed); + + double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY; + double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS); + + INFO("Expected period: " << EXPECTED_PERIOD_DAYS << " days"); + INFO("Measured period: " << measured_period_days << " days"); + INFO("Error: " << period_error_days << " days"); + + REQUIRE(period_error_days < 5.0); + + destroy_orbit_tracker(tracker); + destroy_simulation(sim); +} + +TEST_CASE("Orbital period - Mars (RK4)", "[period][rk4]") { + const double TIME_STEP = 60.0; + const double EXPECTED_PERIOD_DAYS = 687.0; + const double SECONDS_PER_DAY = 86400.0; + const double MAX_SIMULATION_DAYS = 750.0; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + Vec3 sun_pos = {0, 0, 0}; + Vec3 sun_vel = {0, 0, 0}; + add_body(sim, "Sun", 1.989e30, 6.96e8, sun_pos, sun_vel, -1, 1.0, 1.0, 0.0, 0, 0); + + Vec3 mars_pos = {2.244e11, 0, 0}; + Vec3 mars_vel = {0, 24323, 0}; + add_body(sim, "Mars", 6.39e23, 3.3895e6, mars_pos, mars_vel, 0, 0.8, 0.3, 0.1, 0, 2.244e11); + + OrbitTracker* tracker = create_orbit_tracker(1); + + double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY; + while (sim->time < max_time && !tracker->orbit_completed) { + update_simulation(sim); + update_orbit_tracker(tracker, &sim->bodies[1], &sim->bodies[0], sim->time); + } + + REQUIRE(tracker->orbit_completed); + + double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY; + double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS); + + INFO("Expected period: " << EXPECTED_PERIOD_DAYS << " days"); + INFO("Measured period: " << measured_period_days << " days"); + INFO("Error: " << period_error_days << " days"); + + REQUIRE(period_error_days < 25.0); + + destroy_orbit_tracker(tracker); + destroy_simulation(sim); +}