Создаем симулятор солнечной системы

Предисловие
Вечная тяга к новому подтолкнула к изучению такого замечательного языка программирования, как Питон. Как это часто бывает, отсутствие идеи, на реализацию которой не жалко потратить свое время, сильно тормозило процесс.

Волею судьбы на глаза попался замечательный цикл статей о создании игры-платформера на Питоне
тут и тут.
Я решил взяться за один старый проект. За симулятор движения тел под действием сил гравитации.

Что из этого вышло читайте дальше.

Часть первая. Теоритическая

Чтобы решить задачу, нужно сначала четко себе ее представить.
Предположим, всеми правдами и неправдами нам удалось заполучить двумерный участок безвоздушного пространства с находящимися в нем телами. Все тела перемещаются под действием сил гравитации. Внешнего воздействия нет.
Нужно построить процесс их движения относительно друг друга. Простота реализации и красочность конечного результата послужат стимулом и наградой. Освоение Питона будет хорошей инвестицией в будущее.

Введем систему координат.

Пускай наша система состоит из двух тел:
1. массивной звезды массой М и центром (x0, y0)
2. легкой планеты массой m, с центром в точке (x, y), скоростью v = (vx, vy) и ускорением a = (ax, ay).

Когда нам удастся разобрать этот случай, мы легко перейдем к сложным системам со взаимным влиянием звезд и планет друг на друга. Сейчас же речь пойдет о самом простом.

После несложных манипуляций со вторым законом Ньютона, законом всемирного тяготения и подобными треугольниками, я нашел, что:

ax = G * M * (x0-x) / r^3
ay = G * M * (y0-y) / r^3

Это позволяет составить алгоритм перемещения планеты в поле гравитации звезды:

1. Перед началом задаем начальное положение планеты (x, y) и начальную скорость (vx, vy)
2. На каждом шаге вычисляем новое ускорение по формуле выше, после этого пересчитываем скорость и координаты:

vx := vx + T * ax
vy := vy + T * ax

x := x + T * vx
y := y + T * yx

Осталось разобраться с константами G и T. Положим G = 1. Для нашей задачи это не так важно. Параметр T влияет на точность и скорость вычислений. Тоже положим 1 для начала.

Часть вторая. Практическая

Итак, моя первая программа на Питоне. При этом еще раз хочется поблагодарить Velese за практическое руководство.

import pygame, math from pygame import * from math import *  WIN_WIDTH = 800 WIN_HEIGHT = 640 PLANET_WIDTH = 20 PLANET_HEIGHT = 20 DISPLAY = (WIN_WIDTH, WIN_HEIGHT) SPACE_COLOR = "#000022" SUN_COLOR = "yellow" PLANET_COLOR = "blue"  #Sun position X0 = WIN_WIDTH // 2 Y0 = WIN_HEIGHT // 2 #Sun mass M0 = 5000 #Stop conditions CRASH_DIST = 10 OUT_DIST = 1000  def main():     #PyGame init     pygame.init()     screen = pygame.display.set_mode(DISPLAY)     pygame.display.set_caption("Solar Mechanics v0.1")          #Space init     bg = Surface((WIN_WIDTH,WIN_HEIGHT))     bg.fill(Color(SPACE_COLOR))          draw.circle (bg, Color(SUN_COLOR), (X0, Y0), 10)                          #Timer init                          timer = pygame.time.Clock()          #Planet init     planet = Surface((PLANET_WIDTH, PLANET_HEIGHT))     planet.fill(Color(SPACE_COLOR))     draw.circle (planet,         Color(PLANET_COLOR),         (PLANET_WIDTH // 2, PLANET_HEIGHT // 2),         5)      #Planet to Sun distance     r = 0.0     #Initial planet pos, speed and accel     x = 100.0     y = 290.0     vx = 0.1     vy = 1.5     ax = 0.0     ay = 0.0      done = False     while not done:         timer.tick(50)         for e in pygame.event.get():             if e.type == QUIT:                 done = True                 break                  r = sqrt((x - X0)**2 + (y - Y0)**2)                  ax = M0 * (X0 - x) / r**3         ay = M0 * (Y0 - y) / r**3          #New spped based on accel         vx += ax         vy += ay          #New pos based on speed         x += vx         y += vy          screen.blit(bg, (0, 0))               screen.blit(planet, (int(x), int(y)))         pygame.display.update()               if r < CRASH_DIST:             done = True             print("Crashed")             break         if r > OUT_DIST:             done = True             print("Out of system")             break      #Farewell     print (":-)")  if __name__ == "__main__":     main() 

Так выглядит наша система после некоторого времени симуляции

Пока писалась эта заметка, симулятор разросся новой функциональностью: количество объектов в звездной системе не ограничевается, учитывается взаимное их влияние друг на друга, расчетная часть вынесена в свой класс, конфигурация системы задается в отдельном файле и добавлена возможность выбора систем.

Сейчас я занимаюсь поиском интересных сценариев системы и небольшими улучшениями интерфейса.

Вот пример того, что на данный момент в разработке:

Если эта заметка встретит положительные отзывы, обещаю продолжить рассказ о более новой версии.