use maxDepth, no segfault but result incorrect

lib/node.cpp

@@ -89,6 +89,7 @@ Node::Node(double* masses, double* positions, unsigned nParticles, unsigned dept
 Node::~Node() {
     for (unsigned i=0; i<8; ++i) {
         delete children_[i];
+        children_[i] = nullptr;
     }
 }
 
@@ -106,6 +107,9 @@ Node& Node::operator=(const Node& other) {
 }
 
 void Node::update(List localParticles, unsigned nLocalParticles) {
+    double minNParticles = 8;
+    double maxDepth = 10;
+
     localParticles_ = localParticles;
     nLocalParticles_ = nLocalParticles;
     updateMassAndCOM();
@@ -118,11 +122,19 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
     }
 
     isLeaf_ = true;
+    if (nLocalParticles_ <= minNParticles || depth_ >= maxDepth) {
+        // if the node contains less than minNParticles, it is a leaf
+        for (unsigned i=0; i<8; ++i) {
+            delete children_[i];
+            children_[i] = nullptr;
+        }
+        return;
+    }
     for (unsigned i=0; i<8; ++i) {
         if (children_[i] != nullptr) {
-            // if child is not a nullptr and particleDivision[i] is empty, delete the child
+            // if child is not a nullptr and particleDivision[i] is empty or contains a single particle, delete the child
             // otherwise update the child
-            if (particleDivision[i].size() == 0) {
+            if (particleDivision[i].size() <= minNParticles) {
                 delete children_[i];
                 children_[i] = nullptr;
             } else {
@@ -130,8 +142,8 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
                 isLeaf_ = false;
             }
         } else {
-            // if child is a nullptr and particleDivision[i] is not empty, create a new child
-            if (particleDivision[i].size() != 0) {
+            // if child is a nullptr and particleDivision[i] contains more than one particle, create a new child
+            if (particleDivision[i].size() > minNParticles) {
                 double childSize = size_/2;
                 Vector childCenter(3);
                 for (unsigned j=0; j<3; ++j) {