Add eccentric orbit support with unified config format

Extend config format to support both circular and eccentric orbits using
eccentricity and semi-major axis parameters. All bodies now use the same
format: e=0 for circular orbits, e>0 for eccentric orbits.

Velocity calculation uses vis-viva equation for all cases. Add example
comet with highly eccentric orbit (e=0.7) to test configuration.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

configs/test_simple.txt

@@ -1,22 +1,28 @@
 # Simple Test Configuration
-# Two planets in circular orbits for easy position verification
+# Planets in circular and eccentric orbits
 #
 # Orbital periods (for verification):
 #   Earth: ~365 days
 #   Mars: ~687 days
+#   Comet: ~1444 days
 #
 # At t=0:
 #   Earth starts at (1.0 AU, 0, 0) = (1.496e11 m, 0, 0)
 #   Mars starts at (1.5 AU, 0, 0) = (2.244e11 m, 0, 0)
-#   Both orbit counter-clockwise when viewed from +Z
+#   Comet starts at perihelion (0.75 AU, 0, 0) = (1.122e11 m, 0, 0)
+#   All orbit counter-clockwise when viewed from +Z
 #
-# Format: name mass(kg) radius(m) x(m) y(m) z(m) parent_index r g b
+# Format: name mass(kg) radius(m) x(m) y(m) z(m) parent_index r g b eccentricity semi_major_axis(m)
 
 # Sun at origin (index 0)
-Sun 1.989e30 6.96e8 0 0 0 -1 1.0 1.0 0.0
+Sun 1.989e30 6.96e8 0 0 0 -1 1.0 1.0 0.0 0 0
 
-# Earth at 1 AU on positive X axis (index 1, parent Sun)
-Earth 5.972e24 6.371e6 1.496e11 0 0 0 0.0 0.5 1.0
+# Earth at 1 AU - circular orbit (index 1, parent Sun)
+Earth 5.972e24 6.371e6 1.496e11 0 0 0 0.0 0.5 1.0 0 1.496e11
 
-# Mars at 1.5 AU on positive X axis (index 2, parent Sun)
-Mars 6.39e23 3.3895e6 2.244e11 0 0 0 0.8 0.3 0.1
+# Mars at 1.5 AU - circular orbit (index 2, parent Sun)
+Mars 6.39e23 3.3895e6 2.244e11 0 0 0 0.8 0.3 0.1 0 2.244e11
+
+# Comet - eccentric orbit (e=0.7, a=2.5 AU)
+# At perihelion: 0.75 AU, aphelion: 4.25 AU
+Comet 1e14 5e3 1.122e11 0 0 0 0.5 0.8 1.0 0.7 3.74e11

src/config_loader.cpp

@@ -5,20 +5,29 @@
 
 // Parse a single body definition line
 bool parse_body_line(const char* line, char* name, double* mass, double* radius,
-                     Vec3* pos, int* parent_index, float* r, float* g, float* b) {
+                     Vec3* pos, int* parent_index, float* r, float* g, float* b,
+                     double* eccentricity, double* semi_major_axis) {
     // Skip empty lines and comments
     if (line[0] == '\0' || line[0] == '#' || line[0] == '\n') {
         return false;
     }
 
-    // Parse: name mass(kg) radius(m) x(m) y(m) z(m) parent_index r g b
-    int result = sscanf(line, "%63s %lf %lf %lf %lf %lf %d %f %f %f",
+    // Format: name mass radius x y z parent_index r g b eccentricity semi_major_axis
+    int result = sscanf(line, "%63s %lf %lf %lf %lf %lf %d %f %f %f %lf %lf",
                        name, mass, radius, &pos->x, &pos->y, &pos->z,
-                       parent_index, r, g, b);
+                       parent_index, r, g, b, eccentricity, semi_major_axis);
 
-    return result == 10;  // All fields must be present
+    return result == 12;
 }
 
+struct OrbitParams {
+    double eccentricity;
+    double semi_major_axis;
+};
+
+// Forward declaration
+void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params);
+
 // Load system configuration from file
 bool load_system_config(SimulationState* sim, const char* filepath) {
     FILE* file = fopen(filepath, "r");
@@ -33,11 +42,18 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
     Vec3 pos;
     int parent_index;
     float r, g, b;
+    double eccentricity, semi_major_axis;
+
+    OrbitParams orbit_params[100];
 
     while (fgets(line, sizeof(line), file)) {
-        if (parse_body_line(line, name, &mass, &radius, &pos, &parent_index, &r, &g, &b)) {
-            Vec3 vel = {0.0, 0.0, 0.0};  // Velocity will be calculated later
+        if (parse_body_line(line, name, &mass, &radius, &pos, &parent_index, &r, &g, &b,
+                           &eccentricity, &semi_major_axis)) {
+            Vec3 vel = {0.0, 0.0, 0.0};
             add_body(sim, name, mass, radius, pos, vel, parent_index, r, g, b);
+
+            orbit_params[sim->body_count - 1].eccentricity = eccentricity;
+            orbit_params[sim->body_count - 1].semi_major_axis = semi_major_axis;
         }
     }
 
@@ -48,8 +64,8 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
         return false;
     }
 
-    // Calculate initial velocities for circular orbits
-    calculate_initial_velocities(sim);
+    // Calculate initial velocities
+    calculate_velocities(sim, orbit_params);
 
     // Calculate SOI radii
     calculate_soi_radii(sim);
@@ -58,9 +74,8 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
     return true;
 }
 
-// Calculate circular orbit velocity: v = sqrt(G * M / r)
-// Velocity is perpendicular to position vector
-void calculate_initial_velocities(SimulationState* sim) {
+// Calculate initial velocities using vis-viva equation
+void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params) {
     // First, handle multiple root bodies (binary stars, etc.)
     // Find all root bodies and calculate barycentric orbits
     int root_count = 0;
@@ -168,8 +183,16 @@ void calculate_initial_velocities(SimulationState* sim) {
                 continue;
             }
 
-            // Calculate circular orbit speed
-            double speed = sqrt(G * parent->mass / distance);
+            double e = orbit_params[i].eccentricity;
+            double a = orbit_params[i].semi_major_axis;
+
+            // Use vis-viva equation: v = sqrt(G*M*(2/r - 1/a))
+            double speed = sqrt(G * parent->mass * (2.0 / distance - 1.0 / a));
+
+            if (e > 0.01) {
+                printf("  %s: eccentric orbit (e=%.2f, a=%.3e m), speed at r=%.3e m: %.3f km/s\n",
+                       body->name, e, a, distance, speed / 1000.0);
+            }
 
             // Create velocity perpendicular to position vector
             // If position is mostly in XY plane, make velocity in XY plane