From c3e5a77de8d8d034eff689af2e4f5350474d6a54 Mon Sep 17 00:00:00 2001
From: "armindamon.riess" <armindamon.riess@uzh.ch>
Date: Wed, 14 Dec 2022 10:49:21 +0100
Subject: [PATCH] added function to get center of mass
---
lib/nBodySim.cpp | 28 +++++++++++++++++++++++-----
lib/nBodySim.hpp | 2 ++
lib/node.cpp | 4 ++++
lib/node.hpp | 2 ++
lib/tree.cpp | 4 ++++
lib/tree.hpp | 2 ++
6 files changed, 37 insertions(+), 5 deletions(-)
diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 9088f79..03243f6 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -68,25 +68,27 @@ void nBodySim::runSimulation(double dt, unsigned nSteps) {
for (unsigned i=0; i < nSteps; ++i) {
write2file(positions_, file, 3, skip);
+ Vector COM = getCenterOfMass();
// Integrate with leapfrog
auto start = std::chrono::high_resolution_clock::now();
kick(dt/2.0);
drift(dt);
- updateForces();
+ treeUpdateForces();
kick(dt/2.0);
auto stop = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
std::cout << std::right << "Step " << std::setw(5) << i
- << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude()
+ // << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude()
// << ", DEV = " << std::setw(12) << verifyForces()
- << ", MCD = " << std::setw(8) << calculateMeanCenterDistance()
- // << ", MID = " << std::setw(8) << calculateMeanCenterDistance()
+ // << ", MCD = " << std::setw(8) << calculateMeanCenterDistance()
+ << ", MID = " << std::setw(8) << calculateMeanInterparticleDistance()
+ << ", COM = " << std::setw(4) << "(" << COM[0] << ", " << COM[1] << ", " << COM[2] << ") "
<< ", E = " << std::setw(12) << calculateTotalEnergy()
<< ", T = " << std::setw(8) << duration.count()/1000.0 << " s"
<< std::endl;
- // updateTree();
+ updateTree();
}
file.close();
}
@@ -289,6 +291,22 @@ double nBodySim::calculateTotalEnergy() const {
return totalEnergy;
}
+Vector nBodySim::getCenterOfMass() const {
+ /*
+ Vector com = {0.0, 0.0, 0.0};
+ for (unsigned i = 0; i < nParticles_; i++) {
+ com[0] += positions_[3*i];
+ com[1] += positions_[3*i+1];
+ com[2] += positions_[3*i+2];
+ }
+ com[0] /= nParticles_;
+ com[1] /= nParticles_;
+ com[2] /= nParticles_;
+ return com;
+ */
+ return tree_->getCOM();
+}
+
void nBodySim::write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const {
// writes array to file, with N rows and dim columns
for (unsigned i = 0; i < nParticles_; i+=skip) {
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index fdf4120..daef913 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -36,6 +36,8 @@ public:
double calculateMeanCenterDistance() const;
// Loop over all particles and calculate the total energy
double calculateTotalEnergy() const;
+ // get center of mass
+ Vector getCenterOfMass() const;
// Writes data to file, every row is a particle
void write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const;
diff --git a/lib/node.cpp b/lib/node.cpp
index b1bf6c6..22ebc3c 100644
--- a/lib/node.cpp
+++ b/lib/node.cpp
@@ -254,4 +254,8 @@ unsigned Node::getOctant(unsigned particle) {
}
}
}
+}
+
+Vector Node::getCOM() const {
+ return centerOfMass_;
}
\ No newline at end of file
diff --git a/lib/node.hpp b/lib/node.hpp
index 244e530..b1ab8c2 100644
--- a/lib/node.hpp
+++ b/lib/node.hpp
@@ -32,6 +32,8 @@ public:
// Function that takes a particle index in the global particles_ array and returns the octant it is in relative to the center fo the current node
// this function doesn't check if the particle is actually in the node.
unsigned getOctant(unsigned particle);
+ // get center of mass
+ Vector getCOM() const;
private:
// Node* root_;
diff --git a/lib/tree.cpp b/lib/tree.cpp
index 3d20168..04bd078 100644
--- a/lib/tree.cpp
+++ b/lib/tree.cpp
@@ -31,4 +31,8 @@ void Tree::updateForce(unsigned particle) {
void Tree::update() {
root_->update(allParticles_, nParticles_);
+}
+
+Vector Tree::getCOM() const {
+ return root_->getCOM();
}
\ No newline at end of file
diff --git a/lib/tree.hpp b/lib/tree.hpp
index 86cac25..90cf1a8 100644
--- a/lib/tree.hpp
+++ b/lib/tree.hpp
@@ -16,6 +16,8 @@ public:
void updateForce(unsigned particle);
// Update tree
void update();
+ // get center of mass
+ Vector getCOM() const;
private:
Node* root_;
--
GitLab