---

simulation.py:
from body import *
import numpy as np
# from itertools import combinations

RES = WIDTH, HEIGHT = 1400, 750
FPS = 60

running = True
paused = False

pygame.init()
screen = pygame.display.set_mode(RES)
clock = pygame.time.Clock()

scale = 300
camera_x, camera_y = 0, 0
camera_up = 10

AU = 1.495978707e11
DAY = 86400
SUN_MASS = 1.9885e30
G = 6.67430e-11 * SUN_MASS * DAY ** 2 / (AU ** 3)
r_aph = 1.0167
a = 1
v_aph = np.sqrt(G * (2.0 / r_aph - 1.0 / a))

dt = 0.01
visual_speed = 10
speed_per_FPS = visual_speed / FPS
steps_per_frame = int(speed_per_FPS / dt)

# ПАРАМЕТРЫ ЗЕМЛИ
earth_a = 1.0  # большая полуось в AU
earth_e = 0.0167  # эксцентриситет
earth_r_aph = earth_a * (1 + earth_e)  # 1.0167 AU
earth_mass = 3.0027e-6  # масса Земли в массах Солнца (точное значение)
earth_v_aph = np.sqrt(G * (2.0 / earth_r_aph - 1.0 / earth_a))

# ПАРАМЕТРЫ ЛУНЫ
moon_a = 0.00257  # большая полуось орбиты Луны в AU (384,400 км)
moon_period_real = 27.321661  # реальный сидерический период в днях

# Рассчитаем необходимую скорость для этого периода
# Из 3-го закона Кеплера: T = 2π√(a³/μ), где μ = G * (M_earth + M_moon)
moon_mass = earth_mass / 81.3
mu_moon = G * (earth_mass + moon_mass)

# Орбитальная скорость для круговой орбиты с периодом moon_period_real
moon_v_circ = 2 * np.pi * moon_a / moon_period_real

# Альтернативно: из формулы круговой орбиты v = √(μ/a)
moon_v_circ2 = np.sqrt(mu_moon / moon_a)

# ИНИЦИАЛИЗАЦИЯ ТЕЛ
Sun = Body(1, [0, 0, 0], [0, 0, 0], 0.004649,  10, (255, 255, 150))

# Меркурий
Mercury = Body(earth_mass*0.055274,
             [solar(semi_major_axis=0.38709927, eccentricity=0.20563593, g=G)[0], 0, 0],
             [0, solar(semi_major_axis=0.38709927, eccentricity=0.20563593, g=G)[1], 0],
              1.63e-5, 3, (80, 60, 50))

# Венера
Venus = Body(earth_mass*0.815,
             [solar(semi_major_axis=0.723332, eccentricity=0.0068, g=G)[0], 0, 0],
             [0, solar(semi_major_axis=0.723332, eccentricity=0.0068, g=G)[1], 0],
              4.05e-5, 5, (200, 220, 0))

# Земля
Earth = Body(earth_mass,
             [solar(semi_major_axis=1, eccentricity=0.01671123, g=G)[0], 0, 0],
             [0.0, solar(semi_major_axis=1, eccentricity=0.01671123, g=G)[1], 0],
              4.26e-5, 5, (0, 0, 255))

# Луна
Moon = Body(moon_mass,
            [earth_r_aph + moon_a, 0.0, 0.0],  # апогей
            [0.0, earth_v_aph + moon_v_circ, 0.0],  # скорость Земли + орбитальная
             1.16e-5, 2, (255, 255, 255))

# Марс
Mars = Body(earth_mass*0.107,
             [solar(semi_major_axis=1.523662, eccentricity=0.0933941, g=G)[0], 0, 0],
             [0.0, solar(semi_major_axis=1.523662, eccentricity=0.0933941, g=G)[1], 0],
             2.27e-5, 4, (200, 70, 40))

# Юпитер
Jupiter = Body(earth_mass*317.8,
             [solar(semi_major_axis=5.204267, eccentricity=0.048775, g=G)[0], 0, 0],
             [0.0, solar(semi_major_axis=5.204267, eccentricity=0.048775, g=G)[1], 0],
             0.000478, 7, (120, 80, 70))


bodies = []

Sun.orbital.append(Mercury)
Sun.orbital.append(Venus)
Sun.orbital.append(Earth)
Sun.orbital.append(Mars)
Sun.orbital.append(Jupiter)

Earth.orbital.append(Moon)

bodies.append(Sun)
bodies.append(Mercury)
bodies.append(Venus)
bodies.append(Earth)
bodies.append(Moon)
bodies.append(Mars)
bodies.append(Jupiter)


for body in bodies:
    body.acc = np.zeros(3)
    for other in bodies:
        if other != body:
            r_vec = other.pos - body.pos
            r = np.linalg.norm(r_vec)
            if r > 0:
                body.acc += G * other.mass * r_vec / (r ** 3)
    body.acc_prev = body.acc.copy()


a = 0
while running:
    a += 1
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

        if event.type == pygame.MOUSEBUTTONDOWN and event.button == 1:
            mouse_x, mouse_y = pygame.mouse.get_pos()
            world_x = (mouse_x - WIDTH/2 - camera_x) / scale
            world_y = (HEIGHT - mouse_y - HEIGHT/2 - camera_y) / scale

            new = Body(1 / 333030, [world_x, world_y, 0], [0, 0, 0], 4.26e-5, 7, (255, 0, 255))
            bodies.append(new)


        elif event.type == pygame.KEYDOWN:
            if event.key == pygame.K_SPACE:
                paused = not paused

            if event.key == pygame.K_s:
                camera_y += camera_up*scale/300

            if event.key == pygame.K_w:
                camera_y -= camera_up*scale/300

            if event.key == pygame.K_d:
                camera_x -= camera_up*scale/300

            if event.key == pygame.K_a:
                camera_x += camera_up*scale/300

            if event.key == pygame.K_q:
                center_world_x = (-camera_x) / scale
                center_world_y = (-camera_y) / scale

                scale /= 1.5

                camera_x = -center_world_x * scale
                camera_y = -center_world_y * scale

            if event.key == pygame.K_e:
                center_world_x = (-camera_x) / scale
                center_world_y = (-camera_y) / scale

                scale *= 1.5

                camera_x = -center_world_x * scale
                camera_y = -center_world_y * scale

            if event.key == pygame.K_z:
                camera_up /= 2

            if event.key == pygame.K_x:
                camera_up *= 2


    if not paused:
        for i in range(steps_per_frame):
            # Способ 1: Использование отдельных методов leapfrog
            # Шаг 1: Обновление скорости на половину шага
            for body in bodies:
                body.update_velocity_half_step(dt)

            # Шаг 2: Обновление позиций
            for body in bodies:
                body.update_position(dt)

            # Шаг 3: Расчет новых ускорений
            for body in bodies:
                body.update_acceleration(bodies, G)

            # Шаг 4: Обновление скорости на вторую половину шага
            for body in bodies:
                body.update_velocity_full_step(dt)


    screen.fill('black')

    for body in bodies:
        for point in body.trajectory[::1]:
            pygame.draw.circle(screen, body.color,
                                   (int(point[0] * scale + WIDTH/2 + camera_x),
                                    HEIGHT - (point[1] * scale + HEIGHT/2 + camera_y)), 1)


        if scale < 100:
            pygame.draw.circle(screen, body.color,
                               (int(body.pos[0] * scale + WIDTH/2 + camera_x),
                                HEIGHT - int(body.pos[1] * scale + HEIGHT/2 + camera_y)), 3)

        elif scale > 500 and body.radius * scale > 1:
            pygame.draw.circle(screen, body.color,
                               (int(body.pos[0] * scale + WIDTH/2 + camera_x),
                                HEIGHT - int(body.pos[1] * scale + HEIGHT/2 + camera_y)), body.radius * scale)

        else:
            pygame.draw.circle(screen, body.color,
                               (int(body.pos[0] * scale + WIDTH/2 + camera_x),
                                HEIGHT - int(body.pos[1] * scale + HEIGHT/2 + camera_y)), body.display_radius)

    pygame.display.flip()
    clock.tick(FPS)

pygame.quit()

body.py:
import pygame
import numpy as np


class Body:
    def __init__(self, mass, pos, vel, radius, display_radius, color):
        self.mass = mass
        self.pos = np.array(pos, dtype=float)
        self.vel = np.array(vel, dtype=float)
        self.acc = np.zeros(3)  # Текущее ускорение
        self.acc_prev = np.zeros(3)  # Предыдущее ускорение (для leapfrog)
        self.radius = radius
        self.display_radius = display_radius
        self.color = color
        self.trajectory = []
        self.orbital = []

    def update_position(self, dt):
        """Обновление позиции методом leapfrog"""
        self.pos += self.vel * dt + 0.5 * self.acc_prev * dt ** 2
        self.trajectory.append([self.pos[0], self.pos[1], self.pos[2]])

    def update_velocity_half_step(self, dt):
        """Половина шага для скорости (первая половина leapfrog)"""
        self.vel += 0.5 * self.acc_prev * dt

    def update_acceleration(self, bodies, G):
        """Расчет нового ускорения"""
        self.acc_prev = self.acc.copy()  # Сохраняем предыдущее ускорение
        self.acc = np.zeros(3)

        for body in bodies:
            if body != self:
                r_vec = body.pos - self.pos
                r = np.linalg.norm(r_vec)
                if r > 0:
                    self.acc += G * body.mass * r_vec / (r ** 3)

    def update_velocity_full_step(self, dt):
        """Вторая половина шага для скорости"""
        self.vel += 0.5 * self.acc * dt

    # Альтернативный вариант: один метод для leapfrog целиком
    def leapfrog_update(self, bodies, dt, G):
        """Полный шаг leapfrog"""
        # Сохраняем текущее ускорение
        acc_old = self.acc.copy()

        # Обновляем позицию
        self.pos += self.vel * dt + 0.5 * acc_old * dt ** 2
        self.trajectory.append([self.pos[0], self.pos[1], self.pos[2]])

        # Рассчитываем новое ускорение
        self.acc = np.zeros(3)
        for body in bodies:
            if body != self:
                r_vec = body.pos - self.pos
                r = np.linalg.norm(r_vec)
                if r > 0:
                    self.acc += G * body.mass * r_vec / (r ** 3)

        # Обновляем скорость
        self.vel += 0.5 * (acc_old + self.acc) * dt


def solar(semi_major_axis, eccentricity, g):
    r_aph = semi_major_axis * (1 + eccentricity)
    v_aph = np.sqrt(g * (2.0 / r_aph - 1.0 / semi_major_axis))

    return r_aph, v_aph