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]             
Xi=[Y_terre[0],Y_terre[1],Y_terre[2],Y_lune[0],Y_lune[1],Y_lune[2]]
Vi=[Y_terre[3],Y_terre[4],Y_terre[5],Y_lune[3],Y_lune[4],Y_lune[5]]

def accel(X):
    rT3 = numpy.power(X[0]*X[0]+X[1]*X[1]+X[2]*X[2],1.5)
    rL3 = numpy.power(X[3]*X[3]+X[4]*X[4]+X[5]*X[5],1.5)
    xTL = X[3]-X[0]
    yTL = X[4]-X[1]
    zTL = X[5]-X[2]
    rTL3 = numpy.power(xTL*xTL+yTL*yTL+zTL*zTL,1.5)
    aT = -k2*((1.0+mT)*X[0]/rT3+mL*(-xTL/rTL3+X[3]/rL3))
    bT = -k2*((1.0+mT)*X[1]/rT3+mL*(-yTL/rTL3+X[4]/rL3))
    cT = -k2*((1.0+mT)*X[2]/rT3+mL*(-zTL/rTL3+X[5]/rL3))
    aL = -k2*((1.0+mL)*X[3]/rL3+mT*(xTL/rTL3+X[0]/rT3))
    bL = -k2*((1.0+mL)*X[4]/rL3+mT*(yTL/rTL3+X[1]/rT3))
    cL = -k2*((1.0+mL)*X[5]/rL3+mT*(zTL/rTL3+X[2]/rT3))
    return numpy.array([aT,bT,cT,aL,bL,cL])

def pas_eulerA(accel,h,t,Xn,Vn):
    X = Xn+h*Vn
    A = accel(X)
    V = Vn+h*A
    return (X,V)

def eulerA(accel,Xi,Vi,T,h):
    X = numpy.array(Xi)
    V = numpy.array(Vi)
    t = 0.0
    liste_t = [t]
    liste_x = [X]
    liste_v = [V]
    while t<T:
        (X,V) = pas_eulerA(accel,h,t,X,V)
        t += h
        liste_t.append(t)
        liste_x.append(X)
        liste_v.append(V)
    return (numpy.array(liste_t),numpy.array(liste_x),numpy.array(liste_v))


T=365
h=1e-1
(t,X,V) = eulerA(accel,Xi,Vi,T,h)
xT=X[:,0]
yT=X[:,1]
zT=X[:,2]
xL=X[:,3]
yL=X[:,4]
zL=X[:,5]
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,X,V) = eulerA(accel,Xi,Vi,T,h)
xT_sl=X[:,0]
yT_sl=X[:,1]
zT_sl=X[:,2]


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


show()