---

pip install numpy scipy matplotlib
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.integrate import solve_ivp

# Гравитационная постоянная и массы тел (для симуляции берем G=1 и равные массы)
G = 1.0
m = np.array([1.0, 1.0, 1.0])

def n_body_equations(t, y):
    """
    Вычисление производных для системы.
    Вектор состояния y содержит: [x1, y1, x2, y2, x3, y3, vx1, vy1, vx2, vy2, vx3, vy3]
    """
    r = y[:6].reshape((3, 2))
    v = y[6:]
    
    dvdt = np.zeros_like(r)
    # Вычисляем ускорения для каждого тела
    for i in range(3):
        for j in range(3):
            if i != j:
                diff = r[j] - r[i]
                dist = np.linalg.norm(diff)
                # Избегаем деления на ноль при столкновениях
                if dist > 0.05: 
                    dvdt[i] += G * m[j] * diff / dist**3
                    
    return np.concatenate((v, dvdt.flatten()))

# --- Начальные условия ---
# Позиции: Тело 1 (слева), Тело 2 (справа), Тело 3 (по центру)
r0 = np.array([[-1.0, 0.0], [1.0, 0.0], [0.0, 0.0]])

# Скорости для оригинальной системы (подобраны так, чтобы центр масс покоился)
v0_original = np.array([[0.5, 0.5], [-0.5, 0.5], [0.0, -1.0]])

# Скорости для изменённой системы: меняем vx первого тела на 0.001
v0_altered = v0_original.copy()
v0_altered[0, 0] += 0.001

# Собираем начальные векторы состояний (по 12 элементов каждый)
y0_original = np.concatenate((r0.flatten(), v0_original.flatten()))
y0_altered = np.concatenate((r0.flatten(), v0_altered.flatten()))

# --- Интегрирование ---
# Время симуляции: от 0 до 20 секунд, 2000 кадров для плавности
t_span = (0, 20)
t_eval = np.linspace(t_span[0], t_span[1], 2000)

print("Решение уравнений для первой системы...")
sol_orig = solve_ivp(n_body_equations, t_span, y0_original, t_eval=t_eval, rtol=1e-9, atol=1e-11)

print("Решение уравнений для второй (изменённой) системы...")
sol_alt = solve_ivp(n_body_equations, t_span, y0_altered, t_eval=t_eval, rtol=1e-9, atol=1e-11)

# Переформатируем результаты обратно в координаты [тело, координата, время]
r_orig = sol_orig.y[:6].reshape((3, 2, -1))
r_alt = sol_alt.y[:6].reshape((3, 2, -1))

# --- Анимация ---
fig, ax = plt.subplots(figsize=(9, 9))
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect('equal')
ax.set_title("Проблема трёх тел: Расхождение траекторий\nСплошные: Оригинал | Пунктир: $v_x + 0.001$")
ax.grid(True, linestyle='--', alpha=0.5)

colors = ['#FF3333', '#33FF33', '#3333FF'] # Красный, Зеленый, Синий

# Линии и точки для оригинальной системы
lines_orig = [ax.plot([], [], '-', color=colors[i], linewidth=2.5)[0] for i in range(3)]
pts_orig = [ax.plot([], [], 'o', color=colors[i], markersize=8, markeredgecolor='black')[0] for i in range(3)]

# Линии и точки для изменённой системы
lines_alt = [ax.plot([], [], '--', color=colors[i], linewidth=1.5, alpha=0.7)[0] for i in range(3)]
pts_alt = [ax.plot([], [], 'x', color=colors[i], markersize=8)[0] for i in range(3)]

# Длина "хвоста" траектории, чтобы не загромождать экран
tail_length = 250 

def init():
    for element in lines_orig + pts_orig + lines_alt + pts_alt:
        element.set_data([], [])
    return lines_orig + pts_orig + lines_alt + pts_alt

def update(frame):
    start = max(0, frame - tail_length)
    
    for i in range(3):
        # Оригинальная система
        lines_orig[i].set_data(r_orig[i, 0, start:frame], r_orig[i, 1, start:frame])
        pts_orig[i].set_data([r_orig[i, 0, frame]], [r_orig[i, 1, frame]])
        
        # Изменённая система
        lines_alt[i].set_data(r_alt[i, 0, start:frame], r_alt[i, 1, start:frame])
        pts_alt[i].set_data([r_alt[i, 0, frame]], [r_alt[i, 1, frame]])
        
    return lines_orig + pts_orig + lines_alt + pts_alt

print("Отрисовка анимации...")
ani = FuncAnimation(fig, update, frames=len(t_eval), init_func=init, blit=True, interval=10)

plt.tight_layout()
plt.show()