Fixed Bugs

Author: TimGoTheCreator

src/gravity/step.h

@@ -11,26 +11,53 @@ inline void Step(std::vector<Particle>& p, real dt) {
     real theta = 0.5;
     real half  = dt * real(0.5);
 
+    auto buildTree = [&](Octree*& root) {
+        // Compute bounding box
+        real minx=+1e30, miny=+1e30, minz=+1e30;
+        real maxx=-1e30, maxy=-1e30, maxz=-1e30;
+
+        for (auto& a : p) {
+            minx = std::min(minx, a.x);
+            miny = std::min(miny, a.y);
+            minz = std::min(minz, a.z);
+            maxx = std::max(maxx, a.x);
+            maxy = std::max(maxy, a.y);
+            maxz = std::max(maxz, a.z);
+        }
+
+        real cx = (minx + maxx) * 0.5;
+        real cy = (miny + maxy) * 0.5;
+        real cz = (minz + maxz) * 0.5;
+        real size = std::max({maxx-minx, maxy-miny, maxz-minz}) * 0.5;
+
+        if (size <= 0) size = 1; // safety
+
+        root = new Octree(cx, cy, cz, size);
+
+        for (auto& a : p)
+            root->insert(&a);
+
+        root->computeMass();
+    };
+
     // =========================
     // First Kick (dt/2)
     // =========================
     {
-        Octree root(0,0,0, 1e9);
-
-        for (auto& a : p) {
-            root.insert(&a);
-        }
-        root.computeMass();
+        Octree* root = nullptr;
+        buildTree(root);
 
         #pragma omp parallel for schedule(static)
         for (int i = 0; i < (int)p.size(); i++) {
             real ax = 0, ay = 0, az = 0;
-            bhAccel(&root, p[i], theta, ax, ay, az);
+            bhAccel(root, p[i], theta, ax, ay, az);
 
             p[i].vx += ax * half;
             p[i].vy += ay * half;
             p[i].vz += az * half;
         }
+
+        delete root;
     }
 
     // =========================
@@ -47,21 +74,19 @@ inline void Step(std::vector<Particle>& p, real dt) {
     // Second Kick (dt/2)
     // =========================
     {
-        Octree root(0,0,0, 1e9);
-
-        for (auto& a : p) {
-            root.insert(&a);
-        }
-        root.computeMass();
+        Octree* root = nullptr;
+        buildTree(root);
 
         #pragma omp parallel for schedule(static)
         for (int i = 0; i < (int)p.size(); i++) {
             real ax = 0, ay = 0, az = 0;
-            bhAccel(&root, p[i], theta, ax, ay, az);
+            bhAccel(root, p[i], theta, ax, ay, az);
 
             p[i].vx += ax * half;
             p[i].vy += ay * half;
             p[i].vz += az * half;
         }
+
+        delete root;
     }
-}
+}