Revert BH kernel

Author: TimGoTheCreator

src/gravity/octree.h

@@ -2,6 +2,7 @@
 #include "../struct/particle.h"
 #include <vector>
 
+
 // =========================
 // Octree definition
 // =========================
@@ -89,12 +90,6 @@ struct Octree {
     }
 };
 
-inline float rsqrt(float x) {
-    float y = _mm_cvtss_f32(_mm_rsqrt_ss(_mm_set_ss(x)));
-    // optional refinement:
-    return y * (1.5f - 0.5f * x * y * y);
-}
-
 // =========================
 // Barnes–Hut force function
 // =========================
@@ -106,25 +101,21 @@ inline void bhForce(Octree* node, Particle& p, real theta, real dt)
     real dy = node->cy - p.y;
     real dz = node->cz - p.z;
 
-    real r2 = dx*dx + dy*dy + dz*dz + real(1e-6);
+    real dist = sqrt(dx*dx + dy*dy + dz*dz) + 1e-6;
 
-    // Barnes–Hut acceptance
-    real dist = sqrt(r2);
+    // Barnes–Hut acceptance criterion
     if (node->leaf || (node->size / dist) < theta) {
-
-        // division-free inverse distance cubed
-        real invDist = rsqrt(r2);
-        real invDist3 = invDist * invDist * invDist;
-
-        real f = node->m * invDist3 * dt;
-
+        real inv = 1.0 / (dist * dist * dist);
+        real f = node->m * inv * dt;
         p.vx += dx * f;
         p.vy += dy * f;
         p.vz += dz * f;
         return;
     }
 
-    for (int i = 0; i < 8; i++)
+    // Otherwise recurse into children
+    for (int i = 0; i < 8; i++) {
         if (node->child[i])
             bhForce(node->child[i], p, theta, dt);
-}
+    }
+}