Brute Force N Body Implementation in C++

I have wrote the following code in C++ for the n-body problem. This code is sequential as later on I am planning to parallelize it using OpenMP. I want to know whether I have wrote the correct implementation, if there are some errors or bugs, or if this is an invalid approach.

```
#include
#include
#include
using namespace std;
#define N 100
#define G 6.673e-11
#define TIMESTAMP 1e11
struct Particle{
double rx, ry;//position components
double vx, vy;//velocity components
double fx, fy;//force components
double mass;//mass of the particle

};
Particle Update(Particle p, double timestamp)
{
p.vx += timestamp*p.fx / p.mass;
p.vy += timestamp*p.fy / p.mass;
p.rx += timestamp*p.vx;
p.ry += timestamp*p.vy;
return p;
}
void PrintParticle(Particle p)
{
printf("rx == %f ry == %f vx == %f vy == %f mass == %f\n", p.rx,p.ry,p.vx,p.vy,p.mass);
}
//Reset the forces on particle
Particle ResetForce(Particle p)
{
p.fx = 0.0;
p.fy = 0.0;
return p;
}
//Add force to particle a by particle b
Particle AddForce(Particle a,Particle b)
{
double EPS = 3E4; // softening parameter (just to avoid infinities)
double dx = b.rx - a.rx;
double dy = b.ry - a.ry;
double dist = sqrt(dxdx + dydy);
double F = (G a.mass b.mass) / (distdist + EPSEPS);
a.fx += F * dx / dist;
a.fy += F * dy / dist;
return a;

}

int main()
{
Particle particles[N];
srand(time(NULL));
//randomly generating N Particles
for (int i = 0; i < N; i++){
double rx = 1e18exp(-1.8)(.5 - rand());
particles[i].rx = rx;
double ry = 1e18exp(-1.8)(.5 - rand());
particles[i].ry = ry;
double vx = 1e18exp(-1.8)(.5 - rand());
particles[i].vx = vx;
double vy = 1e18exp(-1.8)(.5 - rand());
particles[i].vy = vy;
double mass = 1.98892e30rand()10 + 1e20;
particles[i].mass = mass;

}

int numberofiterations = 10;
int count = 0;
whil

Here are some observations that may help you improve your code.

Make sure to #include all required headers

This program calls printf and srand but does not include the corresponding headers. Fix that by adding these lines:

#include 
#include 

Use objects

You have a Particle structure and then separate functions that operate on Particle data. With only a slight syntax change, you would have a real object instead of C-style code written in C++. For example, you have a loop right now that says this:

for (int i = 0; i < N; i++)
{
    particles[i] = Update(particles[i], TIMESTAMP);
}

With objects, and with C++11 or newer, you could instead write this:

for (auto &p : particles)
    p.update(TIMESTAMP);

All that would be required would be to make Update a member function.

Prefer stream I/O to printf

The printf function was a capable workhorse for many years, but the C++ iostream library is better in a number of regards. Although it's often more typing for the programmer initially, it's better because it has better type checking, less possibility for runtime overhead, and it fits well with the rest of C++. So I would probably change your PrintParticle to this:

friend ostream& operator<<(ostream& out, const Particle &p) {
    return out << 
        "rx == " << p.rx_ <<
        " ry == " << p.ry_ <<
        " vx == " << p.vx_ <<
        " vy == " << p.vy_ <<
        " mass == " << p.mass_ << '\n';
}

Use const where practical

There are a number of things in your code that look and act like constants, such as numberofiterations in main(). Declare them as const for a program more resistant to error. Also, move the constants, such as N into main rather than having them at file scope.

Create a function rather than repeating code

The functions which create the particle use the same rand()-based function several times. Instead of repeating it, make it into a function:

double xrand() {
    return 1e18*exp(-1.8)*(.5 - rand());
}

Don't abuse using namespace std

Putting using namespace std at the top of every program is a bad habit that you'd do well to avoid. It isn't necessarily wrong to use, but be aware of when you absolutely shouldn't do it (such as in header files).

Consider using a better random number generator

If you are using a compiler that supports at least C++11, consider using a better random number generator. In particular, instead of rand, you might want to look at std::uniform_real_distribution and friends in the ` header.

Consider altering your algorithm

At the moment, the force on each particle from every other particle is calculated exhaustively, but Newton tells us that the force actually acts equally (but opposite) on both particles, so your loop could be simplified and calculations speeded if you take advantage of that fact to only calculate force vectors once and then apply them to both.

Double check your constant values

It seems to me that 3e4 is a somewhat large value for EPS. Additionally, in the calculation of the random mass, the math is currently 1.98892e30rand()10 + 1e20, but why not simply write 1.98892e31*rand() + 1e20 instead?

Eliminate return 0 at the end of main

When a C++ program reaches the end of main the compiler will automatically generate code to return 0, so there is no reason to put return 0; explicitly at the end of main`.

StackExchange Code Review Q#87341, answer score: 15