You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
192 lines
8.2 KiB
192 lines
8.2 KiB
#include <catch2/catch_test_macros.hpp> |
|
#include <catch2/matchers/catch_matchers_floating_point.hpp> |
|
#include "../src/physics.h" |
|
#include "../src/simulation.h" |
|
#include "../src/orbital_objects.h" |
|
#include "../src/maneuver.h" |
|
#include "../src/config_loader.h" |
|
#include "../src/test_utilities.h" |
|
#include <cmath> |
|
|
|
using Catch::Matchers::WithinAbs; |
|
|
|
SCENARIO("Spacecraft loading and impulsive burn behavior", "[spacecraft][config][burn]") { |
|
const double TIME_STEP = 60.0; |
|
const double SECONDS_TO_SIMULATE = 3600.0; |
|
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP); |
|
|
|
REQUIRE(load_system_config(sim, "tests/test_maneuvers.toml")); |
|
|
|
Spacecraft* craft = &sim->spacecraft[0]; |
|
CelestialBody* earth = &sim->bodies[1]; |
|
|
|
const double initial_r = vec3_distance(craft->local_position, (Vec3){0,0,0}); |
|
const double initial_a = craft->orbit.semi_major_axis; |
|
const double initial_e = craft->orbit.eccentricity; |
|
|
|
auto simulate = [&]() { |
|
double sim_time = 0.0; |
|
while (sim_time < SECONDS_TO_SIMULATE) { |
|
update_simulation(sim); |
|
sim_time += TIME_STEP; |
|
} |
|
}; |
|
|
|
SECTION("spacecraft loads correctly") { |
|
REQUIRE_THAT(sim->craft_count, WithinAbs(1.0, 0.001)); |
|
REQUIRE(std::string(sim->spacecraft[0].name) == "LEO_Satellite"); |
|
REQUIRE_THAT(sim->spacecraft[0].parent_index, WithinAbs(1.0, 0.001)); |
|
} |
|
|
|
SECTION("prograde burn increases velocity and raises apoapsis") { |
|
double v_before = vec3_magnitude(craft->local_velocity); |
|
apply_impulsive_burn(craft, BURN_PROGRADE, 100.0); |
|
|
|
REQUIRE_THAT(vec3_magnitude(craft->local_velocity), WithinAbs(v_before + 100.0, 0.001)); |
|
|
|
simulate(); |
|
|
|
double final_r = vec3_distance(craft->local_position, (Vec3){0,0,0}); |
|
INFO("Initial r: " << initial_r << " m"); |
|
INFO("Final r: " << final_r << " m"); |
|
INFO("delta_r: " << (final_r - initial_r) << " m"); |
|
|
|
// Prograde burn raises apoapsis; after ~225° craft is past periapsis toward apoapsis |
|
REQUIRE_THAT(final_r, WithinAbs(7085656.0, 320000.0)); |
|
} |
|
|
|
SECTION("retrograde burn decreases velocity and lowers periapsis") { |
|
double v_before = vec3_magnitude(craft->local_velocity); |
|
apply_impulsive_burn(craft, BURN_RETROGRADE, 100.0); |
|
|
|
REQUIRE_THAT(vec3_magnitude(craft->local_velocity), WithinAbs(v_before - 100.0, 0.001)); |
|
|
|
simulate(); |
|
|
|
double final_r = vec3_distance(craft->local_position, (Vec3){0,0,0}); |
|
INFO("Initial r: " << initial_r << " m"); |
|
INFO("Final r: " << final_r << " m"); |
|
INFO("delta_r: " << (final_r - initial_r) << " m"); |
|
|
|
// Retrograde burn lowers periapsis; after ~240° craft is near periapsis |
|
REQUIRE_THAT(final_r, WithinAbs(6525686.0, 250000.0)); |
|
} |
|
|
|
SECTION("normal burn changes orbital plane (z displacement)") { |
|
double initial_z = craft->local_position.z; |
|
apply_impulsive_burn(craft, BURN_NORMAL, 500.0); |
|
|
|
simulate(); |
|
|
|
double z_change = fabs(craft->local_position.z - initial_z); |
|
INFO("Initial z: " << initial_z << " m"); |
|
INFO("Final z: " << craft->local_position.z << " m"); |
|
INFO("|z_change|: " << z_change << " m"); |
|
|
|
// Normal burn introduces inclination; after 3600s z displacement ~348km |
|
REQUIRE_THAT(z_change, WithinAbs(348678.0, 350000.0)); |
|
} |
|
|
|
SECTION("custom burn applies arbitrary delta-v vector") { |
|
Vec3 initial_vel = craft->local_velocity; |
|
Vec3 delta_v = {10.0, 20.0, 30.0}; |
|
|
|
apply_custom_burn(craft, delta_v); |
|
|
|
REQUIRE_THAT(craft->local_velocity.x, WithinAbs(initial_vel.x + 10.0, 0.001)); |
|
REQUIRE_THAT(craft->local_velocity.y, WithinAbs(initial_vel.y + 20.0, 0.001)); |
|
REQUIRE_THAT(craft->local_velocity.z, WithinAbs(initial_vel.z + 30.0, 0.001)); |
|
} |
|
|
|
SECTION("propagation stability over 1 day") { |
|
const double total_time = 86400.0; |
|
double sim_time = 0.0; |
|
while (sim_time < total_time) { |
|
update_simulation(sim); |
|
sim_time += TIME_STEP; |
|
} |
|
|
|
double final_r = vec3_distance(craft->local_position, (Vec3){0,0,0}); |
|
double final_a = craft->orbit.semi_major_axis; |
|
double final_e = craft->orbit.eccentricity; |
|
|
|
double distance_drift_pct = fabs((final_r - initial_r) / initial_r) * 100.0; |
|
double a_drift_pct = fabs((final_a - initial_a) / initial_a) * 100.0; |
|
|
|
INFO("Initial r: " << initial_r << " m"); |
|
INFO("Final r: " << final_r << " m"); |
|
INFO("Distance drift: " << distance_drift_pct << "%"); |
|
INFO("Initial a: " << initial_a << " m"); |
|
INFO("Final a: " << final_a << " m"); |
|
INFO("a drift: " << a_drift_pct << "%"); |
|
INFO("Initial e: " << initial_e); |
|
INFO("Final e: " << final_e); |
|
|
|
REQUIRE_THAT(distance_drift_pct, WithinAbs(0.0, 0.01)); |
|
REQUIRE_THAT(a_drift_pct, WithinAbs(0.0, 0.01)); |
|
REQUIRE_THAT(final_e, WithinAbs(initial_e, 1e-6)); |
|
} |
|
|
|
SECTION("state vectors at orbital quarters") { |
|
const double orbit_radius = vec3_distance(craft->local_position, (Vec3){0,0,0}); |
|
const double earth_mass = earth->mass; |
|
const double orbital_period = 2.0 * M_PI * sqrt(pow(orbit_radius, 3.0) / (G * earth_mass)); |
|
const double quarter_orbit_time = orbital_period / 4.0; |
|
const int steps_per_quarter = (int)(quarter_orbit_time / TIME_STEP); |
|
|
|
INFO("Orbital radius: " << orbit_radius << " m"); |
|
INFO("Expected orbital period: " << orbital_period << " s (" << orbital_period / 3600.0 << " hours)"); |
|
INFO("Steps per quarter: " << steps_per_quarter); |
|
|
|
double previous_angle = atan2(craft->local_position.y, craft->local_position.x); |
|
|
|
for (int quarter = 0; quarter <= 4; quarter++) { |
|
INFO(""); |
|
INFO("=== Point " << quarter << "/4 (" << (quarter * 90) << " deg) ==="); |
|
INFO("Local position: (" << craft->local_position.x << ", " << craft->local_position.y << ", " << craft->local_position.z << ") m"); |
|
INFO("Local velocity: (" << craft->local_velocity.x << ", " << craft->local_velocity.y << ", " << craft->local_velocity.z << ") m/s"); |
|
|
|
double current_radius = vec3_distance(craft->local_position, (Vec3){0,0,0}); |
|
double current_velocity = vec3_magnitude(craft->local_velocity); |
|
double current_angle = atan2(craft->local_position.y, craft->local_position.x); |
|
|
|
INFO("Radius: " << current_radius << " m"); |
|
INFO("Velocity magnitude: " << current_velocity << " m/s"); |
|
INFO("Angular position: " << current_angle << " rad (" << (current_angle * 180.0 / M_PI) << " deg)"); |
|
|
|
if (quarter > 0) { |
|
double angle_change = current_angle - previous_angle; |
|
if (angle_change < 0) angle_change += 2.0 * M_PI; |
|
INFO("Angle change from previous: " << angle_change << " rad (" << (angle_change * 180.0 / M_PI) << " deg)"); |
|
// Each quarter advances ~89.5953° (23 steps of 60s on 5544.93s orbit) |
|
REQUIRE_THAT(angle_change, WithinAbs(1.56373, 0.001)); |
|
} |
|
|
|
if (quarter < 4) { |
|
for (int step = 0; step < steps_per_quarter; step++) { |
|
update_simulation(sim); |
|
} |
|
} |
|
|
|
previous_angle = current_angle; |
|
} |
|
|
|
INFO(""); |
|
INFO("=== Final Summary ==="); |
|
double final_radius = vec3_distance(craft->local_position, (Vec3){0,0,0}); |
|
double final_velocity = vec3_magnitude(craft->local_velocity); |
|
double final_angle = atan2(craft->local_position.y, craft->local_position.x); |
|
double total_rotation = final_angle; |
|
|
|
if (total_rotation < 0) total_rotation += 2.0 * M_PI; |
|
|
|
INFO("Total rotation: " << total_rotation << " rad (" << (total_rotation * 180.0 / M_PI) << " deg)"); |
|
INFO("Radius change: " << ((final_radius - orbit_radius) / orbit_radius * 100.0) << "%"); |
|
INFO("Velocity change: " << ((final_velocity - vec3_magnitude(craft->local_velocity)) / vec3_magnitude(craft->local_velocity) * 100.0) << "%"); |
|
|
|
// 92 steps of 60s = 5520s on 5544.93s orbit; ~358.38° total |
|
REQUIRE_THAT(total_rotation, WithinAbs(6.25493, 0.001)); |
|
} |
|
|
|
destroy_simulation(sim); |
|
}
|
|
|