# N body simulation in C++

You can find the complete code on github. Since I'm rather new to C++ I thought that it might be a good idea to write a simulation of the n-body problem to apply some of the concepts that I learned in class. The relevant part can be found in main.cpp on the above linked github repo, but the imporant part is given by:

#include <string>
#include <eigen3/Eigen/Dense>
#include <utility>
#include <iostream>
#include <chrono>
#include <cmath>
#include <fstream>
#include "integrate.hpp"
#include <vector>

using mass_t = double const;
using name_t = std::string const;
using vector_t = Eigen::Vector3d;
using phase_t = Eigen::VectorXd;
using rhs_t = phase_t (*)(double const &, phase_t const &);
using step_t = phase_t (*)(rhs_t , phase_t const &, double, double);

double square(double x){  // the functor we want to apply
return std::pow(x,2);
}

double norm(vector_t const& v){
return std::sqrt(v.unaryExpr(&square).sum());
}

// Class representing a stellar object
class StellarObject{
public:
explicit StellarObject(  vector_t pos = vector_t(0,0,0),       //Position of the sun relative to itself
vector_t speed = vector_t(0,0,0),     //Velocity of the sun relative to itself
mass_t const& mass = 1.00000597682,         //Mass of the Sun [kg]
name_t const& name = "Sun") : mass_(mass),
name_(name),
position_(std::move(pos)),
velocity_(std::move(speed)) {};
vector_t get_position() const {
return position_;
}
vector_t get_velocity() const {
return velocity_;
}
double get_mass() const {
return mass_;
}
private:
vector_t position_;
vector_t velocity_;
name_t name_;
mass_t mass_;
};

phase_t nbody_prod(double const & t, phase_t const & z,
std::vector<double> const& masses, double const & G){
int const N = masses.size();
phase_t q(3*N);                     //space for positions of the particles
q << z.head(3*N);                   //fill it with the relevant elements from z
phase_t ddx(3*N);                   //space for the rhs of the ODE
ddx.fill(0);                       //make sure that set to zero everywhere
for(std::size_t k = 0; k<N; ++k){
vector_t tmp;                       //tmp object to store ddx
tmp.fill(0);
vector_t qk = q.segment(3*k,3); //get position of k-th object
for(std::size_t i=0; i < N; ++i){
if(i!=k){                       //calculate acceleration of k-th object
vector_t qi = q.segment(3*i,3);
vector_t diff = qi-qk;
double normq = std::pow(norm(diff),3);
tmp += masses[i]/normq * diff;
}
}
ddx.segment(3*k,3) = G * tmp;
}
return ddx;
}

// z = [r_(1,x), r_(1,y), r_(1,z), r_(2,x), .... r_(N,z), v_(1,x), v_(1,y), ..., v_(N, z)]
phase_t nbody_rhs(double const & t, phase_t const & z, std::vector<double> const& masses,
double const G){
int const N = masses.size();
phase_t rhs(6*N);                     //space for the rhs
rhs.head(3*N) = z.tail(3*N);        // fill velocities in first 3*N elements of rhs
rhs.tail(3*N) = nbody_prod(t, z, masses, G); // fill last 3*N elements with nbody_prod
return rhs;
}

static double G;
static std::vector<double> masses;

Eigen::MatrixXd n_body_solver(std::vector<StellarObject> const & planets){
int const N = planets.size();
phase_t z0(6*N);
for(std::size_t k=0; k < N; ++k) {
z0.segment(3*k,3) = planets[k].get_position();
}
for(std::size_t k=N; k < 2*N; ++k) {
z0.segment(3*k,3) = planets[k%N].get_velocity();
}
for(auto & planet : planets){
masses.push_back(planet.get_mass());
}
G = 2.95912208286e-4;
auto reduced_rhs = [](double const & t, phase_t const & z0) {return nbody_rhs(t, z0, masses, G);};
return explicit_midpoint(reduced_rhs, z0, 60000, 40000);
}


Note that I left out the main-function of main.cpp since there is nothing interesting happening there.. Just the setup for a simple test.

The code works, it produces results that are sensible (if you have a working phyton environment and are using a Unix based machine you can run run.sh file inside the repo and should get a nice animation of the stellar movement), but there are several things that bother me about the code.

First of all: I'd like to put all of these different function into one solver class, which I think would make it much easier to use and understand. The goal would be something like the following code (but I can't get it to work...)

class n_body_solver{
public:
n_body_solver(std::vector<StellarObject> const & planets,
double const & G=2.95912208286e-4) : N_(planets.size()) {
G_= G;
for(std::size_t k=0; k < N_; ++k) {
z0_.segment(3*k,3) = planets[k].get_position();
}
for(std::size_t k=N_; k < 2*N_; ++k) {
z0_.segment(3*k,3) = planets[k%N_].get_velocity();
}
for(auto & planet : planets){
masses_.push_back(planet.get_mass());
}
}

phase_t nbody_prod(double const & t, phase_t const & z,
std::vector<double> const& masses, double const & G){
int const N = masses.size();
phase_t q(3*N);                     //space for positions of the particles
q << z.head(3*N);                   //fill it with the relevant elements from z
phase_t ddx(3*N);                   //space for the rhs of the ODE
ddx.fill(0);                       //make sure that set to zero everywhere
for(std::size_t k = 0; k<N; ++k){
vector_t tmp;                       //tmp object to store ddx
tmp.fill(0);
vector_t qk = q.segment(3*k,3); //get position of k-th object
for(std::size_t i=0; i < N; ++i){
if(i!=k){                       //calculate acceleration of k-th object
vector_t qi = q.segment(3*i,3);
vector_t diff = qi-qk;
double normq = std::pow(norm(diff),3);
tmp += masses[i]/normq * diff;
}
}
ddx.segment(3*k,3) = G * tmp;
}
return ddx;
}

// z = [r_(1,x), r_(1,y), r_(1,z), r_(2,x), .... r_(N,z), v_(1,x), v_(1,y), ..., v_(N, z)]
phase_t nbody_rhs(double const & t, phase_t const & z, std::vector<double> const& masses,
double const G){
int const N = masses.size();
phase_t rhs(6*N);                     //space for the rhs
rhs.head(3*N) = z.tail(3*N);        // fill velocities in first 3*N elements of rhs
rhs.tail(3*N) = nbody_prod(t, z, masses, G); // fill last 3*N elements with nbody_prod
return rhs;
}

void solve(double const & time_intervall, unsigned int const & bins){
auto reduced_rhs = [](double const & t, phase_t const & z0) {return nbody_rhs(t, z0, masses_, G_);};
result_ = explicit_midpoint(reduced_rhs, z0_, time_intervall, bins);
}

Eigen::MatrixXd get_result() const {
return result_;
}

private:
phase_t z0_;
unsigned short const N_;
static double G_;
static std::vector<double> masses_;
Eigen::MatrixXd result_;
};


Other than that I'd be happy to receive feedback regarding general C++ coding style, choices of containers, types, etc. Anything that you notice when reading the code.

• You should not assume that each potential reviewer knows what the n-body problem is. A simple link to Wikipedia helps, you can just edit your question. Other than that, it's a good question. – Roland Illig Feb 7 '20 at 21:23
• @RolandIllig That's absolutely fair, edited! – Sito Feb 7 '20 at 21:34
• I'm not very familiar with Eigen, but is it really intended to be used like this? You are taking a Nx3 matrix and iterating over it in a n^2 fashion. I assume Eigen has some way to abstract that so that you're using fewer loops, and possibly get a benefit of hardware support. Is this code actually leveraging Eigen in a way that's better than just using a std::vector of struct {float x,y,z;}; ? – butt Feb 7 '20 at 22:46
• @butt To be honest, I don't know. I'm not that familiar with Eigen either, actually, I did this exercise here precisely to become a bit familiar with it (almost all of the functions I used here from Eigen, I used for the first time...). So if there is a faster (better) way of doing this, I'm more than open to hear about it. – Sito Feb 7 '20 at 23:16
• Here is another n body c++ program you can reference codereview.stackexchange.com/questions/231191/…. A reason to leave an uninteresting main in the program is so that others can test the code and see how the classes are used. – pacmaninbw Feb 8 '20 at 14:40

## 1 Answer

Your code makes me think that you don't understand what exactly you want to model.

explicit StellarObject(
vector_t pos = vector_t(0,0,0),       //Position of the sun relative to itself
vector_t speed = vector_t(0,0,0),     //Velocity of the sun relative to itself
mass_t const& mass = 1.00000597682,   //Mass of the Sun [kg]
name_t const& name = "Sun"
)


When creating a (model of a) stellar object, there is absolutely nothing that should be related to the sun. The sun is just a random star and has no relation to any of the fundamental physical constants.

The comments to the right reference the sun though, which makes the whole code look wrong. I very much doubt that the sun's mass is a single kilogram.[citation needed]

It is wrong to provide default values for any of the parameters of this constructor since it doesn't make sense to have 5 objects with the same name or the same position or the same mass, that's just not a realistic scenario. The caller of this constructor must be forced to think about all these values explicitly.

G = 2.95912208286e-4;


This is a magic number. The symbol G usually stands for the gravitational constant, whose value is approximately $$\ 6.67430\cdot10^{-11} \frac{\text{m}^3}{\text{kg} \cdot \text{s}^2}\$$. Your value of $$\2.9\cdot10^{-4}\$$ is nowhere near that value, therefore you must document where you got that number from and what its dimension is. In physical simulations, it's important to carry the dimensions around in the calculations to prevent typos and other mistakes. For example it doesn't make sense to add seconds to meters and divide by Ampère.

Always use the Internation System of Units, don't accept any other measurement units unless you document exactly why you have to use different units and how you did the conversion.

vector_t tmp;    //tmp object to store ddx


The variable name tmp is terrible, it should be forbidden. You should have named it next_ddx, since it collects the positions after the next simulation step. This naming scheme would also suggest a better name for vector_t qk and q. To avoid confusion, these should be called next_pk and next_p. Sure, the names are a bit longer, but the name q does not really tell (me) much, except that it is the letter following p in several Latin alphabets. If that's a well-known naming convention among physicists, it's ok if the code is ever only read by physicists.

double normq = std::pow(norm(diff),3);
tmp += masses[i]/normq * diff;


It's confusing to see a division by $$\r^3\$$ when I only expect a division by $$\r^2\$$. The code would be easier to understand if you just divided by $$\r^2\$$ first and did the direction calculations afterwards and independently. The code is easy to understand if the commonly known formulas like $$\F = \text{G} \cdot \frac{m_1 \cdot m_2}{r^2}\$$ appear exactly in this form. Every deviation from this makes the code more difficult to read and to verify.

I did not analyze the rest of the code in detail. I saw many helpful comments that explained the short variable names, which is good for understanding the code. I also saw many questionable comments that contradicted the code, and these are bad. You should probably read the code aloud to someone else and while doing that, listen to your words to see whether they make sense. If they don't, the code is wrong. Or the model of the world you are building. Either way, something needs to be fixed in these cases.