diff --git a/src/bodies.cpp b/src/bodies.cpp
index f594ce0..3161066 100644
--- a/src/bodies.cpp
+++ b/src/bodies.cpp
@@ -98,11 +98,36 @@ void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_ax
 
 // Update the entire simulation by one time step
 void update_simulation(SimulationState* sim) {
-    // Update each body's physics (except the root body which is stationary)
+    // First, update root bodies (they interact with each other)
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];
 
-        // Skip if this is a root body (parent_index == -1)
+        if (body->parent_index == -1) {
+            // This is a root body - calculate forces from OTHER root bodies
+            Vec3 total_force = {0.0, 0.0, 0.0};
+
+            for (int j = 0; j < sim->body_count; j++) {
+                if (i == j) continue;  // Don't apply force to itself
+
+                CelestialBody* other = &sim->bodies[j];
+                if (other->parent_index == -1) {
+                    // Other is also a root body - apply gravitational force
+                    Vec3 force = calculate_gravity_force(body, other);
+                    total_force = vec3_add(total_force, force);
+                }
+            }
+
+            // Apply total force from all other root bodies
+            Vec3 acceleration = calculate_acceleration(total_force, body->mass);
+            euler_step(body, acceleration, sim->dt);
+        }
+    }
+
+    // Now update non-root bodies (planets, moons, etc.)
+    for (int i = 0; i < sim->body_count; i++) {
+        CelestialBody* body = &sim->bodies[i];
+
+        // Skip root bodies (already updated above)
         if (body->parent_index == -1) {
             continue;
         }
diff --git a/src/config_loader.cpp b/src/config_loader.cpp
index 2db092f..cb0828e 100644
--- a/src/config_loader.cpp
+++ b/src/config_loader.cpp
@@ -61,12 +61,97 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
 // Calculate circular orbit velocity: v = sqrt(G * M / r)
 // Velocity is perpendicular to position vector
 void calculate_initial_velocities(SimulationState* sim) {
+    // First, handle multiple root bodies (binary stars, etc.)
+    // Find all root bodies and calculate barycentric orbits
+    int root_count = 0;
+    int root_indices[32];  // Max 32 root bodies
+    Vec3 barycenter = {0.0, 0.0, 0.0};
+    double total_mass = 0.0;
+
+    // Find all root bodies and calculate barycenter
+    for (int i = 0; i < sim->body_count; i++) {
+        if (sim->bodies[i].parent_index == -1) {
+            if (root_count < 32) {
+                root_indices[root_count++] = i;
+                // Weighted sum for barycenter
+                Vec3 weighted_pos = vec3_scale(sim->bodies[i].position, sim->bodies[i].mass);
+                barycenter = vec3_add(barycenter, weighted_pos);
+                total_mass += sim->bodies[i].mass;
+            }
+        }
+    }
+
+    // Calculate barycenter position
+    if (total_mass > 0.0) {
+        barycenter = vec3_scale(barycenter, 1.0 / total_mass);
+    }
+
+    // Debug output for multiple root bodies
+    if (root_count > 1) {
+        printf("\nBinary/Multiple star system detected:\n");
+        printf("  Number of root bodies: %d\n", root_count);
+        printf("  Barycenter position: (%.3e, %.3e, %.3e) m\n",
+               barycenter.x, barycenter.y, barycenter.z);
+        printf("  Total system mass: %.3e kg\n", total_mass);
+    }
+
+    // Set velocities for root bodies to orbit barycenter
+    if (root_count > 1) {
+        for (int i = 0; i < root_count; i++) {
+            CelestialBody* body = &sim->bodies[root_indices[i]];
+
+            // Calculate position relative to barycenter
+            Vec3 r = vec3_sub(body->position, barycenter);
+            double distance = vec3_magnitude(r);
+
+            if (distance < 1.0) {
+                body->velocity = {0.0, 0.0, 0.0};
+                continue;
+            }
+
+            // Calculate total mass of OTHER root bodies
+            double other_mass = total_mass - body->mass;
+
+            // Calculate circular orbit speed around barycenter
+            // v = sqrt(G * M_other / r)
+            double speed = sqrt(G * other_mass / distance);
+
+            // Create velocity perpendicular to position vector (same logic as below)
+            Vec3 z_axis = {0.0, 0.0, 1.0};
+            Vec3 vel_dir = {
+                r.y * z_axis.z - r.z * z_axis.y,
+                r.z * z_axis.x - r.x * z_axis.z,
+                r.x * z_axis.y - r.y * z_axis.x
+            };
+
+            // If r is parallel to z-axis, use x-axis instead
+            double cross_mag = vec3_magnitude(vel_dir);
+            if (cross_mag < 0.01) {
+                Vec3 x_axis = {1.0, 0.0, 0.0};
+                vel_dir.x = r.y * x_axis.z - r.z * x_axis.y;
+                vel_dir.y = r.z * x_axis.x - r.x * x_axis.z;
+                vel_dir.z = r.x * x_axis.y - r.y * x_axis.x;
+            }
+
+            // Normalize and scale by orbital speed
+            vel_dir = vec3_normalize(vel_dir);
+            body->velocity = vec3_scale(vel_dir, speed);
+
+            // Debug output
+            printf("  %s: distance from barycenter = %.3e m, orbital speed = %.3e m/s\n",
+                   body->name, distance, speed);
+        }
+    } else if (root_count == 1) {
+        // Single root body stays stationary
+        sim->bodies[root_indices[0]].velocity = {0.0, 0.0, 0.0};
+    }
+
+    // Now handle child bodies (planets, moons, etc.)
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];
 
-        // Skip root body (no parent)
+        // Skip root bodies (already handled)
         if (body->parent_index == -1) {
-            body->velocity = {0.0, 0.0, 0.0};
             continue;
         }