import numpy as np
import matplotlib.pyplot as plt
import time

#----------------- declaracao de variaveis

dt = 0.1*60*60          # delta t (hora)
Mt = 5.9736*(10**24)    # Massa da Terra(kg)
Ml = 7.349*(10**22)     # Massa da Lua (kg)
Ms = 1.9891*(10**30)    # Massa do Sol (kg)
G  = 6.67418*(10**-11)  # Constante Gravitacional (m3kg-1s-2)
kt = 0                  # Variavel de controle afelio e perielio terra
kl = 0                  # Variavel de controle apogeu e perigeu lua
i  = 0                  # Variavel contadora dentro do while
l  = 0                  # Variavel contadora lua kl
tT = 0                  # Periodo da terra em dias

#----------------- condições iniciais (t=0)

t = 0
xt = [147098290000]                                                     #posicao inicial terra em x (m)
yt = [0]                                                                #posicao inicial terra em y (m)
xl = [363104000+xt[0]]                                                  #posicao inicial lua em x (m)
yl = [0]                                                                #posicao inicial lua em y (m)
xlt = [xl[0]-xt[0]]                                                     #posicao inicial lua em relacao a terra x (m)
ylt = [yl[0]-yt[0]]                                                     #posicao inicial lua em relacao a terra y (m)
vxt = [0]                                                               #velocidade inicial terra em x (m/s)
vyt = [30278]                                                           #velocidade inicial terra em y (m/s)
vxl = [0]                                                               #velocidade inicial lua em x (m/s)
vyl = [(3977.82/3.6)+vyt[0]]                                            #velocidade inicial lua em y (m/s)
rt = [((xt[0]**2)+(yt[0]**2))**(1/2)]                                   #raio da orbita terra (m)
rl = [((xl[0]**2)+(yl[0]**2))**(1/2)]                                   #raio da orbita lua (m)
rtl =[((xt[0]-xl[0])**2+(yt[0]-yl[0]))**(1/2)]                          #raio da orbita terra em relação a lua (m)
alx =[(-G*(((Ms/(rl[0]**3))*xl[0])+((Mt/(rtl[0]**3))*(xl[0]-xt[0]))))]  #aceleração da lua em x (m/s2)
aly =[(-G*(((Ms/(rl[0]**3))*yl[0])+((Mt/(rtl[0]**3))*(yl[0]-yt[0]))))]  #aceleração da lua em y (m/s2)
atx =[(-G*(((Ms/(rt[0]**3))*xt[0])+((Ml/(rtl[0]**3))*(xt[0]-xl[0]))))]  #aceleração da terra em x (m/s2)
aty =[(-G*(((Ms/(rt[0]**3))*yt[0])+((Ml/(rtl[0]**3))*(yt[0]-yl[0]))))]  #aceleração da terra em y (m/s2)
Xlt =[xl[0]-xt[0]]                                                      #Vetor X da lua em relacao a terra (para achar perielios e apogeus)
Ylt =[yl[0]-yt[0]]                                                      #Vetor Y da lua em relacao a terra (para achar perielios e apogeus)
Vxlt =[vxt[0]-vxl[0]]                                                   #Vetor velocidade x da lua em relacao a terra
Vylt =[vxt[0]-vyl[0]]                                                   #Vetor velocidade Y da lua em relacao a terra

# ----------------- evolução temporal ----

while i<365:
    tT= tT + (dt/(60*60*24))
    if i==0:
        vxt.append((vxt[i]+atx[i]*dt/2))              #velocidade terra em x1 (m/s)
        vyt.append((vyt[i]+aty[i]*dt/2))              #velocidade terra em y1 (m/s)
        vxl.append((vxl[i]+alx[i]*dt/2))              #velocidade lua em x1 (m/s)
        vyl.append((vyl[i]+aly[i]*dt/2))              #velocidade lua em y1 (m/s)
    else:
        vxt.append((vxt[i]+atx[i]*dt))                #velocidade terra em x (m/s)
        vyt.append((vyt[i]+aty[i]*dt))                #velocidade terra em y (m/s)
        vxl.append((vxl[i]+alx[i]*dt))                #velocidade lua em x (m/s)
        vyl.append((vyl[i]+aly[i]*dt))                #velocidade lua em y (m/s)

    xt.append(xt[i]+vxt[i+1]*dt)                      #posicao terra em x (m)
    yt.append(yt[i]+vyt[i+1]*dt)                      #posicao terra em y (m)
    xl.append(xl[i]+vxl[i+1]*dt)                      #posicao lua em x (m)
    yl.append(yl[i]+vyl[i+1]*dt)                      #posicao lua em y (m)
    xlt.append(xl[i+1]-xt[i+1])                       #posicao inicial lua em relacao a terra x (m)
    ylt.append(yl[i+1]-yt[i+1])                       #posicao inicial lua em relacao a terra x (m)

    rt.append(((xt[i+1]**2)+(yt[i+1]**2))**(1/2))                       #raio da orbita terra (m)
    rl.append(((xl[i+1]**2)+(yl[i+1]**2))**(1/2))                       #raio da orbita lua (m)
    rtl.append((((xt[i+1]-xl[i+1])**2+((yt[i+1]-yl[i+1]))**2))**(1/2))  #raio da orbita terra em relação a lua (m)

    alx.append((-G*(((Ms/(rl[i+1]*3))*xl[i+1])+((Mt/(rtl[i+1]**3))*(xl[i+1]-xt[i+1]))))) #aceleração da lua em x (m/s2)
    aly.append((-G*(((Ms/(rl[i+1]**3))*yl[i+1])+((Mt/(rtl[i+1]**3))*(yl[i+1]-yt[i+1]))))) #aceleração da lua em y (m/s2)
    atx.append((-G*(((Ms/(rt[i+1]**3))*xt[i+1])+((Ml/(rtl[i+1]**3))*(xt[i+1]-xl[i+1]))))) #aceleração da terra em x (m/s2)
    aty.append((-G*(((Ms/(rt[i+1]**3))*yt[i+1])+((Ml/(rtl[i+1]**3))*(yt[i+1]-yl[i+1]))))) #aceleração da terra em y (m/s2)

    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

    plt.show()

    print()
    i += 1