# Simple n-body class in C++

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;
}

• @Samuel If you want to provide more information, I think it'd be great to have an explanation of what is the N body simulation. This way, people could understand what your code is supposed to achieve. – IEatBagels Oct 23 '19 at 13:56
• Do you have a specific version of C++ you are targeting (C++11, C++14, C++17)? – pacmaninbw Oct 23 '19 at 14:31
• @pacmaninbw I'm not sure which version to take. What makes sense? Which version should I use when writing a C++ programm? – recipe_for_disaster Oct 23 '19 at 19:43
• You can add the comments, but at this point don't change the code. The version you should probably learn and use at this point is C++17. – pacmaninbw Oct 23 '19 at 21:17
• @Samuel after considering our review feedback and possibly incorporating it into your code, you are very welcome to post it in a follow-up question. I am interested in how this piece of code develops further. – Roland Illig Oct 24 '19 at 5:54

## The Great!

You're not making the basic beginner error of using using namespace std;! The main() function is only 3 lines of code.

The function declarations in the nbody class that don't change things include const which will help optimization later.

The code utilizes the C++ random number generation rather than the C srand() and rand() functions.

Because Nbody was implemented as a class it is very easy to change main() so that it can accept user input for the values of n, dt and t_max.

The #include <string> is missing from nbody.cpp; this is necessary when compiling the code in most cases.

## The Obsolete

The use of inline function declarations is now only a suggestion to the compiler. Optimizing compilers can and will do a better job of optimizing by inlining code based.

The body of the Nbody constructor uses an obsolete form of initialization, rather than using () as in the following code

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


use braces {}:

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


Putting the initialization on a separate line makes it easier to find.

## Prefer STL Container Classes

Prefer STL container classes such as std::vector or std::array over old C style arrays. The std::array<type, size> class is a fixed size array. The std::vector<type> is a variable sized array. STL container classes provide iterators so that pointers aren't necessary. The use of std::vector<Particle> p; might reduce the number of parameters to the constructor. It would definitely remove the need for the variable n within the Nbody class since p.size() would always contain the number of particles after Nbody::init_data() has run. Also after Nbody::init_data() has run iterators could be used to access the particles in p and would allow the code to use a ranged for loop such as

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 (auto particle : p) {
results << t << " ";
results << particle.m << " ";
for (int d=0; d<DIM; d++) {
results << particle.x[d] << " ";
}
for (int d=0; d<DIM; d++) {
results << particle.v[d] << " ";
}
for (int d=0; d<DIM; d++) {
results << particle.F[d] << " ";
}
results << std::endl;
}
results.close();
}
}


Another benefit of making p an STL container class is that the destructor for the class Nbody can then be a default constructor and the array of particles doesn't need to be allocated in the class declaration.

## Variable Names

It's not really clear by just reading the code what the variables n_, n, dt_, dt, t_max_, t_max, x, F and v and p are. For instance, I assume dt means Delta Time, but it is not clear that this is true. The array p might be renamed particles, if I'm correct about dt than deltaTime might be more appropriate.

Yes there are comments for some of the variable names, but if I had to maintain the code I'd rather work with code that was self-documenting than depending on comments.

Example

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 (auto particle : particles) {
results << t << " ";
results << particle.mass << " ";
for (int d=0; d<DIM; d++) {
results << particle.position[d] << " ";
}
for (int d=0; d<DIM; d++) {
results << particle.velocity[d] << " ";
}
for (int d=0; d<DIM; d++) {
results << particle.Force[d] << " ";
}
results << std::endl;
}
results.close();
}
}


## Style

Some, not all, developers prefer to see the public declarations of a class before the private declarations of a class. This is because it becomes easier to find the public interface of the class.

The function void init_data() is not necessary unless you are planning to have multiple constructors, it might be better to move that code into the constructor.

If the functions print_parameter() and print_data() are debug functions it might be better to put them within #ifdef DEBUG and #endif.

In the current implementation return 0; from main() is not necessary. If error handling code is added and there is a return 1; it might be better to keep it. It might also be better to use return EXIT_SUCCESS; and EXIT_FAILURE which are defined in cstdlib (#include <cstdlib>).

## Suggestions

It might be better to allow the user to name the output file that the results go into, either by input through a user interface or as part of the command line arguments. The name could default to the current file name in case the user doesn't specify one.

It might also be better to have only one output file.

• The variable names are very self descriptive for this kind of domain in my opinion. I think it is also safe to assume that a maintainer of a n-body interaction program has a physics background. – mcocdawc Oct 23 '19 at 21:41
• You never know whos gonna maintain your code next. Might be a person with a totally different background. The cost for the naming the variables properly is neglible, but the potential benefit large. – Hakaishin Oct 24 '19 at 9:00
• @Hakaishin the cost is that anyone with the domain knowledge to maintain the code has to translate from the common notation of the domain, $dt$, each time they read the code. IMO this makes it less maintainable. – Pete Kirkham Oct 24 '19 at 11:15
• @Hakaishin The physics equations use mathematical symbols (variables) and not their names and sentences for a reason. To make everything more readable. By using long names the code that implements mathematical equations (e.g. numerical methods) becomes harder to parse. – Vladimir F Oct 24 '19 at 17:44
• I've seen a lot of code with terse, mathematical-seeming variable names that were totally made-up (or came from one obscure paper), and those codebases are terrible. But d for change, and t for time, m for mass, F for force, etc., are really such basic, well-known notation. – hegel5000 Oct 24 '19 at 19:30

# First Of All

You are doing a great job as a beginner. I have been programming for 10 years and my code for a long time was much, much less readable that what you have written. That said:

# What needs Fixing

I'm not privy of all of the details of the n-body problem but I have an idea of what it does. I'm not an expert on numerical accuracy so I won't comment on the arithmetic you are performing. Here are a few things that I see from a design perspective.

## This class is effectively Impossible to test

Between randomizing the input data upon construction and having one method that does the vast majority of the work, it is very difficult to write meaningful automated tests for this class. This is in part because this class does way too much.

## The Public interface does not reflect its usage

The public interface is much broader than what a client would use. As far as I can tell, the only thing a client would need to do is construct one of these objects and immediately call timeIntegration() upon it, then record the results somehow. More on this later.

## You use Non standard ways to convey standard concepts

You provide a "print_data" and a "write_data" method. The dependency on <iostream> & <fstream> is needless for this class and will make it very difficult to test in an automated (read: unit test) fashion. You should provide a << operator for the particle class instead and allow the client to decide what to do with the results.

## There's no way to get at the raw data for this class

Furthermore, since the print_data() and write_data() methods are seemingly the only way get data from this class, the use of this class in anything other than a simple command prompt program is limited. A method to get the internal data in non-printed form would be helpful.

## What to do

A better design for this class may be a public constructor with the necessary parameters which immediately calls everything necessary to compute the integration, and then a method to get the data that has been processed. Nothing else would be public. This way, it is very difficult for a client to use this class incorrectly. A class with a getter for its only owned data should raise a red flag in an OOP design, so all of this rethinking is really leading to a bigger realization that...

# This shouldn't be a class

My biggest consideration would be to not have this be a class at all. None of the data that it owns are invariant across the useful public interface. More on invariants in class design here on Wikipedia. There is no reason for the state that has been introduced to be owned by this class across its lifetime and there are plenty of opportunities to use this class in ways that produce completely invalid data. This instead should have an interface that consists of one high, level function.

The public interface to the n-body calculator should take in two or three things:

1. A settings struct. This will include all necessary pieces to properly run the calculation other than the "hot" data. this will be initialized by the client. If the struct data is not valid (i.e. something that will be a denominator of zero), the function should exit with a return code of some sort (or exception if that's allowed in your environment and that's your thing). This should probably be taken by const l-value reference
2. A std::vector<Particle> by (possibly const l-value) reference, this is the input data to the n-body calculator
3. a time step to run for. This could be part of the settings struct, but in my mind it's distinctly different than the other concepts that would be in the settings struct.

This function should guarantee to either modify the std::vector<Particle> in place or return a transformed std::vector<Particle>. My personal preference is the latter, however depending on which version of C++ you are using, that can be inhibitive to good performance. In essence, all that this function is doing is transforming a list of particle states. It can (and should) use other helper functions to do its work, and these functions would very likely get reused in other parts of a larger particle framework. All functions should be stateless other than the particle set passed in.

The value add from this multi-fold:

1. It is more obvious how to use this interface correctly. See the principle of least surprise. Wiki article.
2. It is much, much easier to test a set of stateless functions than it is to test a big, entangled class.
3. This will allow much higher reuse of basic operations as this code base expands.

# Other Suggestions

## Names

I would suggest better names for the Particle struct members. If they are used correctly in a larger program they will likely become ubiquitous as the base data types. There is nothing wrong with typing out mass, position, velocity and force. While it's true that people will probably know what you mean when you talk about position as x, they will definitely know what you mean when you type position.

## Strong Types

I would use strong types for the particle members. Jonathan Bocarra has some excellent blog articles on it on cppfluent (e.g. CppFluent Strong types). They can be treated the same as doubles, with the advantage of making it much more difficult to switch arguments around in function calls, and making the code more expressive.

## Get Rid Of The Globals

Globals are a bad thing, and should be avoided. Regardless of whether the object-oriented approach is gotten rid of, these should be incorporated into a settings struct of some kind.

## Use the STL more than you are

A lot of your summing for loops can use std::accumulate(); you should be using std::vectors rather than raw c-style arrays. You should be using range-based for loops where you can't use std::vector or an STL algorithm.

• You totally misunderstand the single responsibility principle. – gnasher729 Oct 24 '19 at 10:28
• You are right! I did. I will edit my answer – Marc Olberding Oct 24 '19 at 13:13
• Thank you! Great review in my opinion. – recipe_for_disaster Oct 25 '19 at 20:25

• Use unsigned integer type for DIM, Nbody.step, and Nbody.n since none of this can be negative;
• Use constexprsince C++11 instead just const for both DIM and EPS;
• Get rid of the unused argc and argv arguments in main;
• Consider more usage of const. For example f in Nbody::force() can be const, and a in Nbody::update_position can be const and so on.
• Comments are not for extended discussion; this conversation has been moved to chat. – Mathieu Guindon Oct 25 '19 at 16:27

Your code is written in a hybrid C/C++ style. For instance your destructor has a delete (I can't find where the corresponding new is) and that is basically never needed. Use a std::vector to store array-like data.

Also you do a lot of parameter passing like void Nbody::update_position(Particle *p). Use references instead, and use const Particle &p if the particle is only read.

Otherwise it looks like an n-body code to me. It's quadratic rather than something more sophisticated/efficient, but that's probably ok.

Oh, I've found the new: you have Particle *p = new Particle[n]; in the class definition, but n is uninitialized. That is probably undefined behavior, definitely extremely dangerous, and most likely completely wrong.

Don't use new to allocate an array! Use std::vector, as follows:

std::vector<Particle> the_particles;
public:
Particles(int n) : the_particles(vector<Particle>(n)) {}
}

• Indeed. As part of writing my answer, I refactored the code to use a std::vector<Particle>, and it felt much nicer. – Roland Illig Oct 23 '19 at 20:13
• I use new in the class definition. Did I use it the wrong way? What exactly is wrong with parameter passing? Don't I change the particle when I update it's velocity or position? Wouldn't it be wrong the use a constant reference then? Can you be more specific how comp_velocity() and update_velocity() would look like then? – recipe_for_disaster Oct 23 '19 at 21:47
• Can you please explain your code in more detail? I am not very familiar with the STL. Do I define std::vector<Particle> the_particles as private? What exactly is ParticleS? Where is it defined? How do I use it in the rest of my code? – recipe_for_disaster Oct 25 '19 at 7:28
• My compiler throws an error: error: expected unqualified-id before ‘int’ Particle(int n) : p(vector<Particle>(n)) {};. Do you see what I am doing wrong? – recipe_for_disaster Oct 25 '19 at 7:33
• 1. This uses C++17 idioms. Make sure to use the appropriate compiler flag! 2. Yes, make the vector private. 3. I used my own naming. It takes too much scrolling to find everything in your code. I'm sure you know what I mean. – Victor Eijkhout Oct 25 '19 at 15:48

The init_data function does not belong in the Nbody class. Nowhere in the definition of the N-body problem will you find the word "random", and using random input data is only connected to your particular situation, therefore this code should be moved into main.cpp.

In the constructor of Nbody, there's no need for the trailing underscore in the parameter names. The following code looks cleaner and is otherwise equivalent to your current code:

Nbody::Nbody(int n, double dt, double t_max)
: n(n), dt(dt), t_max(t_max) {
init_data();  // should be removed, as I said above
}


For debugging purposes it would be good to have not only the timeIntegration method, but also a simple step method that only does a single step. This allows you to write better unit tests. It also makes another of the constructor parameters, namely t_max unnecessary.

Still in timeIntegration, instead of step+=1 you should write ++step. Writing step++ would be equivalent, but that would tell every reader that you don't know C++ well. In C++ the ++ usually comes before the variable, in other languages like Java or C or Go it usually comes after the variable. See this Stack Overflow answer for some more details.

Comparing the code of timeIntegration with update_velocity reveals that you use an inconsistent programming style. You should decide for yourself whether to use camelCase or snake_case identifiers. Then, use that style consistently. Another thing is that you placed spaces around the operators * and /, but not around +. I would have expected it the other way round, since * and / bind the operands more tightly than +. The usual style is to always surround binary operators with spaces. Therefore t < t_max; t += dt; step++.

Your Nbody class does not account for tricky situations where the particles are so close together that dt becomes too large for a realistic simulation. This is something that you must document.

I like it that you separated updated_velocity and update_position into two separate methods. This makes them easy to read. (Plus, it's necessary from an implementation's point of view since you must first update the velocity of all particles before you can update any particle's position, otherwise the result depends on the ordering of the particles.)

The abbreviation comp in comp_position is ambiguous. It could mean compare or compute. You should spell it out.

In Nbody::force you should not name the parameters i and j, since these variable names are reserved for integers, by convention. I'd rather choose p and q. And if you rename Nbody::p to ps since it is plural anyway, there's no naming collision anymore.

In write_data the parameter step is not necessary since Nbody::step is accessible by the same name. You can just remove the parameter.

The method print_parameter should be called print_parameters since it is about all parameters, not just a single one.

At the API level, I would not put dt and t_max in the constructor but rather pass dt as parameter to the step method, and t_max as parameter to the timeIntegration method.

In nbody.h there is the EPS constant, which looks dubious. For a dt of 0.001 it may have an appropriate value of 0.0001, but what if I want to simulate using dt = 1.0e-9? I don't think it should be a global constant. Not even the speed of light should be, because there are so many different speeds of light, depending on the exact experiment.

In Nbody::init_data you wrote 1. without a trailing 0. Sure, it may save a single key stroke, but in my opinion it's not worth it. Just write the canonical 1.0, as you already did in several other places in the same function.

The data you write to the data_*.log files is quite imprecise. The typical double type provides 16 to 17 digits of precision, yet you only write out 6 of them, which is the C++ default. Since 2017, C++ finally supports printing floating point numbers accurately.

• C and Java both have pre-increment and post-increment just like C++ (I don't know Go, but would be surprised if that differed). It's true to say that we generally prefer pre-increment when the value is to be discarded (though you ignore that where writing the for loop's increment expression - why?), and that it may be more efficient for non-inlinable types, but post-increment does have its place in C++, and should be used when it's what's required. Perhaps that paragraph could be re-worded? – Toby Speight Oct 24 '19 at 7:02
• @TobySpeight Thanks, did that. – Roland Illig Oct 24 '19 at 17:37
• Thank you for this nice review! Is it correct to remove the trailing underscore also in the n-body class? Or is that a bad style? What do you mean by adding dt to the step method? Where would you put EPS then? – recipe_for_disaster Oct 25 '19 at 19:51
• Yes, it's ok to remove the underscores everywhere. Regarding the step method: I would define Nbody::step(double dt, double eps). This gives full control to the caller. – Roland Illig Oct 26 '19 at 6:34
• You said that init_data() does not belong in the Nbody class. Can you elaborate on this point a little? How do I initialize the particles then? Does it also make sense to add init_data() to the Particle struct? – recipe_for_disaster Oct 28 '19 at 16:16

# Use a vector math library

Find a suitable library that implements coordinate vectors, so you don't have to implement them as arrays of doubles. Ideally your struct Particle should look like:

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


And a suitable library will provide functions and operator overloads to make working with these types very easy. You should be able to write something like:

void Nbody::update_position(Particle *p) {
double a = dt * 0.5 / p->m;
p->x += dt * (p->v + a * p->F);
p->F_old = p->F;
}


There are many libraries available. I am partial to GLM myself. For a discussion of possible libraries, see https://stackoverflow.com/questions/1380371/what-are-the-most-widely-used-c-vector-matrix-math-linear-algebra-libraries-a.

# Make function manipulating Particles member functions of Particle

You have a lot of functions that mainly manipulate a particle's state, but they are not part of struct Particle itself. For example, update_position() is something that apart from the timestep dt only manipulates a Particle's member variables. If you make it a member function of Particle, it becomes a much cleaner looking function:

struct Particle {
...
void update_position(double dt);
};

void Particle::update_position(double dt) {
double a = dt * 0.5 / m;
x += dt * (v + a * F);
F_old = F;
}


And you call it like so:

void Nbody::comp_position() {
for (auto &p: particles) {
p.update_position(dt);
}
}


You can do the same for update_velocity(), and even force().

• I tried to incorporate first the second part of your answer but get the following errors: nbody.cpp: In member function ‘void Particle::update_position(double)’: nbody.cpp:47:33: error: ‘p’ was not declared in this scope const double a = dt * 0.5 / p->m; ^ nbody.cpp: In member function ‘void Nbody::comp_position()’: nbody.cpp:55:19: error: ‘particles’ was not declared in this scope for (auto &p: particles) {  – recipe_for_disaster Oct 25 '19 at 7:51
• Is it recommended to use for (auto &p: particles) or is for (Particle &p: particles) better? When should I not use auto? – recipe_for_disaster Oct 25 '19 at 15:36
• I recommend you use auto whenever the type is already clear from the context, or when you don't care about the type. Here, it should be obvious that p is a Particle, because it's an element of particles. Using auto in these cases avoids repeating types unnecessarily, and potentially avoids mistakes. – G. Sliepen Oct 25 '19 at 16:14
• Can you please outline how you would compute the force with your approach? I have difficulties since I have in the case of the force computation a nested loop and I am not sure how to do that in this case. – recipe_for_disaster Oct 25 '19 at 16:24
• Once you've written void force(Particle &other) as a member function of struct Particle, change force(&p[i], &p[j]) to p[i].force(p[j]) inside the nested loop. Here you should probably keep using int i and int j as iterator variables, as that's the easiest to ensure each pair is only visited once. – G. Sliepen Oct 25 '19 at 18:28

I can't comment due to being new here, but Roland Illig's assertion that it should be ++step and not step++ and that it shows that you do not understand C++ is incorrect.

In C++, the position of the ++ determines the order of how the expression is evaluated. So in ++step, the variable is incremented before any action with it is performed, while in step++, the action is performed before the value is incremented. Just having a step++ or ++step as a single line of code is basically equivalent, but the difference is apparent in an example like so:

int step = 0;
std::cout << ++step << std::endl; // would print 1
std::cout << step << std::endl; // would print 1


while

int step = 0;
std::cout << step++ << std::endl; // would print 0
std::cout << step << std::endl; // would print 1


Just clarifying this, as you should understand the difference as opposed to preferring one over the other for stylistic/reputation reasons!

• I think that Rolan Illig meant that ++i usually is faster than i++. – eanmos Oct 24 '19 at 6:52
• @eanmos, that is possible, but in the link you provided it states that this performance hit is usually only present if the type being incremented is not a primitive type. However in OP's post, step is of type int, which should be optimized correctly by the compiler. – ashiswin Oct 24 '19 at 7:03
• yes, but it is a good habit to write ++i by default and use post increment only when it is need to be used. – eanmos Oct 24 '19 at 7:11
• @eanmos fair enough. It was more of the assertion that "In C++ the ++ comes before the variable, in other languages like Java or C or Go it comes after the variable." But yes, the convention would be to use ++i. – ashiswin Oct 24 '19 at 7:12
• @Samuel, in short. i++ must create a copy of the old value, "return" it and then increment i. In the same time ++i doesn't have to create a copy and "return" it, it increments i and just "return" it. For build-in types (such as int) the compiler probably be able to optimize your code. But i can be an instance of a C++ class so i++ and ++i making calls to one of the operator++ function. In that case the compiler probably not be able to optimize and then copying of an object i when using i++ can be really slow. So in general case just use ++i. – eanmos Oct 24 '19 at 10:42

I will focus on one thing already addressed by another answer but that I think deserve more attention: the Single Responsibility Principle.

Your NBody class has several functionalities merged into one, which would be advisable to separate. It, as far as I can see:

• it represent a group of N particles
• it provides the algorithm to perform the physics simulation
• it provides the facilities to print the results of the simulation

I think there is enough material to separate these into three separate entities, leaving more flexibility for changing in the future.

Also, some of the methods in your NBody class actually act only on the given Particle, so they could be refactored as methods of the Particle struct.

Another suggestion is to take a look at the Template Method Pattern, which could be a useful starting point for the simulation framework to provide the right flexibility to change the integration method if it ever becomes necessary.

• Do you mean by separating that I write different classes? A particle class, a simulation class and a class that prints the results? – recipe_for_disaster Oct 24 '19 at 17:08

In addition to G. Sliepen’s idea, you could use the STL’s std::valarray<double>. This would let you replace something like

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];
}


with something like

p->F_old = p->F;
p->x += dt * (p->v + a * p->F);


It would also be possible to lay out a structure of arrays rather than an array of structures. If there are more particles than dimensions, this could let you perform wider vector operations on all the x-coordinates, then all the y-coordinates and all the z-coordinates, rather than being limited to the width of the coordinate system. That is, each p might have only two or three parallel computations, but if you have a number of std::array<std::valarray<double>, DIM> with the x-coordinates in x[0], the y-coordinates in x[1] and the z-coordinates in x[2], the velocities in v[0], etc., that might look like:

for (size_t i = 0; i < x.size(); ++i) {
F_old[i] = F[i];
x[i] += dt * (v[i] + a * F[i]);
}
`

and be able to use the full width of your vector registers. This would not, however, work as well if the computations are not so cleanly separable.