corrected tree code force: set to 0, use softening

lib/nBodySim.cpp

@@ -129,7 +129,11 @@ void nBodySim::updateForces() {
 }
 
 void nBodySim::treeUpdateForces() {
-    // #pragma omp parallel for
+    #pragma omp parallel for
+    for (unsigned i=0; i < 3*nParticles_; ++i) {
+        forces_[i] = 0.0;
+    }
+    #pragma omp parallel for
     for (unsigned i=0; i < nParticles_; ++i) {
         tree_->updateForce(i);
     }

lib/node.cpp

@@ -127,12 +127,12 @@ Vector Node::calculateForce(unsigned particle) {
                     dy /= r;
                     dz /= r;
 
-                    double m = mass_;
+                    double factor = mass_*mass_ / (r2 + softening_*softening_);
 
                     // add the force from each particle in the child
-                    force[0] += m*m * dx / (r2 + softening_*softening_);
-                    force[1] += m*m * dy / (r2 + softening_*softening_);
-                    force[2] += m*m * dz / (r2 + softening_*softening_);
+                    force[0] += factor * dx;
+                    force[1] += factor * dy;
+                    force[2] += factor * dz;
                 }
             // if the child exists, get the force from the child
             } else {
@@ -148,10 +148,11 @@ Vector Node::calculateForce(unsigned particle) {
     } else {
         
         // calculate force from this node (from the center of mass)
-        double r3 = lambda*lambda*lambda;
+        double r2 = lambda*lambda;
 
         for (unsigned j=0; j<3; ++j) {
-            force[j] += totalMass_ * (centerOfMass_[j] - positions_[3*particle + j]) / r3;
+            double factor = totalMass_*mass_ / (r2 + softening_*softening_);
+            force[j] += factor * (centerOfMass_[j] - positions_[3*particle + j])/lambda;
         }
 
     }

lib/node.hpp

@@ -51,7 +51,7 @@ private:
     double totalMass_;
     const double theta0_ = 0.4;
     const double minNParticles_ = 8;
-    const double maxDepth_ = 24;
+    const double maxDepth_ = 16;
     bool isLeaf_;
 };