Вступление и подводка
Каюсь, до сего момента я был веб-разработчиком и ничего тяжелее node в руках не держал. Тем страшнее и загадочнее для меня выглядел мир указателей, ссылок и (о ужас) типизированных массивов, да еще и фиксированной длины. Но сегодня вечером я решился наконец-то исследовать этот мир deep dark fantasies. Я джва года мечтал о своей собственной няшной двумерной симуляции движения небесных тел, и я собрался писать её на крестах!
Ни в коем случае не пособие для новичков, просто интересная статья на тему своих проектов.
В качестве библиотеки для отрисовки графики я выбрал sfml, просто потому что он выпал в поиске первым. Я еще не совсем хорошо понимаю, как заливать куда-либо c++ проект вместе с его зависимостями, поэтому вам придется самостоятельно устанавливать либу и подключать ее, если вы захотите это потестить.
Математика и логика
Законы моего мира совсем простые - имеется небесное тело, у него есть декартовы координаты x и y, скорости соответствующих компонент vx и vy, а также физические свойства, такие как масса m, радиус r и некий коэффициент "плотности" d. Радиус я считаю как произведение m * d для простоты, а "плотность" в кавычках, потому что это не плотность в прямом ее понимании, просто коэффициент.
Каждое небесное тело влияет на каждое соразмерно расстоянию - его можно посчитать просто по теореме Пифагора, думаю тут объяснения излишни:dist = sqrt( pow(x1 - x2, 2) + pow(y1 - y2, 2) );
И наконец, сила гравитационного взаимодействия между телами выражается в виде формулы:
где j12 - радиус-вектор, высчитываемый как x2-x1. Это требуется для векторных подсчетов, типа система координат все такое я сам хз почему, в школе нам такого точно не объясняли, с меня взятки гладки, просто пользуемся формулой (буду благодарен, если в комментариях люди с более старым образованием объяснят). В качестве гравитационной постоянной G будем использовать случайно подобранное значение, поскольку наши тела несоразмерны реальному миру.
Скорости vx и vy изменяются благодаря ускорению, которое в свою очередь является частным полученной выше силы и массы тела 1. Таким образом, одна масса сократится, и формула изменения скорости примет вид:vx += G * m2 / dist / dist * (x2 - x1) / dist
vy += G * m2 / dist / dist * (y2 - y1) / dist
На этом, пожалуй все, переходим к самому легкому :D
Код
Константы, которые потребуются нам для работы:
const double PI = 3.1415926536;
const double G = 1; //наша гравитационная постоянная
const int boundX = 1200; //размер окна по ширине
const int boundY = 800; //...и высотеДалее, для небесного тела я создал следующую структуру (зачем нужны структуры если это тоже самое что классы?):
struct Body
{
float x;
float y;
float vx;
float vy;
float m;
float d;
float r = m * d;
};Я переживал за производительность, да и не знал как будет работать это все, поэтому на всякий случай обзавелся таймером на миллисекундах, который стопает главный поток, чтобы дать коду отдышаться:
void rest(int ms)
{
std::this_thread::sleep_for(std::chrono::milliseconds(ms));
}Вот так выглядит main целиком:
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;
}Сердцем всего проекта служит метод void update(RenderWindow* window). Он по описанным выше законам обновляет данные тел в векторе bodies, а затем отдает их экрану для отрисовки:
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;
//я уверен, что можно сделать это лучшее и быстрее, пока не знаю как.
}Про проверку коллизии - все достаточно прозаично, достаточно взглянуть на картинку и понять, что если расстояние между радиусами (dist) меньше сумм радиусов, то такие окружности точно пересекутся, а значит убьются.
Базовый функционал готов, но я хотел бы добавить еще создание планет по клику мыши! Для этого мы обратимся к событиям в sfml - тут нам и пригодится функцияstd::vector<Body> handleEvents(RenderWindow* window), возвращающая вектор тел, которые я бы хотел добавить, и которые я затем в main вставляю в вектор bodies. Вот как она устроена:
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;
}Заключение
Все хорошее имеет свойство кончаться, и мой пост к сожалению не исключение. Больше всего мне хотелось бы реализовать силу и направление начального движения планеты через клик и оттягивание мыши (как в энгри бердс). Без понятия, возможно ли такое сделать в sfml, надеюсь на талантов в комментариях. Был бы рад услышать предложения оптимизации, особенно момента очистки убитых планет. Всем добра, счастья вам ♡
Полный код
... выглядит чуточку иначе, как в примерах. В частности, в нем добавлены приколы типа того, что звезды едят планеты, которые в них врезались, и набирают массу понемногу, и прочее. Базово все так же, приятного просмотра
код
#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;
}
Комментарии (23)
Naf2000
25.11.2022 22:22+4По сути такой же делал в школе на turbo pascal 6.0
IvanPetrof
26.11.2022 05:32+2Аналогично.. Только паскаль был пятый))
Удивительно, что статьи на эту тему стабильно "лайкаются". Такое уже было на Хабре на разных языках. На JavaScript, например. И судя по комментам, каждый раз это вызывает восхищение, как от магии.
Одно время я тоже думал написать в песочницу подобную статью (ещё когда это небыло мейнстримом). Но видимо из-за повышенного уровня синдрома самозванца, мне казалось что такая банальная вещь не достойна целой статьи.
a-tk
25.11.2022 22:27+1Накопление ошибок интегрирования? Не, не слышал.
yatanai
25.11.2022 23:40Как я знаю вроде нету нормальных алгоритмов для такого? Ты либо энергию теряешь, либо фаза сдвигается постоянно, не?
Поясню то что "знаю" - когда ты меняешь состояние, ты пишешь типо "+2м/с" к объекту. Потеря энергии означает что сумма всех изменений не будет равна тому, что мы имеем в процессе интегрирования. А "фаза", как я это понял, то что изменения непостоянны/"слегка смещены" по времени, от чего график с интегрированными vs запись в массив, будет постоянно сдвигаться на какой-то шаг, хотя энергия может сохраняться.
PS за матан вообще не шарю, самоучка
Soft_touch_plastic Автор
26.11.2022 20:28Спасибо за информацию, буду читать что за ошибки интегрирования такие)
Emulyator
25.11.2022 23:54Не уверен, что алгоритм проверки столкновения отработает корректно, если за одну итерацию достаточно большое приращение координат "пропустит" момент, когда "расстояние между радиусами (
dist) меньше сумм радиусов". Особенно на встречных, но смещенных от прямой курсах?
brotchen
25.11.2022 23:56Делал такую симуляцию в детстве. Проблема была в том, что орбиты были неустойчивы. Возможно, дело не только в ошибках интегрирования, но и в том, что такая система дифференциальных уравнений сама по себе не является асимптотически устойчивой.
s_f1
26.11.2022 00:36Больше всего мне хотелось бы реализовать силу и направление начального движения планеты через клик и оттягивание мыши (как в энгри бердс).
По нажатию мыши запоминаешь координаты x1, y1. По отпусканию получаешь координаты x2, y2. Вектор {x1 — x2, y1 — y2} будет и направлением, и скоростью движения, только домножь на нужный коэффициент, чтобы было адекватно остальной картине мира.a-tk
26.11.2022 07:59+1только домножь на нужный коэффициент, чтобы было адекватно остальной картине мира
так и получаются монстры с урановыми ломами в ртути...
AMDmi3
26.11.2022 00:39+1Я еще не совсем хорошо понимаю, как заливать куда-либо c++ проект вместе с его зависимостями, поэтому вам придется самостоятельно устанавливать либу и подключать ее, если вы захотите это потестить.
Вместе с зависимостями никто проекты не заливает. Достаточно выложить на github репозиторий с двумя файлами - исходником и CMakeLists.text/meson.build/Makefile в зависимости от того чем он собирается. По-хорошему ещё с лицензией и readme.
phenik
26.11.2022 06:30+1Может полезно будет для развития проекта — симуляция N-тел.
Soft_touch_plastic Автор
26.11.2022 20:39большое спасибо, я не ожидал если честно такого теплого приема и советов????
AxeFizik
26.11.2022 07:00+3Советую ознакомиться с курсом https://classroom.udacity.com/courses/cs222 там рассматриваются методы численного решения диффуров, и в первом же уроке рассматривается пример гравитационного взаимодействия тел. Это поможет вам сделать симуляцию более корректной
Un_ka
26.11.2022 20:31+1Для полноты картины хотя бы вставили гифку результата симуляции.
Также может быть интересна двумерная симуляция по другим законам тяготения:
Soft_touch_plastic Автор
26.11.2022 20:36По всей видимости буду доделывать проект, учитывая советы выше и ваш комментарий. Спасибо!
TheCalligrapher
26.11.2022 22:21В данном случае, конечно, не критично, учитывая специфику задачи, но тем не менее: использование
powдля возведения в целые степени всегда режет глаз. А уж в квадрат...Также: почему
float, а неdouble?Soft_touch_plastic Автор
26.11.2022 22:37По поводу pow - мне казалось наоборот каким-то оверхедом писать два раза одно и тоже, если можно воспользоваться готовой функцией.
По поводу float - не уловил разницы между ними, в силу своего нетипизированного пока еще мышления)
Как я понял, точности float должно хватить.
a-tk
27.11.2022 13:57+2Почитайте вот это на досуге, получите удовольствие: https://stihi.ru/diary/datomir/2013-06-01
TheGast
Прикол с таймером для "отдышаться" лишний. Это точно никому не нужно. Если что, в крайнем случае у тебя ядро процессора будет долбиться в сотку с максимальным фпс. Но лагать ничего не будет.
"Как в Angry Birds" делается тремя событиями (pressed уже обрабатывается в коде): https://www.sfml-dev.org/tutorials/2.5/window-events.php#the-mousebuttonpressed-and-mousebuttonreleased-events и https://www.sfml-dev.org/tutorials/2.5/window-events.php#the-mousemoved-event . На pressed надо запоминать координату нажатия, на moved - рисовать стрелку между координатой нажатия и текущей координатой мыши, а на released -- создавать объект и давать ускорение.
Soft_touch_plastic Автор
Понял, спасибо за совет)
playermet
Нужно, но с нулевой задержкой и только если выключен vsync. Как раз потому что грузить одно ядро на 100% приложению которое рисует два кружка это моветон. Люди на ноутбуках спасибо скажут.