import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from timeit import default_timer as timer

# Константы
G = 6.674080000*10**(-11)
M = (5*10**6*1.989)*10**30
e = 0.8
a = 20*1.496*10**11
c = 299792458
rg = 2*G*M/c**2 # Гравитационный радиус
Q = a*(1+e) # Апоцентр (с него начинаем)

# Засекаем время выполнения программы
start1 = timer()

# Решаем уравнения
def f(r,t,params):
    r1, r2, r3 = r
    l = params
    dvrdt = -rg*c**2/(2*r1**2) + l**2/r1**3 - 3*l**2*rg/r1**4
    dphidt = l/r1**2
    return [r2, dvrdt, dphidt]

# Список начальных значений
r0 = [Q,0,0]
# Количество шагов по времени
N = 100000

p = a*(1-e**2) # Фокальный параметр

L = np.sqrt(p*G*M) # Момент импульс 
params = L

T = np.sqrt(4*np.pi**2/(G*M))*a**(3/2) # Период обращения 
t = np.linspace(0,86*T, N)

# Решаем уравнение
sol = odeint(f,r0,t,args = (params,))

# Обозначение переменных из решения
radius = sol[:,0]
angle = sol[:,2]

X = radius*np.cos(angle)
Y = radius*np.sin(angle)

# Построение анимации
x_stars = np.random.uniform(-1.2*Q, 1.2*Q, 400)
y_stars = np.random.uniform(-1.2*Q, 1.2*Q, 400)
plt.style.use(['dark_background'])
fig = plt.figure(figsize=(12,12))
plt.plot(x_stars/(149*10**9),y_stars/(149*10**9), 'o', ms=0.1,color= 'white')
plt.plot([0],[0], 'o', ms=11, color ='#b32400')
plt.plot([0],[0], 'o', ms=9, color ='#262626')
plt.plot(X/(149*10**9),Y/(149*10**9), 'r', alpha = 0.05)
plt.xlabel('a.e.')
plt.ylabel('a.e.')

def func_anim(i):
    plt.plot(X[:i]/(149*10**9),Y[:i]/(149*10**9), 'green', linewidth = 0.6)
    
ani = animation.FuncAnimation(fig, func_anim, frames = np.arange(0,N,100), interval=1)

plt.axis('equal')
ani.save('Black_hole_field.gif')

# Вывод времени выполнения программы
duration1 = timer() - start1
print('Время выполнения =',duration1)