|
|
|
@ -6,6 +6,8 @@ |
|
|
|
#include <cmath> |
|
|
|
#include <cmath> |
|
|
|
#include <vector> |
|
|
|
#include <vector> |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#define COMET_TEST_VERBOSE_OUTPUT 0 |
|
|
|
|
|
|
|
|
|
|
|
struct ParentChange { |
|
|
|
struct ParentChange { |
|
|
|
double time_seconds; |
|
|
|
double time_seconds; |
|
|
|
double time_days; |
|
|
|
double time_days; |
|
|
|
@ -15,14 +17,65 @@ struct ParentChange { |
|
|
|
double distance_to_sun_au; |
|
|
|
double distance_to_sun_au; |
|
|
|
}; |
|
|
|
}; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#if COMET_TEST_VERBOSE_OUTPUT |
|
|
|
|
|
|
|
static void print_orbital_snapshots(const std::vector<OrbitalElements>& snapshots, |
|
|
|
|
|
|
|
double expected_sma, double expected_ecc) { |
|
|
|
|
|
|
|
printf("\n=== Comet Orbital Elements Over Time ===\n"); |
|
|
|
|
|
|
|
printf("Expected: a = %.1f AU, e = %.1f\n\n", expected_sma, expected_ecc); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
for (const auto& elem : snapshots) { |
|
|
|
|
|
|
|
printf("Day %.0f:\n", elem.time_days); |
|
|
|
|
|
|
|
printf(" Semi-major axis: %.6f AU (expected: %.1f AU)\n", elem.semi_major_axis_au, expected_sma); |
|
|
|
|
|
|
|
printf(" Eccentricity: %.6f (expected: %.1f)\n", elem.eccentricity, expected_ecc); |
|
|
|
|
|
|
|
printf(" Distance to Sun: %.4f AU\n", elem.distance_to_sun_au); |
|
|
|
|
|
|
|
printf(" Distance to Mars: %.6f AU\n", elem.distance_to_ref_body_au); |
|
|
|
|
|
|
|
printf(" Velocity: %.2f km/s\n", elem.velocity_magnitude / 1000.0); |
|
|
|
|
|
|
|
printf(" Specific energy: %.3e J/kg\n", elem.specific_energy); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
double sma_error = fabs(elem.semi_major_axis_au - expected_sma); |
|
|
|
|
|
|
|
double ecc_error = fabs(elem.eccentricity - expected_ecc); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if (elem.time_days >= 1530 && elem.time_days <= 1540) { |
|
|
|
|
|
|
|
printf(" *** CLOSEST APPROACH TO MARS ***\n"); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
printf(" Error: Δa = %.6f AU, Δe = %.6f\n\n", sma_error, ecc_error); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
static void print_parent_changes(const std::vector<ParentChange>& parent_changes, |
|
|
|
|
|
|
|
int mars_index) { |
|
|
|
|
|
|
|
printf("\n=== Parent Index Changes ===\n"); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if (parent_changes.empty()) { |
|
|
|
|
|
|
|
printf("No parent changes detected\n\n"); |
|
|
|
|
|
|
|
} else { |
|
|
|
|
|
|
|
for (const auto& change : parent_changes) { |
|
|
|
|
|
|
|
printf("Time: %.0f days (%.0f seconds)\n", change.time_days, change.time_seconds); |
|
|
|
|
|
|
|
printf(" Parent: %d -> %d", change.old_parent, change.new_parent); |
|
|
|
|
|
|
|
if (change.new_parent == mars_index) { |
|
|
|
|
|
|
|
printf(" (Sun -> Mars)"); |
|
|
|
|
|
|
|
} else if (change.old_parent == mars_index) { |
|
|
|
|
|
|
|
printf(" (Mars -> Sun)"); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
printf("\n"); |
|
|
|
|
|
|
|
printf(" Distance to Mars: %.6f AU\n", change.distance_to_mars_au); |
|
|
|
|
|
|
|
printf(" Distance to Sun: %.6f AU\n\n", change.distance_to_sun_au); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
#endif |
|
|
|
|
|
|
|
|
|
|
|
TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[comet][orbital-elements][soi]") { |
|
|
|
TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[comet][orbital-elements][soi]") { |
|
|
|
const double TIME_STEP = 60.0; |
|
|
|
const double TIME_STEP = 60.0; |
|
|
|
const double SECONDS_PER_DAY = 86400.0; |
|
|
|
const double SECONDS_PER_DAY = 86400.0; |
|
|
|
const double DAYS_TO_SIMULATE = 5000.0; |
|
|
|
const double DAYS_TO_SIMULATE = 5000.0; |
|
|
|
const double AU = 1.496e11; |
|
|
|
const double AU = 1.496e11; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#if COMET_TEST_VERBOSE_OUTPUT |
|
|
|
const double EXPECTED_SMA = 2.5; |
|
|
|
const double EXPECTED_SMA = 2.5; |
|
|
|
const double EXPECTED_ECC = 0.7; |
|
|
|
const double EXPECTED_ECC = 0.7; |
|
|
|
|
|
|
|
#endif |
|
|
|
|
|
|
|
|
|
|
|
SimulationState* sim = create_simulation(10, TIME_STEP); |
|
|
|
SimulationState* sim = create_simulation(10, TIME_STEP); |
|
|
|
|
|
|
|
|
|
|
|
@ -79,48 +132,11 @@ TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[ |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
#if 1 |
|
|
|
#if COMET_TEST_VERBOSE_OUTPUT |
|
|
|
printf("\n=== Comet Orbital Elements Over Time ===\n"); |
|
|
|
print_orbital_snapshots(snapshots, EXPECTED_SMA, EXPECTED_ECC); |
|
|
|
printf("Expected: a = %.1f AU, e = %.1f\n\n", EXPECTED_SMA, EXPECTED_ECC); |
|
|
|
print_parent_changes(parent_changes, MARS_INDEX); |
|
|
|
|
|
|
|
|
|
|
|
for (const auto& elem : snapshots) { |
|
|
|
|
|
|
|
printf("Day %.0f:\n", elem.time_days); |
|
|
|
|
|
|
|
printf(" Semi-major axis: %.6f AU (expected: %.1f AU)\n", elem.semi_major_axis_au, EXPECTED_SMA); |
|
|
|
|
|
|
|
printf(" Eccentricity: %.6f (expected: %.1f)\n", elem.eccentricity, EXPECTED_ECC); |
|
|
|
|
|
|
|
printf(" Distance to Sun: %.4f AU\n", elem.distance_to_sun_au); |
|
|
|
|
|
|
|
printf(" Distance to Mars: %.6f AU\n", elem.distance_to_ref_body_au); |
|
|
|
|
|
|
|
printf(" Velocity: %.2f km/s\n", elem.velocity_magnitude / 1000.0); |
|
|
|
|
|
|
|
printf(" Specific energy: %.3e J/kg\n", elem.specific_energy); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
double sma_error = fabs(elem.semi_major_axis_au - EXPECTED_SMA); |
|
|
|
|
|
|
|
double ecc_error = fabs(elem.eccentricity - EXPECTED_ECC); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if (elem.time_days >= 1530 && elem.time_days <= 1540) { |
|
|
|
|
|
|
|
printf(" *** CLOSEST APPROACH TO MARS ***\n"); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
printf(" Error: Δa = %.6f AU, Δe = %.6f\n\n", sma_error, ecc_error); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
printf("\n=== Parent Index Changes ===\n"); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if (parent_changes.empty()) { |
|
|
|
|
|
|
|
printf("No parent changes detected\n\n"); |
|
|
|
|
|
|
|
} else { |
|
|
|
|
|
|
|
for (const auto& change : parent_changes) { |
|
|
|
|
|
|
|
printf("Time: %.0f days (%.0f seconds)\n", change.time_days, change.time_seconds); |
|
|
|
|
|
|
|
printf(" Parent: %d -> %d", change.old_parent, change.new_parent); |
|
|
|
|
|
|
|
if (change.new_parent == MARS_INDEX) { |
|
|
|
|
|
|
|
printf(" (Sun -> Mars)"); |
|
|
|
|
|
|
|
} else if (change.old_parent == MARS_INDEX) { |
|
|
|
|
|
|
|
printf(" (Mars -> Sun)"); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
printf("\n"); |
|
|
|
|
|
|
|
printf(" Distance to Mars: %.6f AU\n", change.distance_to_mars_au); |
|
|
|
|
|
|
|
printf(" Distance to Sun: %.6f AU\n\n", change.distance_to_sun_au); |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
} |
|
|
|
|
|
|
|
#endif |
|
|
|
#endif |
|
|
|
|
|
|
|
|
|
|
|
REQUIRE(parent_changes.size() > 0); |
|
|
|
REQUIRE(parent_changes.size() > 0); |
|
|
|
|
|
|
|
|
|
|
|
bool found_mars_transition = false; |
|
|
|
bool found_mars_transition = false; |
|
|
|
|