@ -30,57 +30,57 @@ double calculate_orbital_period(CelestialBody* body, CelestialBody* parent) {
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 ;
}
@ -88,13 +88,13 @@ 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 ;
@ -103,10 +103,10 @@ double find_max_stable_dt(SimulationState* sim, const TestBody& test_body) {
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 ;
}
@ -116,33 +116,33 @@ void print_summary(const StabilityResult* results, int num_results, double min_s
printf ( " =============================================================================== \n " ) ;
printf ( " TIME STEP STABILITY TEST RESULTS \n " ) ;
printf ( " =============================================================================== \n \n " ) ;
printf ( " STABILITY CRITERIA: \n " ) ;
printf ( " - Energy drift < %.1f%% over %d orbits \n " , ENERGY_TOLERANCE , NUM_ORBITS ) ;
printf ( " - Distance drift < 5.0%% \n " ) ;
printf ( " - No SOI transitions (parent changes) \n \n " ) ;
printf ( " PER-BODY RESULTS: \n " ) ;
printf ( " +----------------------+----------------+------------------+----------------+ \n " ) ;
printf ( " | Body | Period (days) | Max Stable dt (s) | Stability Status | \n " ) ;
printf ( " +----------------------+----------------+------------------+----------------+ \n " ) ;
for ( int i = 0 ; i < num_results ; i + + ) {
const StabilityResult & r = results [ i ] ;
double ratio = default_dt / r . max_stable_dt ;
const char * status = ratio < 0.5 ? " Very Stable " : ratio < 0.8 ? " Stable " : " Limited Margin " ;
printf ( " | %-20s | %14.2f | %16.0f | %-14s | \n " , r . name , r . period_days , r . max_stable_dt , status ) ;
}
printf ( " +----------------------+----------------+------------------+----------------+ \n \n " ) ;
printf ( " OVERALL ANALYSIS: \n " ) ;
printf ( " Minimum stable time step: %.0f seconds \n " , min_stable_dt ) ;
printf ( " Recommended safe dt: %.0f seconds (0.7x safety margin) \n " , min_stable_dt * 0.7 ) ;
printf ( " Current default dt: %.0f seconds \n " , default_dt ) ;
printf ( " Current dt stability: %.0fx \n " , default_dt / min_stable_dt ) ;
printf ( " \n " ) ;
if ( default_dt < min_stable_dt * 0.7 ) {
printf ( " STATUS: Current time step (60s) is VERY STABLE with good margin. \n " ) ;
printf ( " Can be increased significantly if needed. \n \n " ) ;
@ -153,7 +153,7 @@ void print_summary(const StabilityResult* results, int num_results, double min_s
printf ( " STATUS: Current time step (60s) is near stability limit. \n " ) ;
printf ( " Consider reducing for safety. \n \n " ) ;
}
printf ( " RECOMMENDATIONS: \n " ) ;
printf ( " - For MESSENGER-like close orbits: Keep dt <= %.0f seconds \n " , min_stable_dt ) ;
printf ( " - For planetary missions: Current dt=60s is excellent \n " ) ;
@ -165,60 +165,60 @@ void print_summary(const StabilityResult* results, int num_results, double min_s
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 }
} ;
StabilityResult results [ 3 ] ;
double max_dt = MAX_DT ;
printf ( " \n === Finding minimum stable dt across all bodies === \n \n " ) ;
for ( int i = 0 ; i < 3 ; i + + ) {
results [ i ] . name = bodies [ i ] . name ;
results [ i ] . period_days = bodies [ i ] . expected_period_days ;
@ -227,70 +227,70 @@ TEST_CASE("Find minimum stable time step across all bodies", "[timestep][stabili
max_dt = results [ i ] . max_stable_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 " ) ;
print_summary ( results , 3 , max_dt , BASE_DT ) ;
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 ) ;
}