diff --git a/CLAUDE.md b/CLAUDE.md
index 321ec53..b2e49ae 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -7,7 +7,6 @@
 ## Architecture
 - C-style C++ (structs/functions, NO classes/templates)
 - Raylib (git submodule) for 3D - chose over SFML (no 3D support)
-- See docs/verbose_project_overview.md for detailed technical architecture
 - See docs/implementation_plan.md for data structures reference
 
 ## Coding Rules
diff --git a/Makefile b/Makefile
index 27a7e09..9f10a77 100644
--- a/Makefile
+++ b/Makefile
@@ -14,7 +14,7 @@ TEST_TARGET = orbit_test
 # Source files
 CPP_SOURCES = $(SRC_DIR)/main.cpp \
                $(SRC_DIR)/physics.cpp \
-               $(SRC_DIR)/bodies.cpp \
+               $(SRC_DIR)/simulation.cpp \
                $(SRC_DIR)/config_loader.cpp \
                $(SRC_DIR)/renderer.cpp
 
@@ -86,12 +86,12 @@ test-build:
 	$(CXX) $(CXXFLAGS) \
 		-c $(SRC_DIR)/physics.cpp -o build/test_physics.o
 	$(CXX) $(CXXFLAGS) \
-		-c $(SRC_DIR)/bodies.cpp -o build/test_bodies.o
+		-c $(SRC_DIR)/simulation.cpp -o build/test_simulation.o
 	$(CXX) $(CXXFLAGS) \
 		-c $(SRC_DIR)/config_loader.cpp -o build/test_config_loader.o
 	$(CC) -Wall -Wextra -g -ggdb3 -std=c99 -I./src -I./ext/tomlc17/src \
 		-c $(SRC_DIR)/../ext/tomlc17/src/tomlc17.c -o build/test_tomlc17.o
-	$(CXX) build/test_main.o build/test_integration.o build/test_energy.o build/test_orbital_period.o build/test_comet_orbit.o build/test_utilities.o build/test_physics.o build/test_bodies.o build/test_config_loader.o build/test_tomlc17.o \
+	$(CXX) build/test_main.o build/test_integration.o build/test_energy.o build/test_orbital_period.o build/test_comet_orbit.o build/test_utilities.o build/test_physics.o build/test_simulation.o build/test_config_loader.o build/test_tomlc17.o \
 		-o $(TEST_TARGET) -lCatch2Main -lCatch2 -lm
 
 # Run automated test suite
diff --git a/docs/implementation_plan.md b/docs/implementation_plan.md
index 5f699b1..9ec771d 100644
--- a/docs/implementation_plan.md
+++ b/docs/implementation_plan.md
@@ -18,7 +18,7 @@ struct Vec3 {
 };
 ```
 
-### CelestialBody (bodies.h)
+### CelestialBody (simulation.h)
 ```cpp
 struct CelestialBody {
     char name[64];
@@ -34,7 +34,7 @@ struct CelestialBody {
 };
 ```
 
-### SimulationState (bodies.h)
+### SimulationState (simulation.h)
 ```cpp
 struct SimulationState {
     CelestialBody* bodies;
@@ -55,7 +55,7 @@ struct RenderState {
 };
 ```
 
-### OrbitalElements (bodies.h)
+### OrbitalElements (simulation.h)
 ```cpp
 struct OrbitalElements {
     double time_days;
@@ -82,7 +82,7 @@ struct AccelerationContext {
 ### Physics (physics.cpp/h)
 Vector math and gravity calculations. RK4 (Runge-Kutta 4th order) integration with `rk4_step()`.
 
-### Bodies (bodies.cpp/h)
+### Simulation (simulation.cpp/h)
 Simulation state management and updates. SOI detection using Hill sphere: `r_soi = a * (m/M)^(2/5)`.
 
 **Key functions:**
diff --git a/docs/session_summaries/2026-01-04-refactor-orbital-elements.md b/docs/session_summaries/2026-01-04-refactor-orbital-elements.md
index 592951c..799f582 100644
--- a/docs/session_summaries/2026-01-04-refactor-orbital-elements.md
+++ b/docs/session_summaries/2026-01-04-refactor-orbital-elements.md
@@ -3,8 +3,8 @@
 ## Changes Made
 
 Refactored orbital elements calculation from test-specific to reusable code:
-- Moved `OrbitalElements` struct from test file to `src/bodies.h`
-- Moved `calculate_orbital_elements()` function to `src/bodies.cpp`
+- Moved `OrbitalElements` struct from test file to `src/simulation.h`
+- Moved `calculate_orbital_elements()` function to `src/simulation.cpp`
 - Created `tests/test_comet_orbit.cpp` with SOI transition tracking
 - Added `ParentChange` tracking to detect when comet switches gravitational parent
 - Removed `tests/test_soi_transitions.cpp` (not useful for detecting actual transitions)
diff --git a/src/config_loader.cpp b/src/config_loader.cpp
index cb7a9d5..749fe17 100644
--- a/src/config_loader.cpp
+++ b/src/config_loader.cpp
@@ -91,14 +91,6 @@ bool parse_toml_body(toml_datum_t body_table, CelestialBody* body) {
     return true;
 }
 
-struct OrbitParams {
-    double eccentricity;
-    double semi_major_axis;
-};
-
-// Forward declaration
-void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params);
-
 // Load system configuration from TOML file
 bool load_system_config(SimulationState* sim, const char* filepath) {
     toml_result_t result = toml_parse_file_ex(filepath);
@@ -131,9 +123,6 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
         return false;
     }
 
-    // Parse each body and store orbit params
-    OrbitParams orbit_params[100];
-
     for (int i = 0; i < body_count; i++) {
         toml_datum_t body_table = bodies.u.arr.elem[i];
         if (!parse_toml_body(body_table, &sim->bodies[i])) {
@@ -141,189 +130,14 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
             toml_free(result);
             return false;
         }
-
-        // Store orbit parameters for velocity calculation
-        orbit_params[i].eccentricity = sim->bodies[i].eccentricity;
-        orbit_params[i].semi_major_axis = sim->bodies[i].semi_major_axis;
     }
 
     sim->body_count = body_count;
     toml_free(result);
 
-    // Calculate initial velocities
-    calculate_velocities(sim, orbit_params);
-
-    // Calculate SOI radii
+    calculate_initial_velocities(sim);
     calculate_soi_radii(sim);
 
     printf("Loaded %d bodies from %s\n", body_count, filepath);
     return true;
 }
-
-// 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;
-    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
-            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 bodies (already handled)
-        if (body->parent_index == -1) {
-            continue;
-        }
-
-        // Get parent body
-        if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
-
-            // Calculate relative position
-            Vec3 r = vec3_sub(body->position, parent->position);
-            double distance = vec3_magnitude(r);
-
-            if (distance < 1.0) {
-                body->velocity = {0.0, 0.0, 0.0};
-                continue;
-            }
-
-            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
-            // Cross product of r with z-axis gives perpendicular vector in XY plane
-            Vec3 z_axis = {0.0, 0.0, 1.0};
-
-            // Calculate cross product: r x z_axis
-            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);
-
-            // Add parent's velocity for absolute reference frame
-            body->velocity = vec3_add(body->velocity, parent->velocity);
-        }
-    }
-}
-
-// Calculate SOI radii for all bodies
-void calculate_soi_radii(SimulationState* sim) {
-    for (int i = 0; i < sim->body_count; i++) {
-        CelestialBody* body = &sim->bodies[i];
-
-        if (body->parent_index == -1) {
-            // Root body has very large SOI
-            body->soi_radius = 1e15;  // ~1000 AU
-        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
-
-            // Update SOI using Hill sphere approximation
-            update_soi(body, parent, body->semi_major_axis);
-        }
-    }
-}
diff --git a/src/config_loader.h b/src/config_loader.h
index ee3aecd..ffbb3e0 100644
--- a/src/config_loader.h
+++ b/src/config_loader.h
@@ -1,7 +1,7 @@
 #ifndef CONFIG_LOADER_H
 #define CONFIG_LOADER_H
 
-#include "bodies.h"
+#include "simulation.h"
 #include "../ext/tomlc17/src/tomlc17.h"
 
 // Load a system configuration from a file (TOML format)
@@ -12,10 +12,4 @@ bool extract_vec3_from_table(toml_datum_t table, Vec3* pos);
 bool extract_color_from_table(toml_datum_t table, float* color);
 bool parse_toml_body(toml_datum_t body_table, CelestialBody* body);
 
-// Calculate initial circular orbit velocities for all bodies
-void calculate_initial_velocities(SimulationState* sim);
-
-// Calculate SOI radii for all bodies
-void calculate_soi_radii(SimulationState* sim);
-
 #endif
diff --git a/src/main.cpp b/src/main.cpp
index 7071925..57b656e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,5 +1,5 @@
 #include "physics.h"
-#include "bodies.h"
+#include "simulation.h"
 #include "config_loader.h"
 #include "renderer.h"
 #include <cstdio>
diff --git a/src/physics.cpp b/src/physics.cpp
index 6e85768..0fd45e1 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -1,5 +1,5 @@
 #include "physics.h"
-#include "bodies.h"
+#include "simulation.h"
 #include <cmath>
 
 // Vector addition
@@ -12,6 +12,15 @@ Vec3 vec3_sub(Vec3 a, Vec3 b) {
     return {a.x - b.x, a.y - b.y, a.z - b.z};
 }
 
+// Cross product
+Vec3 vec3_cross(Vec3 a, Vec3 b) {
+    return {
+        a.y * b.z - a.z * b.y,
+        a.z * b.x - a.x * b.z,
+        a.x * b.y - a.y * b.x
+    };
+}
+
 // Scalar multiplication
 Vec3 vec3_scale(Vec3 v, double s) {
     return {v.x * s, v.y * s, v.z * s};
diff --git a/src/physics.h b/src/physics.h
index 2656917..eb8ecc9 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -15,6 +15,7 @@ const double G = 6.67430e-11;
 // Vector math functions
 Vec3 vec3_add(Vec3 a, Vec3 b);
 Vec3 vec3_sub(Vec3 a, Vec3 b);
+Vec3 vec3_cross(Vec3 a, Vec3 b);
 Vec3 vec3_scale(Vec3 v, double s);
 double vec3_magnitude(Vec3 v);
 double vec3_distance(Vec3 a, Vec3 b);
diff --git a/src/renderer.h b/src/renderer.h
index 7f1ed73..7519c4f 100644
--- a/src/renderer.h
+++ b/src/renderer.h
@@ -1,7 +1,7 @@
 #ifndef RENDERER_H
 #define RENDERER_H
 
-#include "bodies.h"
+#include "simulation.h"
 #include "raylib.h"
 
 // Rendering state
diff --git a/src/bodies.cpp b/src/simulation.cpp
similarity index 52%
rename from src/bodies.cpp
rename to src/simulation.cpp
index 2446cc0..63bbff1 100644
--- a/src/bodies.cpp
+++ b/src/simulation.cpp
@@ -1,4 +1,4 @@
-#include "bodies.h"
+#include "simulation.h"
 #include <cstdlib>
 #include <cstring>
 #include <cstdio>
@@ -139,6 +139,170 @@ 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);
+    double distance = vec3_magnitude(r);
+
+    if (distance < 1.0) {
+        body->velocity = {0.0, 0.0, 0.0};
+        return;
+    }
+
+    double e = body->eccentricity;
+    double a = body->semi_major_axis;
+    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);
+    }
+
+    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);
+    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];
+
+        if (body->parent_index == -1) {
+            continue;
+        }
+
+        if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
+            CelestialBody* parent = &sim->bodies[body->parent_index];
+            compute_orbital_velocity_from_vis_viva(body, parent);
+        }
+    }
+}
+
+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 compute_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};
+    }
+
+    set_child_bodies_velocity(sim);
+}
+
+void calculate_initial_velocities(SimulationState* sim) {
+    compute_initial_velocities(sim);
+}
+
+// Calculate SOI radii for all bodies
+void calculate_soi_radii(SimulationState* sim) {
+    for (int i = 0; i < sim->body_count; i++) {
+        CelestialBody* body = &sim->bodies[i];
+
+        if (body->parent_index == -1) {
+            body->soi_radius = 1e15;
+        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
+            CelestialBody* parent = &sim->bodies[body->parent_index];
+            update_soi(body, parent, body->semi_major_axis);
+        }
+    }
+}
+
 OrbitalElements calculate_orbital_elements(CelestialBody* body, CelestialBody* primary,
                                            CelestialBody* optional_ref_body, double current_time) {
     const double AU = 1.496e11;
diff --git a/src/bodies.h b/src/simulation.h
similarity index 90%
rename from src/bodies.h
rename to src/simulation.h
index d97e83a..739680f 100644
--- a/src/bodies.h
+++ b/src/simulation.h
@@ -38,6 +38,13 @@ int find_dominant_body(SimulationState* sim, int body_index);
 void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_axis);
 void update_simulation(SimulationState* sim);
 
+// Velocity initialization
+void compute_initial_velocities(SimulationState* sim);
+void calculate_initial_velocities(SimulationState* sim);
+
+// SOI helpers
+void calculate_soi_radii(SimulationState* sim);
+
 // Orbital elements calculation
 struct OrbitalElements {
     double time_days;
diff --git a/src/test_utilities.h b/src/test_utilities.h
index d351c72..cd532d5 100644
--- a/src/test_utilities.h
+++ b/src/test_utilities.h
@@ -1,7 +1,7 @@
 #ifndef TEST_UTILITIES_H
 #define TEST_UTILITIES_H
 
-#include "bodies.h"
+#include "simulation.h"
 #include "physics.h"
 
 struct OrbitalMetrics {
diff --git a/tests/test_comet_orbit.cpp b/tests/test_comet_orbit.cpp
index fd27c73..46e593b 100644
--- a/tests/test_comet_orbit.cpp
+++ b/tests/test_comet_orbit.cpp
@@ -1,6 +1,6 @@
 #include <catch2/catch_test_macros.hpp>
 #include "../src/physics.h"
-#include "../src/bodies.h"
+#include "../src/simulation.h"
 #include "../src/config_loader.h"
 #include "../src/test_utilities.h"
 #include <cmath>
diff --git a/tests/test_energy.cpp b/tests/test_energy.cpp
index 19e9fae..a773374 100644
--- a/tests/test_energy.cpp
+++ b/tests/test_energy.cpp
@@ -1,6 +1,6 @@
 #include <catch2/catch_test_macros.hpp>
 #include "../src/physics.h"
-#include "../src/bodies.h"
+#include "../src/simulation.h"
 #include "../src/config_loader.h"
 #include "../src/test_utilities.h"
 #include <cmath>
diff --git a/tests/test_orbital_period.cpp b/tests/test_orbital_period.cpp
index fa491c6..7ebf2f7 100644
--- a/tests/test_orbital_period.cpp
+++ b/tests/test_orbital_period.cpp
@@ -1,6 +1,6 @@
 #include <catch2/catch_test_macros.hpp>
 #include "../src/physics.h"
-#include "../src/bodies.h"
+#include "../src/simulation.h"
 #include "../src/config_loader.h"
 #include "../src/test_utilities.h"
 #include <cmath>