From a64e1ab17c10d4da77c933d511508eb013d475ed Mon Sep 17 00:00:00 2001
From: cinnaboot <cinnaboot@protonmail.com>
Date: Fri, 9 Jan 2026 09:17:05 -0500
Subject: [PATCH] Remove multiple root bodies support from simulation.cpp

Removed barycentric orbit calculations for binary/multiple star systems.
Simplified calculate_initial_velocities() to only support single root body.

- Removed: compute_system_barycenter()
- Removed: compute_total_root_mass()
- Removed: set_root_bodies_velocity()
- Removed: print_system_info_if_multiple_roots()
- Simplified: calculate_initial_velocities() now sets root body velocity to zero

Remove unused multiple root body support
---
 src/simulation.cpp | 97 ++--------------------------------------------
 1 file changed, 4 insertions(+), 93 deletions(-)

diff --git a/src/simulation.cpp b/src/simulation.cpp
index 8e2575e..5c3aa7c 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -113,29 +113,6 @@ void update_simulation(SimulationState* sim) {
     sim->time += sim->dt;
 }
 
-static void compute_perpendicular_orbital_velocity(CelestialBody* body, Vec3 center,
-                                                   double orbiting_mass) {
-    Vec3 r = vec3_sub(body->position, center);
-    double distance = vec3_magnitude(r);
-
-    if (distance < 1.0) {
-        body->velocity = {0.0, 0.0, 0.0};
-        return;
-    }
-
-    double speed = sqrt(G * orbiting_mass / distance);
-    Vec3 z_axis = {0.0, 0.0, 1.0};
-    Vec3 vel_dir = vec3_cross(r, z_axis);
-
-    if (vec3_magnitude(vel_dir) < 0.01) {
-        Vec3 x_axis = {1.0, 0.0, 0.0};
-        vel_dir = vec3_cross(r, x_axis);
-    }
-
-    vel_dir = vec3_normalize(vel_dir);
-    body->velocity = vec3_scale(vel_dir, speed);
-}
-
 static void compute_orbital_velocity_from_vis_viva(CelestialBody* body,
                                                      CelestialBody* parent) {
     Vec3 r = vec3_sub(body->position, parent->position);
@@ -168,56 +145,6 @@ static void compute_orbital_velocity_from_vis_viva(CelestialBody* body,
     body->velocity = vec3_add(body->velocity, parent->velocity);
 }
 
-static Vec3 compute_system_barycenter(SimulationState* sim, int* root_indices,
-                                      int* root_count) {
-    Vec3 barycenter = {0.0, 0.0, 0.0};
-    *root_count = 0;
-
-    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;
-                Vec3 weighted_pos = vec3_scale(sim->bodies[i].position, sim->bodies[i].mass);
-                barycenter = vec3_add(barycenter, weighted_pos);
-            }
-        }
-    }
-
-    double total_mass = 0.0;
-    for (int i = 0; i < *root_count; i++) {
-        total_mass += sim->bodies[root_indices[i]].mass;
-    }
-
-    if (total_mass > 0.0) {
-        barycenter = vec3_scale(barycenter, 1.0 / total_mass);
-    }
-
-    return barycenter;
-}
-
-static double compute_total_root_mass(SimulationState* sim, int* root_indices,
-                                      int root_count) {
-    double total_mass = 0.0;
-    for (int i = 0; i < root_count; i++) {
-        total_mass += sim->bodies[root_indices[i]].mass;
-    }
-    return total_mass;
-}
-
-static void set_root_bodies_velocity(SimulationState* sim, int* root_indices,
-                                     int root_count, Vec3 barycenter, double total_mass) {
-    for (int i = 0; i < root_count; i++) {
-        CelestialBody* body = &sim->bodies[root_indices[i]];
-        double other_mass = total_mass - body->mass;
-        compute_perpendicular_orbital_velocity(body, barycenter, other_mass);
-
-        double distance = vec3_magnitude(vec3_sub(body->position, barycenter));
-        double speed = vec3_magnitude(body->velocity);
-        printf("  %s: distance from barycenter = %.3e m, orbital speed = %.3e m/s\n",
-               body->name, distance, speed);
-    }
-}
-
 static void set_child_bodies_velocity(SimulationState* sim) {
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];
@@ -233,27 +160,11 @@ static void set_child_bodies_velocity(SimulationState* sim) {
     }
 }
 
-static void print_system_info_if_multiple_roots(int root_count, Vec3 barycenter) {
-    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);
-    }
-}
-
 void calculate_initial_velocities(SimulationState* sim) {
-    int root_indices[32];
-    int root_count;
-    Vec3 barycenter = compute_system_barycenter(sim, root_indices, &root_count);
-
-    print_system_info_if_multiple_roots(root_count, barycenter);
-
-    if (root_count > 1) {
-        double total_mass = compute_total_root_mass(sim, root_indices, root_count);
-        set_root_bodies_velocity(sim, root_indices, root_count, barycenter, total_mass);
-    } else if (root_count == 1) {
-        sim->bodies[root_indices[0]].velocity = {0.0, 0.0, 0.0};
+    for (int i = 0; i < sim->body_count; i++) {
+        if (sim->bodies[i].parent_index == -1) {
+            sim->bodies[i].velocity = {0.0, 0.0, 0.0};
+        }
     }
 
     set_child_bodies_velocity(sim);