Fix test_analytical_propagation_timesteps: tolerance, division by zero, and 2π wrapping

tests/test_analytical_propagation_timesteps.cpp

@@ -94,13 +94,15 @@ TEST_CASE("Very small timestep - dt less than 1 second", "[analytical][timestep]
     INFO("Position change: " << pos_change << " m");
     INFO("Velocity change: " << vel_change << " m/s");
 
-    double expected_pos_change = vel_change * small_dt;
+    double v_before_mag = vec3_magnitude(vel_before);
+    double expected_pos_change = v_before_mag * small_dt;
     double pos_error = fabs(pos_change - expected_pos_change);
 
-    INFO("Expected position change: " << expected_pos_change << " m");
+    INFO("Expected position change (v·dt): " << expected_pos_change << " m");
     INFO("Position error: " << pos_error << " m");
+    INFO("Relative position error: " << (pos_error / expected_pos_change * 100.0) << "%");
 
-    REQUIRE(pos_change < VELOCITY_TOLERANCE * small_dt * 10.0);
+    REQUIRE(pos_error < VELOCITY_TOLERANCE * small_dt * 10.0);
     REQUIRE(vel_change < VELOCITY_TOLERANCE);
 
     destroy_simulation(sim);
@@ -149,12 +151,15 @@ TEST_CASE("Accuracy vs timestep size relationship", "[analytical][timestep][accu
         INFO("  Expected error (phase): " << expected_pos_error << " m");
         INFO("  Number of periods: " << num_periods);
 
-        if (expected_num_orbits > 0) {
+        if (expected_num_orbits > 0 && expected_pos_error > 1.0e-6) {
             double relative_error = pos_error / expected_pos_error;
 
             INFO("  Relative error: " << relative_error);
 
             REQUIRE(relative_error < 0.5);
+        } else if (expected_num_orbits > 0) {
+            INFO("  Expected error is zero, skipping relative error check");
+            REQUIRE(pos_error < 1.0e-3);
         }
     }
 
@@ -192,7 +197,8 @@ TEST_CASE("Mean anomaly accumulation over long propagation", "[analytical][times
 
     INFO("Expected final anomaly: " << expected_final_anomaly << " rad");
 
-    double anomaly_error = fabs(final_true_anomaly - expected_final_anomaly);
+    double raw_error = fabs(final_true_anomaly - expected_final_anomaly);
+    double anomaly_error = fmin(raw_error, 2.0 * M_PI - raw_error);
 
     INFO("True anomaly error: " << anomaly_error << " rad (" << anomaly_error * 180.0 / M_PI << "°)");