import numpy
from matplotlib.pyplot import *

k = 0.01720209895
k2 = k*k
mT = 1.0/332946
mL = 1.0/27068620.9


data = numpy.loadtxt("Terre.txt",skiprows=5,delimiter=",",usecols=(2,3,4,5,6,7))
Y_terre = data[0]               
data = numpy.loadtxt("Lune.txt",skiprows=5,delimiter=",",usecols=(2,3,4,5,6,7))
Y_lune = data[0]               
Yi = numpy.append(Y_terre,Y_lune)

def systeme(Y,t):
    rT3 = numpy.power(Y[0]*Y[0]+Y[1]*Y[1]+Y[2]*Y[2],1.5)
    rL3 = numpy.power(Y[6]*Y[6]+Y[7]*Y[7]+Y[8]*Y[8],1.5)
    xTL = Y[6]-Y[0]
    yTL = Y[7]-Y[1]
    zTL = Y[8]-Y[2]
    rTL3 = numpy.power(xTL*xTL+yTL*yTL+zTL*zTL,1.5)
    aT = -k2*((1.0+mT)*Y[0]/rT3+mL*(-xTL/rTL3+Y[6]/rL3))
    bT = -k2*((1.0+mT)*Y[1]/rT3+mL*(-yTL/rTL3+Y[7]/rL3))
    cT = -k2*((1.0+mT)*Y[2]/rT3+mL*(-zTL/rTL3+Y[8]/rL3))
    aL = -k2*((1.0+mL)*Y[6]/rL3+mT*(xTL/rTL3+Y[0]/rT3))
    bL = -k2*((1.0+mL)*Y[7]/rL3+mT*(yTL/rTL3+Y[1]/rT3))
    cL = -k2*((1.0+mL)*Y[8]/rL3+mT*(zTL/rTL3+Y[2]/rT3))
    return numpy.array([Y[3],Y[4],Y[5],aT,bT,cT,Y[9],Y[10],Y[11],aL,bL,cL])

def pas_euler(systeme,h,t,Yn):
    deriv = systeme(Yn,t)
    return Yn+h*deriv

def euler(systeme,Yi,T,h):
    Y = Yi
    t = 0.0
    liste_y = [Y]
    liste_t = [t]
    while t<T:
        Y = pas_euler(systeme,h,t,Y)
        t += h
        liste_y.append(Y)
        liste_t.append(t)
    return (numpy.array(liste_t),numpy.array(liste_y))

T=365
h=1e-1
(t,Y)=euler(systeme,Yi,T,h)
xT=Y[:,0]
yT=Y[:,1]
zT=Y[:,2]
xL=Y[:,6]
yL=Y[:,7]
zL=Y[:,8]
rTL = numpy.sqrt((xL-xT)**2+(yL-yT)**2+(zL-zT)**2)
rTS = numpy.sqrt(xT**2+yT**2+zT**2)

figure()
axes().set_aspect('equal')
plot(xT,yT,"k-")
xlabel("x (ua)")
ylabel("y (ua)")
grid()
xlim(-1.2,1.2)
ylim(-1.2,1.2)

figure()
axes().set_aspect('equal')
plot(xL-xT,yL-yT,"k-")
xlabel("x (ua)")
ylabel("y (ua)")
grid()

figure()
plot(t,rTS,"k-")
xlabel("t (j)")
ylabel("rTS")
grid()

figure()
plot(t,rTL,"k-")
xlabel("t (j)")
ylabel("rTL")
grid()

# mouvement de la Terre sans la Lune
mL=0
(t,Y)=euler(systeme,Yi,T,h)
xT_sl=Y[:,0]
yT_sl=Y[:,1]
zT_sl=Y[:,2]

figure()
axes().set_aspect('equal')
plot(xT-xT_sl,yT-yT_sl,"k-")
xlabel("x (ua)")
ylabel("y (ua)")
grid()


show()