As part of my training, I implemented a n-body class in C++ to simulate gravitational interaction of bodies and to get more familiar with features that C++ offers such as object oriented programming.

This implementation uses a direct integration (Verlet integration) of the differential equations which results in a time complexity of \$\mathcal{O}(n^2)\$, where \$n\$ is the number of particles.

Please be as hard as possible with this implementation and give me constructive feedback. 

I would appreciate advice especially in the following areas:

1. Code style (readability, naming conventions)
2. Class design
3. Efficieny (how to avoid unnecessary complexity)
4. Reinventing the wheel (does the STL offer functionality I should use in my code?)
5. Memory usage

**main.cpp**

```
#include "nbody.h"

int main(int argc, char* argv[]) {
    Nbody nbody(16, 0.001, 1);
    nbody.timeIntegration();
    return 0;
}

```

**nbody.h**

```
#ifndef NBODY_H
#define NBODY_H

// Parameters
const int DIM = 2;          // dimensions
const double EPS = 1e-4;    // smoothing parameter

// Function prototypes
inline double sqr(double);

struct Particle{
    double m;           // mass
    double x[DIM];      // position
    double v[DIM];      // velocity
    double F[DIM];      // force 
    double F_old[DIM];  // force past time step
};

// Nbody class 
class Nbody {
    private:
        int step = 0;
        double t = 0;
        const int n;         // number of particles
        const double dt;           // step size
        const double t_max;        // max simulation time
        Particle *p = new Particle[n]; // allocate memory 
        void init_data();
    public:
        ~Nbody();
        Nbody(int n_, double dt_, double t_max_);
        inline void print_parameter() const;
        inline void print_data() const;
        inline void write_data(int step) const;
        void timeIntegration();
        void comp_force();
        void force(Particle*, Particle*);
        void comp_position();
        void comp_velocity();
        void update_position(Particle*);
        void update_velocity(Particle*);
};

#endif

```

**nbody.cpp**

```
#include <iostream>
#include <fstream> 
#include <cmath>
#include <random>

#include "nbody.h"

// Class methods
Nbody::Nbody(int n_, double dt_, double t_max_) : n(n_), dt(dt_), t_max(t_max_) {
    init_data();
}

Nbody::~Nbody() {
    delete[] p; 
    p = 0; 
}

void Nbody::timeIntegration() {
    comp_force();
    for(; t<t_max; t+=dt, step+=1) {
        comp_position();
        comp_force();
        comp_velocity();
        if (step % 10 == 0) {
            write_data(step);
            //print_data();
        }
    }
}

void Nbody::update_velocity(Particle *p) {
    double a = dt * 0.5 / p->m;
    for (int d=0; d<DIM; d++) {
        p->v[d] += a * (p->F[d] + p->F_old[d]);
    }
}

void Nbody::update_position(Particle *p) {
    double a = dt * 0.5 / p->m;
    for (int d=0; d<DIM; d++) {
        p->x[d] += dt * (p->v[d] + a * p->F[d]);
        p->F_old[d] = p->F[d];
    }
}

void Nbody::comp_velocity() {
    for (int i=0; i<n; i++) {
        update_velocity(&p[i]);
    }
}

void Nbody::comp_position() {
    for (int i=0; i<n; i++) {
        update_position(&p[i]);
    }
}

void Nbody::comp_force() {
    for (int i=0; i<n; i++) {
        for (int d=0; d<DIM; d++) {
            p[i].F[d] = 0;
        }
    }
    for (int i=0; i<n; i++) {
        for (int j=i+1; j<n; j++) {
            force(&p[i], &p[j]);
        }
    }
}

void Nbody::force(Particle *i, Particle *j) {
    double r=EPS; // smoothing
    for (int d=0; d<DIM; d++) {
        r += sqr(j->x[d] - i->x[d]);
    }
    double f = i->m * j->m / (sqrt(r) * r);
    for (int d=0; d<DIM; d++) {
        i->F[d] += f * (j->x[d] - i->x[d]);
        j->F[d] -= f * (j->x[d] - i->x[d]);
    }
}

void Nbody::write_data(int step) const {
    std::ofstream results;
    std::string file_name = "data_" + std::to_string(step) + ".log";
    results.open(file_name);
    if (results.fail()) { // or (!results) ?
        std::cerr << "Error\n" << std::endl;
    } else {
        for (int i=0; i<n; i++) {
            results << t << " ";
            results << p[i].m << " ";
            for (int d=0; d<DIM; d++) {
                results << p[i].x[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << p[i].v[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << p[i].F[d] << " ";
            }
            results << std::endl;
        }
        results.close();
    }
}

void Nbody::print_data() const {
    std::cout.setf(std::ios_base::scientific);
    std::cout.precision(5);
    for (int i=0; i<n; i++) {
        std::cout << t << " ";
        std::cout << p[i].m << " ";
        for (int d=0; d<DIM; d++) {
            std::cout << p[i].x[d] << " ";
        }
        for (int d=0; d<DIM; d++) {
            std::cout << p[i].v[d] << " ";
        }
        for (int d=0; d<DIM; d++) {
            std::cout << p[i].F[d] << " ";
        }
        std::cout << std::endl;
    }
}

void Nbody::init_data() {
    std::random_device rd;          
    std::mt19937 generator(rd()); 
    std::uniform_real_distribution<double> distribution_x(0.0,1.0);
    std::uniform_real_distribution<double> distribution_v(-1.0,1.0);
    for (int i=0; i<n; i++) {
        p[i].m = 1./n;
        for (int d=0; d<DIM; d++) {
            p[i].x[d] = distribution_x(generator);
            p[i].v[d] = distribution_v(generator);
            p[i].F[d] = 0.0;
            p[i].F_old[d] = 0.0;
        }
    }
}

inline void Nbody::print_parameter() const {
    std::cout << n << " " << dt << " " << t_max << std::endl;
}

// Other Functions

inline double sqr(double x) {
    return x*x;
}

```