preparing for implementation of full simulation

lib/nBodySim.cpp

@@ -8,40 +8,40 @@
 nBodySim::nBodySim(std::string& datafile) {
     std::string dummy;
     std::ifstream file(datafile);
-    file >> dummy >> dummy >> N_;
+    file >> dummy >> dummy >> nParticles_;
 
-    masses_ = new double[N_];
-    for (unsigned i=0; i < N_; ++i) {
+    masses_ = new double[nParticles_];
+    for (unsigned i=0; i < nParticles_; ++i) {
         file >> masses_[i];
     }
 
-    positions_ = new double[3*N_];
+    positions_ = new double[3*nParticles_];
     for (unsigned i=0; i < 3; ++i) {
-        for (unsigned j=0; j < N_; ++j) {
+        for (unsigned j=0; j < nParticles_; ++j) {
             file >> positions_[i+3*j];
         }
     }
     
-    velocities_ = new double[3*N_];
+    velocities_ = new double[3*nParticles_];
     for (unsigned i=0; i < 3; ++i) {
-        for (unsigned j=0; j < N_; ++j) {
+        for (unsigned j=0; j < nParticles_; ++j) {
             file >> velocities_[i+3*j];
         }
     }
 
-    softening_ = new double[N_];
-    for (unsigned i=0; i < N_; ++i) {
+    softening_ = new double[nParticles_];
+    for (unsigned i=0; i < nParticles_; ++i) {
         file >> softening_[i];
     }
     
-    potential_ = new double[N_];
-    for (unsigned i=0; i < N_; ++i) {
+    potential_ = new double[nParticles_];
+    for (unsigned i=0; i < nParticles_; ++i) {
         file >> potential_[i];
     }
 
-    forces_ = new double[3*N_];
+    forces_ = new double[3*nParticles_];
     #pragma omp parallel for
-    for (unsigned i=0; i < 3*N_; ++i) {
+    for (unsigned i=0; i < 3*nParticles_; ++i) {
         forces_[i] = 0.0;
     }
 }
@@ -55,15 +55,22 @@ nBodySim::~nBodySim() {
     delete[] forces_;
 }
 
+void nBodySim::runSimulation(double dt, unsigned nSteps) {
+    for (unsigned i=0; i < nSteps; ++i) {
+        calculateForces();
+        updatePositions(dt);
+    }
+}
+
 void nBodySim::calculateForces() {
     #pragma omp parallel for
-    for (unsigned i=0; i < 3*N_; ++i) {
+    for (unsigned i=0; i < 3*nParticles_; ++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++) {
+    for (unsigned i = 0; i < nParticles_; i++) {
+        for (unsigned j = 0; j < nParticles_; j++) {
             if (i == j) continue;
             // calculate distance between particles
             double dx = positions_[3*j]   - positions_[3*i];
@@ -85,10 +92,20 @@ void nBodySim::calculateForces() {
     }
 }
 
+void nBodySim::updatePositions(double dt) {
+    #pragma omp parallel for
+    for (unsigned i=0; i < nParticles_; ++i) {
+        for (unsigned j=0; j < 3; ++j) {
+            // update velocity
+            // update position
+        }
+    }
+}
+
 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 < nParticles_; i++) {
         for (unsigned j = 0; j < dim; j++) {
             file << std::right << std::setw(16) << array[dim*i+j] << " ";
         }

lib/nBodySim.hpp

@@ -9,17 +9,23 @@ public:
     nBodySim(std::string& datafile);
     // destructor deletes arrays 
     ~nBodySim();
+    // run simulation
+    void runSimulation(double dt, unsigned nSteps);
     // loops over all pairs of particles and calculates forces
     void calculateForces();
+    // updates positions and velocities using the calculated forces
+    void updatePositions(double dt);
+    // writes data to file, every row is a particle
+    void write2file(const double* array, std::string filename, unsigned dim) const;
+
+    // functions to get private variables
     double* getMasses() { return masses_; }
     double* getPositions() { return positions_; }
     double* getForces() { return forces_; }
-    // writes data to file, every row is a particle
-    void write2file(const double* array, std::string filename, unsigned dim) const;
-    unsigned getN() const { return N_; }
+    unsigned getNParticles() const { return nParticles_; }
 private:
-    unsigned N_;
-    // 3d vectors are stored as a 1d array of length 3*N_ in the format x1, y1, z1, x2, y2, z2, ...
+    unsigned nParticles_;
+    // 3d vectors are stored as a 1d array of length 3*nParticles_ in the format x1, y1, z1, x2, y2, z2, ...
     double* masses_;
     double* positions_;
     double* velocities_;

main.cpp

@@ -18,8 +18,11 @@ int main(int argc, char** argv) {
             return 1;
     }
 
+    double dt = 0.01;
+    unsigned nSteps = 1000;
+
     nBodySim sim(datafile);
-    unsigned N = sim.getN();
+    unsigned N = sim.getNParticles();
     std::cout << "N = " << N << std::endl;
 
     /*