diff --git a/src/config_loader.cpp b/src/config_loader.cpp
index efa90fc..d16a96e 100644
--- a/src/config_loader.cpp
+++ b/src/config_loader.cpp
@@ -132,21 +132,21 @@ void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params) {
             // v = sqrt(G * M_other / r)
             double speed = sqrt(G * other_mass / distance);
 
-            // Create velocity perpendicular to position vector (counter-clockwise from +Z)
+            // Create velocity perpendicular to position vector
             Vec3 z_axis = {0.0, 0.0, 1.0};
             Vec3 vel_dir = {
-                z_axis.y * r.z - z_axis.z * r.y,
-                z_axis.z * r.x - z_axis.x * r.z,
-                z_axis.x * r.y - z_axis.y * r.x
+                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 = x_axis.y * r.z - x_axis.z * r.y;
-                vel_dir.y = x_axis.z * r.x - x_axis.x * r.z;
-                vel_dir.z = x_axis.x * r.y - x_axis.y * r.x;
+                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
@@ -197,23 +197,23 @@ void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params) {
 
             // Create velocity perpendicular to position vector
             // If position is mostly in XY plane, make velocity in XY plane
-            // Cross product z_axis x r gives counter-clockwise orbit when viewed from +Z
+            // 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: z_axis x r
+            // Calculate cross product: r x z_axis
             Vec3 vel_dir = {
-                z_axis.y * r.z - z_axis.z * r.y,
-                z_axis.z * r.x - z_axis.x * r.z,
-                z_axis.x * r.y - z_axis.y * r.x
+                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 = x_axis.y * r.z - x_axis.z * r.y;
-                vel_dir.y = x_axis.z * r.x - x_axis.x * r.z;
-                vel_dir.z = x_axis.x * r.y - x_axis.y * r.x;
+                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