[Pet] Two-dimensional simulation of the interaction of celestial bodies in C ++

Introduction and introduction

I confess, until now I was a web developer and did not hold anything heavier than node in my hands. The world of pointers, references and (horror) typed arrays, and even fixed length, looked all the more terrible and mysterious for me. But tonight I decided to finally explore this world of deep dark fantasies. I’ve been dreaming for two years about my own cute 2D simulation of the movement of celestial bodies, and I’m about to write it on crosses!

By no means a guide for beginners, just an interesting article on the topic of their projects.

As a library for rendering graphics, I chose sfml, simply because it fell out in the search first. I still don’t quite understand how to upload a c++ project somewhere along with its dependencies, so you will have to install the lib and include it yourself if you want to test it.

Mathematics and Logic

The laws of my world are very simple – there is a celestial body, it has Cartesian coordinates x and ythe velocities of the corresponding components vx and vyas well as physical properties such as mass mradius r and some “density” factor d. I consider the radius as a product m * d for simplicity, and “density” in quotation marks, because it is not density in its direct sense, just a coefficient.

Each celestial body affects each in proportion to the distance – it can be calculated simply by the Pythagorean theorem, I think there is no need for explanations:
dist = sqrt( pow(x1 - x2, 2) + pow(y1 - y2, 2) );

And finally, the force of gravitational interaction between bodies is expressed as a formula:

${ F = G \cdot \frac{m_{1} \cdot m_{2}}{r^{2}} \cdot \frac{j_{12}}{r}}$

where j12 is the radius vector, calculated as x2-x1. This is required for vector calculations, such as the coordinate system, all that, I don’t know why, they didn’t explain this to us at school, bribes are smooth from me, we just use the formula (I would be grateful if people with an older education explain in the comments). We will use a randomly selected value as the gravitational constant G, since our bodies are disproportionate to the real world.

Speeds vx and vy change due to acceleration, which in turn is a quotient of the force obtained above and body mass 1. Thus, one mass will be reduced, and the formula for changing the speed will take the form:
vx += G * m2 / dist / dist * (x2 - x1) / dist vy += G * m2 / dist / dist * (y2 - y1) / dist

That’s about it, let’s move on to the easiest 😀

The code

The constants we need to work with are:

const double PI = 3.1415926536;
const double G = 1; //наша гравитационная постоянная
const int boundX = 1200; //размер окна по ширине
const int boundY = 800; //...и высоте

Further, for the celestial body, I created the following structure (why are structures needed if this is the same as classes?):

struct Body
{
	float x;
	float y;
	float vx;
	float vy;
	float m;
    float d;
	float r = m * d;
};

I was worried about performance, and I didn’t know how it would all work, so just in case, I got a timer on milliseconds that stops the main thread to give the code a breather:

void rest(int ms)
{
	std::this_thread::sleep_for(std::chrono::milliseconds(ms));
}

This is what main looks like in its entirety:

int main() {
	sf::ContextSettings settings;
	settings.antialiasingLevel = 8.0; //уровень сглаживания
	RenderWindow window(VideoMode(boundX, boundY), "Planet Simulator", sf::Style::Close, settings);

	int delay = 50; //50ms отдыха для кода

	while (window.isOpen())
	{
		std::vector<Body> bodiesToAdd = handleEvents(&window); //про обработку событий чуть позже

		window.clear(Color(0, 0, 0, 0)); //очищаем экран

		update(&window); //обновление экрана
		bodies.insert(bodies.end(), bodiesToAdd.begin(), bodiesToAdd.end()); //про это тоже позже :3

		window.display(); //отрисовка экрана

		rest(delay); //отдых
	}

	return 0;
}

The heart of the whole project is the method void update(RenderWindow* window). According to the laws described above, it updates the data of the bodies in the vector bodiesand then gives them to the screen for rendering:

void update(RenderWindow* window)
{
	std::set<int> deleteBodies; //номера небесных тела, которые столкнулись друг с другом и должны умереть
	int size = bodies.size(); //количество небесных тел

    //используем два цикла для фиксации влияния каждого тела на каждое
	for (int i = 0; i < size; i++)
	{
		Body& p0 = bodies[i]; //ссылка на текущее тело
		for (int j = 0; j < size; j++)
		{
            //помним, что под одинаковыми индексами лежит одно и тоже тело, а сами на
            //себя тела не влияют, поэтому пропускаем
			if (i == j)
				continue;

			Body& p = bodies[j]; //ссылка на второе тело
          
			double dist = sqrt(pow(p0.x - p.x, 2) + pow(p0.y - p.y, 2));

            //проверка коллизии тел
			if (dist > p0.r + p.r)
			{
                //собственно, изменение скоростей покомпонентно
				p0.vx += G * p.m / dist / dist * (p.x - p0.x) / dist;
				p0.vy += G * p.m / dist / dist * (p.y - p0.y) / dist;
			}
			else {
				deleteBodies.insert(i);
				deleteBodies.insert(j);
			}
		}

        //изменение координат покомпонентно
		p0.x += p0.vx;
		p0.y += p0.vy;
      
		CircleShape circle = renderBody(p0); //функция, которая на основе структуры просто создает кружок
		window->draw(circle); //рисуем кружок
	}

    //мой способ очистки вектора тел от тех, индексы которых помечены как уничтоженные
    //и содержатся в наборе deleteBodies
	std::vector<Body> copy_bodies;
	for (int i = 0; i < bodies.size(); ++i)
	{
		if (deleteBodies.find(i) == deleteBodies.end()) //если индекса в наборе нет
		{
			copy_bodies.push_back(bodies[i]);
		}
	}
	bodies = copy_bodies;
    //я уверен, что можно сделать это лучшее и быстрее, пока не знаю как.
}

About the collision check – everything is quite prosaic, just look at the picture and understand that if the distance between the radii (dist) is less than the sum of the radii, then such circles will exactly intersect, which means they will be killed.

The basic functionality is ready, but I would like to add the creation of planets on a mouse click! To do this, we will turn to events in sfml – this is where the function comes in handystd::vector<Body> handleEvents(RenderWindow* window), which returns a vector of bodies that I would like to add, which I then insert into the bodies vector in main. Here’s how it’s set up:

std::vector<Body> handleEvents(RenderWindow* window) {
	std::vector<Body> newBodies; //вектор с новыми телами
	Event event;

	while (window->pollEvent(event)) //пока в пуле есть новые события
	{
		if (event.type == Event::Closed) //обработка выхода из программы
			window->close();
		if (event.type == Event::MouseButtonPressed) //если случилось нажатие мыши
		{
			sf::Vector2i position = sf::Mouse::getPosition(*window);
			if (event.mouseButton.button == sf::Mouse::Left)
			{
                //по нажатию на левую кнопку добавляем просто планету
				newBodies.push_back(Body{ static_cast<float>(position.x), static_cast<float>(position.y), 0.0, 0.0, 100.0, 0.1});
			}
			else 
			{
                //по нажатию на правую кнопку добавляем целую тяжелую звезду
				newBodies.push_back(Body{ static_cast<float>(position.x), static_cast<float>(position.y), 0.0, 0.0, 1000.0, 0.05});
			}
		}
	}

	return newBodies;
}

Conclusion

All good things come to an end, and my post is no exception. Most of all, I would like to implement the force and direction of the initial movement of the planet through clicking and dragging the mouse (like in angry birds). I have no idea if this is possible in sfml, I hope for talent in the comments. I would be glad to hear suggestions for optimization, especially the moment of cleaning up dead planets. All the best, happiness to you ♡

Full code

… looks a little different than in the examples. In particular, it added jokes such as the fact that the stars eat the planets that crashed into them, and gain a little mass, and so on. Basically the same, enjoy watching

the code

#include <SFML/Graphics.hpp>

#include <chrono>
#include <thread>
#include <vector>
#include <set>

using namespace sf;

const double PI = 3.1415926536;
const double coef = 1;
const int boundX = 1200;
const int boundY = 800;
const float sunMass = 1000;
const float radiusCoef = 0.1;

struct Body
{
	float x;
	float y;
	float vx;
	float vy;
	float m;
	float r = m * radiusCoef;
};

std::vector<Body> bodies = { Body{500.0, 300.0, 0.0, 0.0, 1000.0} };

void sleep(int ms)
{
	std::this_thread::sleep_for(std::chrono::milliseconds(ms));
}

std::vector<Body> handleEvents(RenderWindow* window) {
	std::vector<Body> newBodies;
	Event event;

	while (window->pollEvent(event))
	{
		if (event.type == Event::Closed)
			window->close();
		if (event.type == Event::MouseButtonPressed)
		{
			sf::Vector2i position = sf::Mouse::getPosition(*window);
			if (event.mouseButton.button == sf::Mouse::Left)
			{
				newBodies.push_back(Body{ static_cast<float>(position.x), static_cast<float>(position.y), 0.0, 0.0, 100.0});
			}
			else 
			{
				newBodies.push_back(Body{ static_cast<float>(position.x), static_cast<float>(position.y), 0.0, 0.0, 2000.0});
			}
		}
	}

	return newBodies;
}

CircleShape renderBody(Body& b)
{
	b.r = b.m * radiusCoef;
	CircleShape circle(b.r);
	circle.setOrigin(b.r, b.r);
	circle.setPosition(b.x, b.y);
	if (b.m >= sunMass)
	{
		circle.setFillColor(Color(246, 222, 1));
	}
	return circle;
}

void update(RenderWindow* window)
{
	std::set<int> deleteBodies;
	int size = bodies.size();

	for (int i = 0; i < size; i++)
	{
		Body& p0 = bodies[i];
		for (int j = 0; j < size; j++)
		{
			if (i == j)
				continue;

			Body& p = bodies[j];
			double d = sqrt(pow(p0.x - p.x, 2) + pow(p0.y - p.y, 2));

			if (d > p0.r + p.r)
			{
				p0.vx += coef * p.m / d / d * (p.x - p0.x) / d;
				p0.vy += coef * p.m / d / d * (p.y - p0.y) / d;
			}
			else {
				if (p0.m >= sunMass && p.m >= sunMass)
				{
					deleteBodies.insert(i);
					deleteBodies.insert(j);
				}
				else
				{
					if (p0.m < sunMass)
					{
						deleteBodies.insert(i);
					}
					else {
						p0.m += p.m * 0.1;
					}

					if (p.m < sunMass)
					{
						deleteBodies.insert(j);
					}
					else {
						p.m += p0.m * 0.1;
					}
				}
			}
		}
		p0.x += p0.vx;
		p0.y += p0.vy;

		CircleShape circle = renderBody(p0);
		window->draw(circle);
	}

	std::vector<Body> copy_bodies;

	for (int i = 0; i < bodies.size(); ++i)
	{
		if (deleteBodies.find(i) == deleteBodies.end())
		{
			copy_bodies.push_back(bodies[i]);
		}
	}

	bodies = copy_bodies;
}

int main() {
	sf::ContextSettings settings;
	settings.antialiasingLevel = 8.0;
	RenderWindow window(VideoMode(boundX, boundY), "Planet Simulator", sf::Style::Close, settings);

	int msForFrame = 50;

	while (window.isOpen())
	{
		std::vector<Body> bodiesToAdd = handleEvents(&window);

		window.clear(Color(0, 0, 0, 0));

		update(&window);
		bodies.insert(bodies.end(), bodiesToAdd.begin(), bodiesToAdd.end());
		bodiesToAdd.clear();

		window.display();

		sleep(msForFrame);
	}

	return 0;
}

108