#Importando as Bibliotecas necessárias
import numpy as np # Biblioteca de vetores e matrizes
import matplotlib.pyplot as plt # Biblioteca de plotar gráficos
import matplotlib.pyplot as pltx # Biblioteca de plotar gráficos

#Constantes
nt = 100000 # Números de intervalos
dt = 360 # Delta_t(s) => 1h
mt = 5.9742*(10**24)  # Massa da Terra(Kg)
ml = 7.36*(10**22)  # Massa da Lua(Kg)
ms = 1.9891*(10**30)  # Massa do Sol(Kg)
g = 6.67428*(10**-20) # Constante gravitacional(SI) => km3/kg*s2

#tt = 31.557.600s #Período da terra em relação ao sol (s) => 365,25 dias => 31.557.600s
#tl= 2.360.880s #Período da lua em relação a terra (s) => 27,325 dias => 2.360.880s

#Criando vetores variáveis e zerando OBS: Sol referenciado fixo
t = np.zeros(nt)   # Tempo discretizado

xt = np.zeros(nt)  # Posição horizontal da Terra
yt = np.zeros(nt)  # Posição vertical da Terra
vxt = np.zeros(nt) # Velocidade horizontal da Terra
vyt = np.zeros(nt) # Velocidade vertical da Terra
axt = np.zeros(nt) # Aceleração horizontal da Terra
ayt = np.zeros(nt) # Aceleração vertical da Terra

rt = np.zeros(nt) # distância do Sol à Terra
rl = np.zeros(nt) # distância do Sol à Lua
rlt = np.zeros(nt) # distância da Lua à Terra

xl = np.zeros(nt)  # Posição horizontal da Lua
yl = np.zeros(nt)  # Posição vertical da Lua
vxl = np.zeros(nt) # Velocidade horizontal da Lua
vyl = np.zeros(nt) # Velocidade vertical da Lua
axl = np.zeros(nt) # Aceleração horizontal da Lua
ayl = np.zeros(nt) # Aceleração vertical da Lua

xlt = np.zeros(nt) # Posição horizontal da Lua em relação à Terra
ylt = np.zeros(nt) # Posição vertical da Lua em relação à Terra

#Condições iniciais -------------------
t[0] = 0

xt[0] = 147098290 #Km => periélio
yt[0] = 0

xl[0] = xt[0] + 363104 #Km (rt + rlt)
yl[0] = 0

xlt[0] = 363104 #Km
ylt[0] = 0

vxt[0] = 0
vyt[0] = 30.2787  #Km/s
vyl[0] = vyt[0] + 1.10495 #31,38365 #Km/s


#calculando distâncias em t=0
rt[0] = ( xt[0]**2 + yt[0]**2 )**0.5
rl[0] = ( xl[0]**2 + yl[0]**2 )**0.5
rlt[0] = ( (xt[0]-xl[0])**2 + (yt[0]-yl[0])**2 )**0.5

#calculando acelerações em t=0
axt[0] = -g*( (ms*xt[0]/rt[0]**3) + ml*(xt[0]-xl[0])/rlt[0]**3 )
axl[0] = -g*( (ms*xl[0]/rl[0]**3) + mt*(xl[0]-xt[0])/rlt[0]**3 )
ayt[0] = 0
ayl[0] = 0

#Fórmulas para a evolução temporal
for it in range (1,nt):
  t[it] = it*dt

  if it==1:
    #calculando velocidades
    vxt[1] = vxt[0] + (0.5*dt*axt[0])
    vyt[1] = vyt[0] + (0.5*dt*ayt[0])

    vxl[1] = vxl[0] + (0.5*dt*axl[0])
    vyl[1] = vyl[0] + (0.5*dt*ayl[0])
  else:
    #calculando velocidades
    vxt[it] = vxt[it-1] + (dt*axt[it-1])
    vyt[it] = vyt[it-1] + (dt*ayt[it-1])

    vxl[it] = vxl[it-1] + (dt*axl[it-1])
    vyl[it] = vyl[it-1] + (dt*ayl[it-1])

  #calculando posição (x,y)
  xt[it] = xt[it-1] + (vxt[it] * dt)
  xl[it] = xl[it-1] + (vxl[it] * dt)
  yt[it] = yt[it-1] + (vyt[it] * dt)
  yl[it] = yl[it-1] + (vyl[it] * dt)

  xlt[it] = xl[it] - xt[it]
  ylt[it] = yl[it] - yt[it]

  #calculando distâncias
  rt[it] = ( xt[it]**2 + yt[it]**2 )**0.5
  rl[it] = ( xl[it]**2 + yl[it]**2 )**0.5
  rlt[it] = ( (xt[it]-xl[it])**2 + (yt[it]-yl[it])**2 )**0.5

  #calculando aceleração
  axt[it] = -g*( (ms*xt[it]/rt[it]**3) + ml*(xt[it]-xl[it])/rlt[it]**3 )
  axl[it] = -g*( (ms*xl[it]/rl[it]**3) + mt*(xl[it]-xt[it])/rlt[it]**3 )

  ayt[it] = -g*( (ms*yt[it]/rt[it]**3) + ml*(yt[it]-yl[it])/rlt[it]**3 )
  ayl[it] = -g*( (ms*yl[it]/rl[it]**3) + mt*(yl[it]-yt[it])/rlt[it]**3 )

#---------------
print('afélio  : ', round(np.max(rt),0))
print('periélio: ', round(np.min(rt),0))
print('apogeu  : ', round(np.max(rlt),0))
print('perigeu : ', round(np.min(rlt),0))
print()

#determinando período da Lua => 2 * tempo para ylt tornar-se negativo
for it in range (1,nt):
  if( ylt[it] < 0 ):
    print('Período da Lua: ', round(2*it*dt /(24*60*60),2))
    break

#determinando período da Terra => 2 * tempo para yt tornar-se negativo
for it in range (1,nt):
  if( yt[it] < 0 ):
    print('Período da Terra: ', round(2*it*dt /(24*60*60),2))
    break
print()

# Plotando os Gráficos

# colocando um estilo para uso
plt.style.use('fivethirtyeight')

# criando uma figura
fig = plt.figure()

# definindo subplots e suas posições na figura
plt1 = fig.add_subplot(121)

plt1.set_title('$Terra/Lua - Sol$')
plt1.set_xlabel('x(km)')
plt1.set_ylabel('y(km)')

for i in range(0, nt, 20):
    # plotting points on each subplot
    plt1.scatter(0., 0., color='r', s=40) #sol = vermelho
    plt1.scatter(xt[i], yt[i], color='b', s=3) #terra = azul
    plt1.scatter(xl[i], yl[i], color='g', s=3) #lua = verde

# função para mostrar a plotagem
plt.show()

print()

# colocando um estilo para uso
pltx.style.use('fivethirtyeight')

# criando uma figura
fig = pltx.figure()

# definindo subplots e suas posições na figura
plt2 = fig.add_subplot(121)

plt2.set_title('$Lua - Terra$')
plt2.set_xlabel('xlt(km)')
plt2.set_ylabel('ylt(km)' )

for i in range(0, nt, 500):
    # plotting points on each subplot
    plt2.scatter(0., 0., color='b', s=40) #terra = azul
    plt2.scatter(xlt[i], ylt[i], color='g', s=3) #lua = verde

# função para mostrar a plotagem
pltx.show()