diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 179a270daa45609b4e4ea0f324ff6b3edda9d530..256152ec54eb488c0ab98eafe4f9b5d187c7554f 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 0ab77d565e4b97e1512e1e9aeaa087242ffabbff..2d66c680267f21a763d68da4f1e4cbc2272efd54 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 ad7ab8561e43cfd9094a2ae59db6bf56a26b99e3..f6dd47cdaa243f7c7cd2ad225b41469b352ae461 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")