import sympy as sy
import numpy as np
import matplotlib.pyplot as plt

# Exercício 1. Método de Euler orbita da Terra 
# EULER:
#for i in range(N):

#   x[i+1] = x[i] + vx[i] * dt
#   vx[i+1] = vx[i] + aceler * dt
#   t[i+1] = t[i] + dt

#consts
AU = 1.498 * 10**11  # dm da Terra ao Sol (m)
ano = 3.15 * 10**7  # segundos num ano
M = 1.989 * 10**30  # massa do Sol (kg)
m = 5.9722 * 10**24  # massa da Terra (kg)
G = 6.67408 * 10**-11  # constante gravitacional (m^3/(kg*s^2))

r0 = (1 * AU, 0)    # posição inicial da Terra
v0 = (0, 2 * np.pi * AU / ano)  # velocidade inicial da Terra

# Fg = -G * (m * M * r) / (modR**3)  # aceleração gravitacional

########
#1
#Simule a órbita da Terra á volta do sol, usando o método de Euler
#trabalhe no sistema astronómico de unidades

G = 4 * np.pi**2  # (AU^3/(M*ano^2))
M = 1  # massa do sol (M)


r0 = np.array([1.0, 0.0])          
v0 = np.array([0.0, 2 * np.pi])

t0 = 0.0          
tf = 4.0           
dt = 0.01
n = int(tf / dt)

t = np.linspace(t0, tf, n)
v = np.empty((n,2))
r = np.empty((n,2))

r[0] = r0
v[0] = v0

#euler ##nao conserva energia
for i in range(n - 1):

    modR = np.sqrt(r[i][0]**2 + r[i][1]**2)
    Fg = -G * (m * M * r[i]) / (modR**3)
    a = Fg / m
    
    r[i+1] = r[i] + v[i] * dt
    v[i+1] = v[i] + a * dt
    t[i+1] = t[i] + dt


# Plot da orbita
plt.figure()
plt.plot(r[:, 0], r[:, 1], label="Órbita (Euler)")
plt.xlabel("x (UA)")
plt.ylabel("y (UA)")
plt.title("Simulação da Órbita da Terra em torno do Sol")
plt.grid()
plt.show()

#utilizando o metodo de euler, a orbita da terra nao é fechada, nem uma elipse