@ -13,6 +13,12 @@ struct TestBody {
double expected_period_days ;
} ;
struct StabilityResult {
const char * name ;
double max_stable_dt ;
double period_days ;
} ;
const double SECONDS_PER_DAY = 86400.0 ;
const double MIN_DT = 30.0 ;
const double MAX_DT = 600.0 ;
@ -105,6 +111,57 @@ double find_max_stable_dt(SimulationState* sim, const TestBody& test_body) {
return max_stable ;
}
void print_summary ( const StabilityResult * results , int num_results , double min_stable_dt , double default_dt ) {
printf ( " \n " ) ;
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 " ) ;
} else if ( default_dt < min_stable_dt ) {
printf ( " STATUS: Current time step (60s) is STABLE with adequate margin. \n " ) ;
printf ( " Moderate increases possible. \n \n " ) ;
} else {
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 " ) ;
printf ( " - For Moon-scale orbits: Could use dt=120s+ safely \n " ) ;
printf ( " \n " ) ;
printf ( " =============================================================================== \n \n " ) ;
}
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 ) ;
@ -157,14 +214,17 @@ TEST_CASE("Find minimum stable time step across all bodies", "[timestep][stabili
{ " 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 + + ) {
double body_max_dt = find_max_stable_dt ( sim , bodies [ i ] ) ;
if ( body_max_dt < max_dt ) {
max_dt = body_max_dt ;
results [ i ] . name = bodies [ i ] . name ;
results [ i ] . period_days = bodies [ i ] . expected_period_days ;
results [ i ] . max_stable_dt = find_max_stable_dt ( sim , bodies [ i ] ) ;
if ( results [ i ] . max_stable_dt < max_dt ) {
max_dt = results [ i ] . max_stable_dt ;
}
}
@ -174,6 +234,8 @@ TEST_CASE("Find minimum stable time step across all bodies", "[timestep][stabili
INFO ( " Minimum stable dt: " < < max_dt < < " seconds " ) ;
print_summary ( results , 3 , max_dt , BASE_DT ) ;
REQUIRE ( max_dt > = MIN_DT ) ;
destroy_simulation ( sim ) ;