import numpy
from scipy.integrate import odeint
from matplotlib.pyplot import *

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


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]
data = numpy.loadtxt("Jupiter.txt",skiprows=5,delimiter=",",usecols=(2,3,4,5,6,7))
Y_jupiter = data[0]
Yi = numpy.append(Y_terre,Y_lune)
Yi = numpy.append(Yi,Y_jupiter)

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)
    rJ3 = numpy.power(Y[12]*Y[12]+Y[13]*Y[13]+Y[14]*Y[14],1.5)
    xTL = Y[6]-Y[0]
    yTL = Y[7]-Y[1]
    zTL = Y[8]-Y[2]
    xTJ = Y[12]-Y[0]
    yTJ = Y[13]-Y[1]
    zTJ = Y[14]-Y[2]
    xJL = Y[6]-Y[12]
    yJL = Y[7]-Y[13]
    zJL = Y[8]-Y[14]
    rTL3 = numpy.power(xTL*xTL+yTL*yTL+zTL*zTL,1.5)
    rTJ3 = numpy.power(xTJ*xTJ+yTJ*yTJ+zTJ*zTJ,1.5)
    rJL3 = numpy.power(xJL*xJL+yJL*yJL+zJL*zJL,1.5)
    xrT = Y[0]/rT3
    yrT = Y[1]/rT3
    zrT = Y[2]/rT3
    xrL = Y[6]/rL3
    yrL = Y[7]/rL3
    zrL = Y[8]/rL3
    xrJ = Y[12]/rJ3
    yrJ = Y[13]/rJ3
    zrJ = Y[14]/rJ3
    aT = -k2*((1.0+mT)*xrT+mL*(-xTL/rTL3+xrL)+mJ*(-xTJ/rTJ3+xrJ))
    bT = -k2*((1.0+mT)*yrT+mL*(-yTL/rTL3+yrL)+mJ*(-yTJ/rTJ3+yrJ))
    cT = -k2*((1.0+mT)*zrT+mL*(-zTL/rTL3+zrL)+mJ*(-zTJ/rTJ3+zrJ))
    aL = -k2*((1.0+mL)*xrL+mT*(xTL/rTL3+xrT)+mJ*(xJL/rJL3+xrJ))
    bL = -k2*((1.0+mL)*yrL+mT*(yTL/rTL3+yrT)+mJ*(yJL/rJL3+yrJ))
    cL = -k2*((1.0+mL)*zrL+mT*(zTL/rTL3+zrT)+mJ*(zJL/rJL3+zrJ))
    aJ = -k2*((1.0+mJ)*xrJ+mT*(xTJ/rTJ3+xrT)+mL*(-xJL/rJL3+xrL))
    bJ = -k2*((1.0+mJ)*yrJ+mT*(yTJ/rTJ3+yrT)+mL*(-yJL/rJL3+yrL))
    cJ = -k2*((1.0+mJ)*zrJ+mT*(zTJ/rTJ3+zrT)+mL*(-zJL/rJL3+zrL))
    return numpy.array([Y[3],Y[4],Y[5],aT,bT,cT,Y[9],Y[10],Y[11],aL,bL,cL,Y[15],Y[16],Y[17],aJ,bJ,cJ])

T=365*5
N=5000
t=numpy.linspace(0,T,N)
Y=odeint(systeme,Yi,t,rtol=1e-9,atol=1e-9)
xT=Y[:,0]
yT=Y[:,1]
zT=Y[:,2]
xL=Y[:,6]
yL=Y[:,7]
zL=Y[:,8]
xJ=Y[:,12]
yJ=Y[:,13]
zJ=Y[:,14]

figure()
axes().set_aspect('equal')
plot(xT,yT,"b-")
plot(xJ,yJ,"r-")
xlabel("x (ua)")
ylabel("y (ua)")
grid()
xlim(-6,6)
ylim(-6,6)

mJ=0
t=numpy.linspace(0,T,N)
Y=odeint(systeme,Yi,t,rtol=1e-9,atol=1e-9)
xT_sJ=Y[:,0]
yT_sJ=Y[:,1]
zT_sJ=Y[:,2]
xL_sJ=Y[:,6]
yL_sJ=Y[:,7]
zL_sJ=Y[:,8]

figure()
axes().set_aspect('equal')
plot(xT-xT_sJ,yT-yT_sJ,"b-")
xlabel("x (ua)")
ylabel("y (ua)")
title("Terre")
grid()

figure()
axes().set_aspect('equal')
plot(xL-xL_sJ,yL-yL_sJ,"k-")
xlabel("x (ua)")
ylabel("y (ua)")
title("Lune")
grid()


show()