correct force calculation, implement OpenMP

lib/nBodySim.cpp

 #include "nBodySim.hpp"
+#include <omp.h>
+#include <cmath>
 #include <string>
 #include <fstream>
 #include <iomanip>
@@ -17,14 +19,14 @@ nBodySim::nBodySim(std::string& datafile) {
     positions_ = new double[3*N_];
     for (unsigned i=0; i < 3; ++i) {
         for (unsigned j=0; j < N_; ++j) {
-            file >> positions_[i*N_+j];
+            file >> positions_[i+3*j];
         }
     }
     
     velocities_ = new double[3*N_];
     for (unsigned i=0; i < 3; ++i) {
         for (unsigned j=0; j < N_; ++j) {
-            file >> velocities_[i*N_+j];
+            file >> velocities_[i+3*j];
         }
     }
 
@@ -39,10 +41,8 @@ nBodySim::nBodySim(std::string& datafile) {
     }
 
     forces_ = new double[3*N_];
-    for (unsigned i=0; i < 3; ++i) {
-        for (unsigned j=0; j < N_; ++j) {
-            file >> forces_[i*N_+j];
-        }
+    for (unsigned i=0; i < 3*N_; ++i) {
+        forces_[i] = 0.0;
     }
 }
 
@@ -56,29 +56,41 @@ nBodySim::~nBodySim() {
 }
 
 void nBodySim::calculateForces() {
+    #pragma omp parallel for
+    for (unsigned i=0; i < 3*N_; ++i) {
+        forces_[i] = 0.0;
+    }
     // loop over each pair of particles and calculate the force between them
+    #pragma omp parallel for collapse(2)
     for (unsigned i = 0; i < N_; i++) {
         for (unsigned j = 0; j < N_; j++) {
             if (i == j) continue;
             // calculate distance between particles
-            double dx = positions_[j]   - positions_[i];
-            double dy = positions_[j+1] - positions_[i+1];
-            double dz = positions_[j+2] - positions_[i+2];
+            double dx = positions_[3*j]   - positions_[3*i];
+            double dy = positions_[3*j+1] - positions_[3*i+1];
+            double dz = positions_[3*j+2] - positions_[3*i+2];
             double r2 = dx*dx + dy*dy + dz*dz;
+            double r = std::sqrt(r2);
+            dx /= r;
+            dy /= r;
+            dz /= r;
             // calculate force
-            forces_[i]   += masses_[j] * dx / (r2 + softening_[i]);
-            forces_[i+1] += masses_[j] * dy / (r2 + softening_[i]);
-            forces_[i+2] += masses_[j] * dz / (r2 + softening_[i]);
+            double mi = masses_[i];
+            double mj = masses_[j];
+            double s = softening_[i];
+            forces_[3*i]   += mi*mj * dx / (r2 + s*s);
+            forces_[3*i+1] += mi*mj * dy / (r2 + s*s);
+            forces_[3*i+2] += mi*mj * dz / (r2 + s*s);
         }
     }
 }
 
-void nBodySim::write2file(const double* array, std::string filename, unsigned N, unsigned dim) {
+void nBodySim::write2file(const double* array, std::string filename, unsigned dim) const {
     // writes array to file, with N rows and dim columns
     std::ofstream file(filename);
-    for (unsigned i = 0; i < N; i++) {
+    for (unsigned i = 0; i < N_; i++) {
         for (unsigned j = 0; j < dim; j++) {
-            file << std::right << std::setw(16) << array[i+j] << " ";
+            file << std::right << std::setw(16) << array[dim*i+j] << " ";
         }
         file << std::endl;
     }

lib/nBodySim.hpp

@@ -10,9 +10,11 @@ public:
     ~nBodySim();
     // methods
     void calculateForces();
+    double* getMasses() { return masses_; }
+    double* getPositions() { return positions_; }
     double* getForces() { return forces_; }
-    void write2file(const double* array, std::string filename, unsigned N, unsigned dim);
-    unsigned getN() { return N_; }
+    void write2file(const double* array, std::string filename, unsigned dim) const;
+    unsigned getN() const { return N_; }
 private:
     unsigned N_;
     // 3d vectors are stored as a 1d array of length 3*N_ in the format x1, y1, z1, x2, y2, z2, ...