minor formatting, function for center of mass

lib/nBodySim.hpp

@@ -10,20 +10,25 @@ public:
     nBodySim(std::string datafile);
     // destructor deletes arrays 
     ~nBodySim();
-    // run simulation
+
+    // run simulation for nSteps time steps of size dt, integration with leapfrog algorithm
     void runSimulation(double dt, unsigned nSteps);
+
     // loops over all pairs of particles and calculates forces, calculates current mean interparticle distance 
     void calculateForces();
-    // use treecode to calculate forces
+    // use treecode and multipole expansion to calculate forces
     void treeCalculateForces();
+
     // loop over all pairs of particles and calculate mean interparticle distance
     double calculateMeanInterparticleDistance();
     // loop over all force vectors and calculate mean force magnitude
     double calculateMeanForceMagnitude();
+
     // integrate using leapfrog algorithm. doesn't update forces, doesn't update tree
     void doTimeStep(double dt);
     // update tree
     void updateTree();
+
     // writes data to file, every row is a particle
     void write2file(const double* array, std::string filename, unsigned dim) const;
     // write current state of simulation to file
@@ -34,10 +39,13 @@ public:
     double* getPositions() { return positions_; }
     double* getForces() { return forces_; }
     unsigned getNParticles() const { return nParticles_; }
+    
     double getMeanInterparticleDistance() const { return meanInterparticleDistance_; }
     double getMeanForceMagnitude() const { return meanForceMagnitude_; }
+
 private:
     Tree* tree_;
+
     unsigned nParticles_;
     // 3d vectors are stored as a array of length 3*nParticles_ in the format x1, y1, z1, x2, y2, z2, ...
     double* masses_;
@@ -46,6 +54,7 @@ private:
     double softening_;
     double* potential_;
     double* forces_;
+
     double meanInterparticleDistance_;
     double meanForceMagnitude_;
 };

lib/node.cpp

@@ -104,4 +104,13 @@ double* Node::calculateForce(unsigned particle) {
         // TODO
     }
     return nullptr;
+}
+
+void Node::calculateCenterOfMass() {
+    for (unsigned j=0; j<3; ++j)
+        centerOfMass_[j] = 0;
+    
+    for (unsigned i : localParticles_)
+        for (unsigned j=0; j<3; ++j)
+            centerOfMass_[j] += particles_[3*i + j] / nLocalParticles_;
 }
\ No newline at end of file

lib/node.hpp

@@ -13,7 +13,7 @@ public:
         double size, double center[3], std::vector<unsigned>& localParticles, unsigned nLocalParticles);
     // destuctor
     ~Node();
-    // copy constructor
+    // copy constructor (default)
     Node(const Node&) = default;
     // copy assignment
     Node& operator=(const Node&);
@@ -21,18 +21,26 @@ public:
     void update(std::vector<unsigned>& allParticles);
     // calculate force on a particle recursively
     double* calculateForce(unsigned particle);
+    // calculate center of mass
+    void calculateCenterOfMass();
+
 private:
     Node* root_;
     Node* parent_;
     Node* children_[8];
+
     double* masses_;
     double* particles_;
     unsigned nParticles_;
+
     unsigned depth_;
     double size_;
     double center_[3];
+    double centerOfMass_[3];
+
     std::vector<unsigned> localParticles_;
     unsigned nLocalParticles_;
+
     double mass_;
     const double theta0_ = 0.4;
     bool isLeaf_;

lib/tree.hpp

@@ -16,14 +16,18 @@ public:
     void calculateForce(unsigned particle);
     // update tree
     void update();
+
 private:
     Node* root_;
+
     double* masses_;
     double* particles_;
     double* forces_;
     unsigned nParticles_;
+
     double size_;
     double center_[3];
+    
     std::vector<unsigned> allParticles_;
 };