From 8c89deabb8e58626b66fe7fdbff717f52b0e2aa1 Mon Sep 17 00:00:00 2001
From: "armindamon.riess" <armindamon.riess@uzh.ch>
Date: Wed, 14 Dec 2022 14:40:23 +0100
Subject: [PATCH] clean up, increase number of steps
---
lib/nBodySim.cpp | 32 +++++++++++++-------------------
main.cpp | 2 +-
plots/animation.py | 2 +-
3 files changed, 15 insertions(+), 21 deletions(-)
diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 179a270..256152e 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -63,11 +63,18 @@ nBodySim::~nBodySim() {
void nBodySim::runSimulation(double dt, unsigned nSteps) {
const unsigned skip = 10;
+
std::ofstream positionsFile("../out/positions.dat");
+ std::ofstream energyFile("../out/energy.dat");
+
positionsFile << nParticles_ << " " << skip << " " << dt << std::endl;
+ energyFile << nParticles_ << " " << skip << " " << dt << std::endl;
+
for (unsigned i=0; i < nSteps; ++i) {
write2file(positions_, positionsFile, 3, skip);
+ energyFile << calculateTotalEnergy() << std::endl;
+
Vector COM = getCenterOfMass();
// Integrate with leapfrog
@@ -79,13 +86,13 @@ void nBodySim::runSimulation(double dt, unsigned nSteps) {
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()
- << ", DEV = " << std::setw(12) << verifyForces()
+ std::cout << std::right << "Step " << std::setw(5) << i
+ // << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude()
+ // << ", DEV = " << std::setw(12) << verifyForces()
// << ", MCD = " << std::setw(8) << calculateMeanCenterDistance()
- << ", MID = " << std::setw(8) << calculateMeanInterparticleDistance()
- << ", COM = " << "(" << COM[0] << ", " << COM[1] << ", " << COM[2] << ")"
- << ", E = " << std::setw(12) << calculateTotalEnergy()
+ // << ", MID = " << std::setw(8) << calculateMeanInterparticleDistance()
+ // << ", COM = " << "(" << COM[0] << ", " << COM[1] << ", " << COM[2] << ")"
+ // << ", E = " << std::setw(12) << calculateTotalEnergy()
<< ", T = " << std::setw(8) << duration.count()/1000.0 << " s"
<< std::endl;
updateTree();
@@ -292,18 +299,6 @@ double nBodySim::calculateTotalEnergy() const {
}
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();
}
@@ -344,7 +339,6 @@ void nBodySim::saveState2file(unsigned step, std::ofstream& file) const {
}
}
-
void nBodySim::saveTree2file(std::string filename) const {
std::ofstream file(filename);
tree_->save2file(file);
diff --git a/main.cpp b/main.cpp
index 0ab77d5..2d66c68 100644
--- a/main.cpp
+++ b/main.cpp
@@ -19,7 +19,7 @@ int main(int argc, char** argv) {
}
double dt = 0.000005;
- unsigned nSteps = 120;
+ unsigned nSteps = 240;
std::cout << "dt = " << dt << std::endl;
std::cout << "nSteps = " << nSteps << std::endl;
diff --git a/plots/animation.py b/plots/animation.py
index ad7ab85..f6dd47c 100644
--- a/plots/animation.py
+++ b/plots/animation.py
@@ -20,7 +20,7 @@ def animate(i):
yi = y[i*N:(i+1)*N]
plt.cla()
plt.scatter(xi, yi, s=0.1)
- plt.scatter(np.sum(xi)/N, np.sum(yi)/N, s=1, c="red")
+ # plt.scatter(np.sum(xi)/N, np.sum(yi)/N, s=1, c="red")
plt.xlabel("x")
plt.ylabel("y")
plt.axis("equal")
--
GitLab