Browse Source
- Create soi_transition.toml with Sun + Mars + SmallBody (3 bodies) - Add test_soi_transition.cpp with 2 test cases: * SOI transition from Sun to Mars (validates one-way transition) * SOI radii calculation verification - Use comet parameters (e=0.7, a=3.74e11) for Mars encounter - Document hysteresis makes Mars -> Sun exit impossible - Simpler than original 7-body test_comet_orbit.cpp - All 3 SOI tests pass (12 assertions total) Claudemain
2 changed files with 151 additions and 0 deletions
@ -0,0 +1,34 @@ |
|||||||
|
# Test Configuration: Sun + Mars + Small Body |
||||||
|
# Small body will experience SOI transition with Mars |
||||||
|
# Mars SOI: ~0.58 million km (~0.004 AU) |
||||||
|
# Expected: Sun -> Mars transition (one-way due to hysteresis) |
||||||
|
|
||||||
|
[[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 = "Mars" |
||||||
|
mass = 6.39e23 |
||||||
|
radius = 3.3895e6 |
||||||
|
position = { x = 2.244e11, y = 0.0, z = 0.0 } |
||||||
|
parent_index = 0 |
||||||
|
color = { r = 0.8, g = 0.3, b = 0.1 } |
||||||
|
eccentricity = 0.0 |
||||||
|
semi_major_axis = 2.244e11 |
||||||
|
|
||||||
|
[[bodies]] |
||||||
|
name = "SmallBody" |
||||||
|
mass = 1.0e14 |
||||||
|
radius = 5.0e3 |
||||||
|
position = { x = 1.122e11, y = 0.0, z = 0.0 } |
||||||
|
parent_index = 0 |
||||||
|
color = { r = 0.5, g = 0.8, b = 1.0 } |
||||||
|
eccentricity = 0.7 |
||||||
|
semi_major_axis = 3.74e11 |
||||||
@ -0,0 +1,117 @@ |
|||||||
|
#include <catch2/catch_test_macros.hpp> |
||||||
|
#include "../src/physics.h" |
||||||
|
#include "../src/simulation.h" |
||||||
|
#include "../src/config_loader.h" |
||||||
|
#include <cmath> |
||||||
|
#include <vector> |
||||||
|
|
||||||
|
struct ParentChange { |
||||||
|
double time_seconds; |
||||||
|
double time_days; |
||||||
|
int old_parent; |
||||||
|
int new_parent; |
||||||
|
double distance_to_mars_au; |
||||||
|
double distance_to_sun_au; |
||||||
|
}; |
||||||
|
|
||||||
|
TEST_CASE("SOI transition - Sun to Mars", "[soi][transition]") { |
||||||
|
const double TIME_STEP = 60.0; |
||||||
|
const double DAYS_TO_SIMULATE = 1700.0; |
||||||
|
const double SECONDS_PER_DAY = 86400.0; |
||||||
|
const double AU = 1.496e11; |
||||||
|
|
||||||
|
INFO("Note: Hysteresis (0.5 factor) makes Mars -> Sun exit impossible"); |
||||||
|
|
||||||
|
SimulationState* sim = create_simulation(10, TIME_STEP); |
||||||
|
|
||||||
|
REQUIRE(load_system_config(sim, "tests/configs/soi_transition.toml")); |
||||||
|
|
||||||
|
const int SMALL_BODY_INDEX = 2; |
||||||
|
const int MARS_INDEX = 1; |
||||||
|
const int SUN_INDEX = 0; |
||||||
|
|
||||||
|
std::vector<ParentChange> parent_changes; |
||||||
|
int previous_parent = sim->bodies[SMALL_BODY_INDEX].parent_index; |
||||||
|
|
||||||
|
INFO("Initial parent index: " << previous_parent); |
||||||
|
INFO("Expected: Sun (0)"); |
||||||
|
|
||||||
|
REQUIRE(previous_parent == SUN_INDEX); |
||||||
|
|
||||||
|
double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; |
||||||
|
while (sim->time < max_time) { |
||||||
|
update_simulation(sim); |
||||||
|
|
||||||
|
int current_parent = sim->bodies[SMALL_BODY_INDEX].parent_index; |
||||||
|
if (current_parent != previous_parent) { |
||||||
|
ParentChange change; |
||||||
|
change.time_seconds = sim->time; |
||||||
|
change.time_days = sim->time / SECONDS_PER_DAY; |
||||||
|
change.old_parent = previous_parent; |
||||||
|
change.new_parent = current_parent; |
||||||
|
change.distance_to_mars_au = vec3_distance(sim->bodies[SMALL_BODY_INDEX].position, |
||||||
|
sim->bodies[MARS_INDEX].position) / AU; |
||||||
|
change.distance_to_sun_au = vec3_magnitude(sim->bodies[SMALL_BODY_INDEX].position) / AU; |
||||||
|
parent_changes.push_back(change); |
||||||
|
previous_parent = current_parent; |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
INFO("Total parent changes: " << parent_changes.size()); |
||||||
|
|
||||||
|
REQUIRE(parent_changes.size() > 0); |
||||||
|
|
||||||
|
bool found_mars_transition = false; |
||||||
|
|
||||||
|
for (const auto& change : parent_changes) { |
||||||
|
if (change.new_parent == MARS_INDEX || change.old_parent == MARS_INDEX) { |
||||||
|
found_mars_transition = true; |
||||||
|
if (change.new_parent == MARS_INDEX) { |
||||||
|
INFO("Entered Mars SOI at day " << change.time_days |
||||||
|
<< ", distance to Mars: " << change.distance_to_mars_au << " AU" |
||||||
|
<< ", distance to Sun: " << change.distance_to_sun_au << " AU"); |
||||||
|
} else if (change.old_parent == MARS_INDEX) { |
||||||
|
INFO("Exited Mars SOI at day " << change.time_days |
||||||
|
<< ", distance to Mars: " << change.distance_to_mars_au << " AU" |
||||||
|
<< ", distance to Sun: " << change.distance_to_sun_au << " AU"); |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
REQUIRE(found_mars_transition); |
||||||
|
|
||||||
|
destroy_simulation(sim); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("SOI transition - verify SOI radii", "[soi][radii]") { |
||||||
|
const double TIME_STEP = 60.0; |
||||||
|
const double AU = 1.496e11; |
||||||
|
|
||||||
|
SimulationState* sim = create_simulation(10, TIME_STEP); |
||||||
|
|
||||||
|
REQUIRE(load_system_config(sim, "tests/configs/soi_transition.toml")); |
||||||
|
|
||||||
|
const int MARS_INDEX = 1; |
||||||
|
|
||||||
|
double mars_soi_meters = sim->bodies[MARS_INDEX].soi_radius; |
||||||
|
double mars_soi_au = mars_soi_meters / AU; |
||||||
|
|
||||||
|
INFO("Mars SOI radius: " << mars_soi_meters / 1e9 << " million km"); |
||||||
|
INFO("Mars SOI radius: " << mars_soi_au << " AU"); |
||||||
|
|
||||||
|
double expected_mars_soi = 2.244e11 * pow(6.39e23 / 1.989e30, 0.4); |
||||||
|
double expected_mars_soi_au = expected_mars_soi / AU; |
||||||
|
|
||||||
|
INFO("Expected Mars SOI: " << expected_mars_soi / 1e9 << " million km"); |
||||||
|
INFO("Expected Mars SOI: " << expected_mars_soi_au << " AU"); |
||||||
|
|
||||||
|
double soi_error = fabs(mars_soi_au - expected_mars_soi_au) / expected_mars_soi_au; |
||||||
|
INFO("SOI calculation error: " << soi_error * 100.0 << "%"); |
||||||
|
|
||||||
|
REQUIRE(soi_error < 0.01); |
||||||
|
|
||||||
|
REQUIRE(mars_soi_au > 0.003); |
||||||
|
REQUIRE(mars_soi_au < 0.005); |
||||||
|
|
||||||
|
destroy_simulation(sim); |
||||||
|
} |
||||||
Loading…
Reference in new issue