Modelling solar system, Runge-kutta 4 (n-Body problem)

Question

I have a program that simulations all planets, where the forces due to each planet is considered during each time step.I'm looking to see where I could make some improvements.

main.cpp

// PlanetarySimulator.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <iostream>
#include <fstream>
#include <map>
#include "Planet.h"
#include "MotionNBody.h"
const double days = 400; // 365 days in seconds
const double TimeStep = 24* 60 * 60; // 1 day in seconds




int main()
{
    std::map<std::string, Planet> planets;

    // Initialise Planets

    Planet sun(1.989E30, 695510E3, "Sun");
    sun.x.push_back(0);
    sun.vy.push_back(0);
    sun.y.push_back(0);
    sun.vx.push_back(0);

    Planet mercury(0.3301E24, 2439700, "Mercury");
    mercury.x.push_back(-46000000000);
    mercury.vy.push_back(-58980);
    mercury.y.push_back(0);
    mercury.vx.push_back(0);
    planets.emplace(mercury.getName(), mercury);

    Planet venus(4.8675E24, 6051800, "Venus");
    venus.x.push_back(-107480000000);
    venus.vy.push_back(-35260);
    venus.y.push_back(0);
    venus.vx.push_back(0);
    planets.emplace(venus.getName(), venus);

    Planet earth(5.972E24, 6371E3, "Earth");
    earth.x.push_back(-147095000000);
    earth.vy.push_back(-30300.0);
    earth.y.push_back(0);
    earth.vx.push_back(0);
    planets.emplace(earth.getName(), earth);

    Planet mars(6.4171E23, 3389500, "Mars");
    mars.x.push_back(-206620000000);
    mars.vy.push_back(-26500);
    mars.y.push_back(0);
    mars.vx.push_back(0);
    planets.emplace(mars.getName(), mars);

    Planet jupiter(1898.19E24, 71492000, "Jupiter");
    jupiter.x.push_back(-740520000000);
    jupiter.vy.push_back(-13720);
    jupiter.y.push_back(0);
    jupiter.vx.push_back(0);
    planets.emplace(jupiter.getName(), jupiter);

    Planet saturn(568.34E24, 54364000, "Saturn");
    saturn.x.push_back(-1352550000000);
    saturn.vy.push_back(-10180);
    saturn.y.push_back(0);
    saturn.vx.push_back(0);
    planets.emplace(saturn.getName(), saturn);

    Planet uranus(86.813E24, 24973000, "Uranus");
    uranus.x.push_back(-2741300000000);
    uranus.vy.push_back(-7110);
    uranus.y.push_back(0);
    uranus.vx.push_back(0);
    planets.emplace(uranus.getName(), uranus);

    Planet neptune(102.413E24, 24341000, "Neptune");
    neptune.x.push_back(-4444450000000);
    neptune.vy.push_back(-5500);
    neptune.y.push_back(0);
    neptune.vx.push_back(0);
    planets.emplace(neptune.getName(), neptune);






    MotionNBody motion(sun, planets);
    motion.setSimulationDays(days);
    motion.setSimulationTimeStep(TimeStep);
    motion.simulate();
    for (auto & planet : planets)
    {
        std::ofstream myfile;
        myfile.open("C:\\Users\\matkinso\\Documents\\Training_Projects\\PlanetarySimulatior\\PlanetarySimulatior\\Files\\" + planet.first + ".csv");
        myfile << planet.first << std::endl << "x,y" << std::endl;
        for (int i = 0; i < planet.second.x.size(); i++)
        {
            myfile << planet.second.x[i] << "," << planet.second.y[i] << std::endl;
        }
        myfile.close();
    }
}

Planet.h

#pragma once
#include <string>
#include <vector>

class Planet
{
public:
    Planet(double mass, double radius, std::string  name);
    ~Planet();

    double getMass();
    double getRadius();
    std::string getName();

    std::vector<double> x;
    std::vector<double> y;
    std::vector<double> vx;
    std::vector<double> vy;

private:
    double mass;            // In kg
    double radius;          // In m
    std::string name;   
};

planet.cpp

#include "Planet.h"

Planet::Planet(double mass, double radius, std::string name) :
    mass(mass),
    radius(radius),
    name(name)
{
}


Planet::~Planet()
{
}

double Planet::getMass()
{
    return mass;
}

double Planet::getRadius()
{
    return radius;
}

std::string Planet::getName()
{
    return name;
}

MotionNBody.h

#pragma once
#include <map>
#include <vector>
#include "Planet.h"

class MotionNBody
{
    // take a map of planet name to planet object, and for each time step
    // calculate the new position and velocity for each planet at the same time
    // taking into account all the other planets.
public:
    MotionNBody(Planet & sun, std::map<std::string, Planet> & planets);
    ~MotionNBody();
    void setSimulationDays(double days);
    void setSimulationTimeStep(double time);
    void simulate();

private:
    void doRK4( const std::string & currentPlanet, int i);
    void calcalateNewPositions(const std::string & planetName, int i);
    double calculateAccelerationX(double x, double y, int i, const std::string & currentPlanet);
    double calculateAccelerationY(double x, double y, int i, const std::string & currentPlanet);

    std::map<std::string,Planet> & planets;
    Planet & sun;
    double timeStep;
    double simDays;
};

MotionNBody.cpp

#include "MotionNBody.h"
#include <math.h>


MotionNBody::MotionNBody(Planet & sun, std::map<std::string, Planet> & planets) :
    sun(sun),
    planets(planets)
{
}


MotionNBody::~MotionNBody()
{
}

void MotionNBody::setSimulationDays(double days)
{
    simDays = days;
}

void MotionNBody::setSimulationTimeStep(double time)
{
    timeStep = time;
}


void MotionNBody::simulate()
{
    // For each time step
    // Loop through all planets.
    // summing the acceleration due to all other planets (i!=j) 
    // and calculating the new position.

    for (int i = 0; i < simDays; i++)
    {
        for (auto & planet : planets)
        {
            calcalateNewPositions(planet.first, i);
        }
    }
}

void MotionNBody::calcalateNewPositions(const std::string & planetName, int i)
{
    doRK4(planetName, i);
}

void MotionNBody::doRK4(const std::string & currentPlanet, int i)
{
    double cX  = planets.find(currentPlanet)->second.x[i];
    double cY  = planets.find(currentPlanet)->second.y[i];
    double cVX = planets.find(currentPlanet)->second.vx[i];
    double cVY = planets.find(currentPlanet)->second.vy[i];

    double k1x  = cVX;
    double k1y  = cVY;
    double k1vx = calculateAccelerationX(cX, cY, i, currentPlanet);
    double k1vy = calculateAccelerationY(cX, cY, i, currentPlanet);

    double k2x  = cVX + timeStep / 2 * k1vx;
    double k2y  = cVY + timeStep / 2 * k1vy;
    double k2vx = calculateAccelerationX(cX + timeStep / 2 *k1x, cY + timeStep / 2 * k1y, i, currentPlanet);
    double k2vy = calculateAccelerationY(cX + timeStep / 2 * k1x, cY + timeStep / 2 * k1y, i, currentPlanet);

    double k3x  = cVX + timeStep / 2 * k2vx;
    double k3y  = cVY + timeStep / 2 * k2vy;
    double k3vx = calculateAccelerationX(cX + timeStep / 2 * k2x, cY + timeStep / 2 * k2y, i, currentPlanet);
    double k3vy = calculateAccelerationY(cX + timeStep / 2 * k2x, cY + timeStep / 2 * k2y, i, currentPlanet);

    double k4x = cVX + timeStep * k3vx;
    double k4y = cVY + timeStep * k3vy;
    double k4vx = calculateAccelerationX(cX + timeStep * k3x, cY + timeStep * k3y, i, currentPlanet);
    double k4vy = calculateAccelerationY(cX + timeStep * k3x, cY + timeStep * k3y, i, currentPlanet);

    double newX  = cX + timeStep / 6 * (k1x, + 2 * k2x + 2 * k3x + k4x);
    double newY  = cY + timeStep / 6 * (k1y, + 2 * k2y + 2 * k3y + k4y);
    double newVX = cVX + timeStep / 6 * (k1vx, +2 * k2vx + 2 * k3vx + k4vx);
    double newVY = cVY + timeStep / 6 * (k1vy, +2 * k2vy + 2 * k3vy + k4vy);

    planets.find(currentPlanet)->second.x.push_back(newX);
    planets.find(currentPlanet)->second.y.push_back(newY);
    planets.find(currentPlanet)->second.vx.push_back(newVX);
    planets.find(currentPlanet)->second.vy.push_back(newVY);
}

double MotionNBody::calculateAccelerationX(double x, double y, int i, const std::string & currentPlanet)
{
    double ax = 0;
    for (auto & planet : planets)
    {
        if (planet.first != currentPlanet)
        {
            auto xo = planet.second.x[i];
            auto yo = planet.second.y[i];
            auto m = planet.second.getMass();

            double firstTerm = (6.67408E-11 * m) / ((xo - x)*(xo - x) + (yo - y)*(yo - y));
            double secondTerm = ((xo - x)) / (std::sqrt((xo - x)*(xo - x) + (yo - y)*(yo - y)));
            ax += firstTerm * secondTerm;
        }
    }
    auto xo = sun.x[0];
    auto yo = sun.y[0];
    auto m = sun.getMass();

    double firstTerm = (6.67408E-11 * m) / ((xo - x)*(xo - x) + (yo - y)*(yo - y));
    double secondTerm = ((xo - x)) / (std::sqrt((xo - x)*(xo - x) + (yo - y)*(yo - y)));
    ax += firstTerm * secondTerm;


    return ax;
}

double MotionNBody::calculateAccelerationY(double x, double y, int i, const std::string & currentPlanet)
{
    double ay = 0;
    for (auto & planet : planets)
    {
        if (planet.first != currentPlanet)
        {
            auto xo = planet.second.x[i];
            auto yo = planet.second.y[i];
            auto m = planet.second.getMass();

            double firstTerm = (6.67408E-11 * m) / ((xo - x)*(xo - x) + (yo - y)*(yo - y));
            double secondTerm = ((yo - y)) / (std::sqrt((xo - x)*(xo - x) + (yo - y)*(yo - y)));
            ay += firstTerm * secondTerm;
        }
    }
    auto xo = sun.x[0];
    auto yo = sun.y[0];
    auto m = sun.getMass();

    double firstTerm = (6.67408E-11 * m) / ((xo - x)*(xo - x) + (yo - y)*(yo - y));
    double secondTerm = ((yo - y)) / (std::sqrt((xo - x)*(xo - x) + (yo - y)*(yo - y)));
    ay += firstTerm * secondTerm;

    return ay;
}
```

Answer 1

Welcome to the Code Review website, nice first question. You did a pretty good job with the programming, I can't comment on whither Runge-kutta 4 (n-Body problem) was properly implemented because I don't have the background in astronomy.

First I have one question about the symbolic constant days, the comment 365 days in seconds confuses me.

Magic Numbers

The code would be more readable/understandable if there were more numeric constants, for instance :

const double SolarMass = 1.989E30;
const double SolarRadius = 695510E3;
const double MercuryMass = 0.3301E24;
const double MercuryRadius = 2439700;

including the mass and radius for all the bodies involved.

Numeric constants in code are sometimes referred to as Magic Numbers, because there is no obvious meaning for them. There is a discussion of this on stackoverflow.

DRY Code

There is a programming principle called the Don't Repeat Yourself Principle sometimes referred to as DRY code. If you find yourself repeating the same code multiple times it is better to encapsulate it in a function. If it is possible to loop through the code that can reduce repetition as well.

Complexity

The function main() is too complex (does too much). As programs grow in size the use of main() should be limited to calling functions that parse the command line, calling functions that set up for processing, calling functions that execute the desired function of the program, and calling functions to clean up after the main portion of the program.

There is also a programming principle called the Single Responsibility Principle that applies here. The Single Responsibility Principle states:

that every module, class, or function should have responsibility over a single part of the functionality provided by the software, and that responsibility should be entirely encapsulated by that module, class or function.

By adding 2 functions to the class MotionNBody, changing the reference to the map planets to an actual map in the class MotionNBody and modifying the constructors of the class Planet and the class MotionNBody main is greatly simplified.

// PlanetarySimulator.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <iostream>
#include "Planet.h"
#include "MotionNBody.h"
const double days = 400; // 365 days in seconds
const double TimeStep = 24* 60 * 60; // 1 day in seconds
const double SolarMass = 1.989E30;
const double SolarRadius = 695510E3;
const double MercuryMass = 0.3301E24;
const double MercuryRadius = 2439700;

int main()
{
    // Initialise Planets

    Planet sun(SolarMass , SolarRadius, "Sun", 0, 0, 0, 0);
    MotionNBody motion(sun);
    motion.AddPlanet(MercuryMass , MercuryRadius, "Mercury", -46000000000, -58980);
    motion.AddPlanet(4.8675E24, 6051800, "Venus", -107480000000, -35260);
    motion.AddPlanet(5.972E24, 6371E3, "Earth", -147095000000, -30300.0);
    motion.AddPlanet(6.4171E23, 3389500, "Mars", -206620000000,-26500);
    motion.AddPlanet(1898.19E24, 71492000, "Jupiter", -740520000000, -13720);
    motion.AddPlanet(568.34E24, 54364000, "Saturn", -1352550000000, -10180);
    motion.AddPlanet(86.813E24, 24973000, "Uranus", -2741300000000, -7110);
    motion.AddPlanet(102.413E24, 24341000, "Neptune", -4444450000000, -5500);

    motion.setSimulationDays(days);
    motion.setSimulationTimeStep(TimeStep);
    motion.simulate();
    motion.OutputResults(".\\");

}

Here are the changes to the class Planet :

class Planet {
public:
    Planet(double mass, double radius, std::string name, double x, double vy, double y, double vx);
    ~Planet() = default;

    double getMass() { return mass; };
    double getRadius() { return radius; };
    std::string getName() { return name; };

    std::vector<double> x;
    std::vector<double> y;
    std::vector<double> vx;
    std::vector<double> vy;

private:
    double mass;            // In kg
    double radius;          // In m
    std::string name;
};


Planet::Planet(double mass, double radius, std::string name, double xIn, double vyIn, double yIn, double vxIn) :
        mass(mass),
        radius(radius),
        name(name)
{
    x.push_back(xIn);
    vy.push_back(vyIn);
    y.push_back(yIn);
    vx.push_back(vxIn);
}

Here are the changes to the class MotionNBody

class MotionNBody {
    // take a map of planet name to planet object, and for each time step
    // calculate the new position and velocity for each planet at the same time
    // taking into account all the other planets.
public:
    MotionNBody(Planet & sun);
    ~MotionNBody() = default;
    void AddPlanet(double mass, double radius, std::string name, double x = 0, double vy = 0, double y = 0, double vx = 0);
    void OutputResults(std::string filePath);
    void setSimulationDays(double days) { simDays = days; };
    void setSimulationTimeStep(double time) { timeStep = time; };
    void simulate();

private:
    void doRK4( const std::string & currentPlanet, int i);
    void calcalateNewPositions(const std::string & planetName, int i);
    double calculateAccelerationX(double x, double y, int i, const std::string & currentPlanet);
    double calculateAccelerationY(double x, double y, int i, const std::string & currentPlanet);

    std::map<std::string,Planet> planets;
    Planet & sun;
    double timeStep;
    double simDays;

};

MotionNBody::MotionNBody(Planet & sun) :
        sun(sun)
{
}

void MotionNBody::AddPlanet(double mass, double radius, std::string name, double x, double vy, double y, double vx)
{
    planets.emplace(name, Planet(mass, radius, name, x, vy, y,vx));
}

void MotionNBody::OutputResults(std::string filePath)
{
    for (auto & planet : planets)
    {
        std::ofstream myfile;
        myfile.open(filePath + planet.first + ".csv");
        myfile << planet.first << std::endl << "x,y" << std::endl;
        for (int i = 0; i < planet.second.x.size(); i++)
        {
            myfile << planet.second.x[i] << "," << planet.second.y[i] << std::endl;
        }
        myfile.close();
    }

}

Simplifying Classes and Allowing the Compiler to do Some Work for You

For empty destructors in classes you can allow the compiler to do some work for you by declaring the destructors as default as shown above. For one line functions it is easier to put the code into the header file, there is no reason to define them in the .cpp file.

Answer 2

double newX  = cX + timeStep / 6 * (k1x, + 2 * k2x + 2 * k3x + k4x);
double newY  = cY + timeStep / 6 * (k1y, + 2 * k2y + 2 * k3y + k4y);
double newVX = cVX + timeStep / 6 * (k1vx, +2 * k2vx + 2 * k3vx + k4vx);
double newVY = cVY + timeStep / 6 * (k1vy, +2 * k2vy + 2 * k3vy + k4vy);

Those mid-line commas seem like a bug.

Your compiler should warn you about it - you may need to turn up your compiler warning level.

---

double cX  = planets.find(currentPlanet)->second.x[i];
double cY  = planets.find(currentPlanet)->second.y[i];
double cVX = planets.find(currentPlanet)->second.vx[i];
double cVY = planets.find(currentPlanet)->second.vy[i];

...

planets.find(currentPlanet)->second.x.push_back(newX);
planets.find(currentPlanet)->second.y.push_back(newY);
planets.find(currentPlanet)->second.vx.push_back(newVX);
planets.find(currentPlanet)->second.vy.push_back(newVY);

We're finding the same planet many times over. We could just find it once and keep a reference to it. More to the point, we already have the planet available in simulate (in planet.second):

for (int i = 0; i < simDays; i++)
{
    for (auto & planet : planets)
    {
        calculateNewPositions(planet.first, i);
    }
}

We could pass the planet by reference instead of the name.

---

For some reason we're using a map of planets by name... when each planet already stores the name (so we might as well put them in a vector). So we'd have:

std::vector<Planet> planets;

...

for (int i = 0; i != simDays; ++i)
    for (auto & planet : planets)
        calculateNewPositions(planet, i);

...

void MotionNBody::doRK4(Planet& planet, int i)
{
    double cX  = planet.x[i]; // ta da! no finding

...

---

        double firstTerm = (6.67408E-11 * m) / ((xo - x)*(xo - x) + (yo - y)*(yo - y));
        double secondTerm = ((xo - x)) / (std::sqrt((xo - x)*(xo - x) + (yo - y)*(yo - y)));
        ax += firstTerm * secondTerm;

Isn't it possible for us to be dividing by zero here?

There's also lot of duplicate calculations (xo - x and (xo - x)*(xo - x) + (yo - y)*(yo - y)) that could be put in named variables.

---

...

double k2x  = cVX + timeStep / 2 * k1vx;
double k2y  = cVY + timeStep / 2 * k1vy;
double k2vx = calculateAccelerationX(cX + timeStep / 2 *k1x, cY + timeStep / 2 * k1y, i, currentPlanet);
double k2vy = calculateAccelerationY(cX + timeStep / 2 * k1x, cY + timeStep / 2 * k1y, i, currentPlanet);

...

There's a lot of repeated code for the x and y coordinates here. We should use an existing vector math library or write our own struct Vector2{ double x, y; }; with appropriate operator overloads to get something more like:

Vector2 k2  = cV + timeStep / 2.0 * k1v;
Vector2 k2v = calculateAcceleration(c + timeStep / 2.0 * k1, i, planet);