# -*- coding: utf-8 -*-
"""
Created on Thu Feb 28 15:26:00 2019

@author: Rui
"""
import sys
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(threshold=sys.maxsize)
##########################
minput=input('Do you want the moon to exist in this model? y/n: ')
##########################
m=1
h=100
n=int(1000000/h)
ptno=list(range(0,n+1))
x=np.zeros(n+1)
y=np.zeros(n+1)
ycomb=np.zeros(n+1)
vx=np.zeros(n+1)
vy=np.zeros(n+1)
mx=np.zeros(n+1)
my=np.zeros(n+1)

xm=0
ym=384400000
Mm=7.35E22
Rm=1737100

x[0]=0
y[0]=-42164000
vx[0]=3075
vy[0]=0

for i in range(0,n):
    ycomb[i]=y[i]+ym

R=6371
M=5.972E24
G=6.67408E-11
################################################################################
f3=(-G*M*x)/((x**2+y**2)**(3/2))
f4=(-G*M*y)/((x**2+y**2)**(3/2))
f3m=((-G*Mm*x)/(((x)**2+np.power(ycomb,2))**(3/2)))+((-G*M*x)/((x**2+y**2)**(3/2)))
f4m=((-G*Mm*(ycomb))/(((x)**2+np.power(ycomb,2))**(3/2)))+((-G*M*y)/((x**2+y**2)**(3/2)))

###############################################################################
k1x=np.zeros(n+1)
k1y=np.zeros(n+1)
k1vx=np.zeros(n+1)
k1vy=np.zeros(n+1)
k2x=np.zeros(n+1)
k2y=np.zeros(n+1)
k2vx=np.zeros(n+1)
k2vy=np.zeros(n+1)
k3x=np.zeros(n+1)
k3y=np.zeros(n+1)
k3vx=np.zeros(n+1)
k3vy=np.zeros(n+1)
k4x=np.zeros(n+1)
k4y=np.zeros(n+1)
k4vx=np.zeros(n+1)
k4vy=np.zeros(n+1)
###############################################################################
if minput=='n':
    for i in range (0,n):
        k1x[i]=vx[i]
        k1y[i]=vy[i]
        k1vx[i]=(-G*M*x[i])/((x[i]**2+y[i]**2)**(3/2))
        k1vy[i]=(-G*M*y[i])/((x[i]**2+y[i]**2)**(3/2))

        k2x[i]=(vx[i]+(h*k1vx[i]/2))
        k2y[i]=(vy[i]+(h*(k1vy[i]/2)))
        k2vx[i]=(-G*M*(x[i]+(h*(k1x[i]/2)))/((x[i]+h*(k1x[i]/2))**2+(y[i]+(h*(k1y[i]/2)))**2)**(3/2))
        k2vy[i]=(-G*M*(y[i]+(h*(k1y[i]/2)))/((x[i]+h*(k1x[i]/2))**2+(y[i]+(h*(k1y[i]/2)))**2)**(3/2))

        k3x[i]=(vx[i]+(h*k2vx[i]/2))
        k3y[i]=(vy[i]+(h*k2vy[i]/2))
        k3vx[i]=(-G*M*(x[i]+h*(k2x[i]/2))/(((x[i]+h*(k2x[i]/2))**2+(y[i]+h*(k2y[i]/2))**2)**(3/2)))
        k3vy[i]=(-G*M*(y[i]+h*(k2y[i]/2))/(((x[i]+h*(k2x[i]/2))**2+(y[i]+h*(k2y[i]/2))**2)**(3/2)))

        k4x[i]=(vx[i]+(h*k3vx[i]))
        k4y[i]=(vy[i]+(h*k3vy[i]))
        k4vx[i]=(-G*M*(x[i]+(h*(k3x[i]))/(((x[i]+h*(k3x[i])**2+(y[i]+(h*(k3y[i])**2)**(3/2)))))))
        k4vy[i]=(-G*M*(y[i]+(h*(k3y[i]))/(((x[i]+h*(k3x[i])**2+(y[i]+(h*(k3y[i])**2)**(3/2)))))))

        x[i+1]=x[i]+(h/6)*(k1x[i]+2*k2x[i]+2*k3x[i]+k4x[i])
        y[i+1]=y[i]+(h/6)*(k1y[i]+2*k2y[i]+2*k3y[i]+k4y[i])
        vx[i+1]=vx[i]+(h/6)*(k1vx[i]+2*k2vx[i]+2*k3vx[i]+k4vx[i])
        vy[i+1]=vy[i]+(h/6)*(k1vy[i]+2*k2vy[i]+2*k3vy[i]+k4vy[i])

        if (np.sqrt(x[i]**2+y[i]**2)<R):
            break

    en=int((2*np.pi)/0.1)
    ex=np.zeros(en)
    ey=np.zeros(en)

    for i in range(en):
        ex[i]=R*np.cos(i*h)
        ey[i]=R*np.sin(i*h)

    plt.plot(ex,ey)
    plt.plot(x,y)
    plt.axis('equal')
    plt.show()

    v=np.sqrt((vx**2)+(vy**2))
    r=np.sqrt((x**2)+(y**2))

    ke=(1/2)*m*(v**2)
    pe=(-G*M*m)/r
    etot=ke+pe

    plt.plot(ptno,ke)
    plt.plot(ptno,pe)
    plt.plot(ptno,etot)
#############################################################################
#if minput=='y':
#    for i in range (0,n):
#        k1x=f1(x[i],y[i],vx[i],vy[i])
#        k1y=f2(x[i],y[i],vx[i],vy[i])
#        k1vx=f3m(x[i],y[i],vx[i],vy[i])
#        k1vy=f4m(x[i],y[i],vx[i],vy[i])
#
#        k2x=f1(x[i],y[i],(vx[i]+(h*k1vx/2)),vy[i])
#        k2y=f2(x[i],y[i],vx[i],(vy[i]+(h*(k1vy/2))))
#        k2vx=f3m((x[i]+h*(k1x/2)),(y[i]+(h*(k1y/2))),vx[i],vy[i])
#        k2vy=f4m((x[i]+h*(k1x/2)),(y[i]+(h*(k1y/2))),vx[i],vy[i])
#
#        k3x=f1(x[i],y[i],(vx[i]+(h*k2vx/2)),vy[i])
#        k3y=f2(x[i],y[i],vx[i],(vy[i]+(h*k2vy/2)))
#        k3vx=f3m((x[i]+h*(k2x/2)),(y[i]+h*(k2y/2)),vx[i],vy[i])
#        k3vy=f4m((x[i]+h*(k2x/2)),(y[i]+h*(k2y/2)),vx[i],vy[i])
#
#        k4x=f1(x[i],y[i],(vx[i]+(h*k3vx)),vy[i])
#        k4y=f2(x[i],y[i],(vy[i]+(h*k3vy)),vy[i])
#        k4vx=f3m((x[i]+h*(k3x)),(y[i]+(h*(k3y))),vx[i],vy[i])
#        k4vy=f4m((x[i]+h*(k2x)),(y[i]+(h*(k2y))),vx[i],vy[i])
#
#        x[i+1]=x[i]+(h/6)*(k1x+2*k2x+2*k3x+k4x)
#        y[i+1]=y[i]+(h/6)*(k1y+2*k2y+2*k3y+k4y)
#        vx[i+1]=vx[i]+(h/6)*(k1vx+2*k2vx+2*k3vx+k4vx)
#        vy[i+1]=vy[i]+(h/6)*(k1vy+2*k2vy+2*k3vy+k4vy)
#
#        if (np.sqrt(x[i]**2+y[i]**2)<R):
#            break
#
#    for i in range(en):
#        mx[i]=xm+R*np.cos(i*h)
#        my[i]=ym+R*np.sin(i*h)
#
#    plt.plot(ex,ey)
#    plt.plot(x,y)
#    plt.plot(mx,my)
#    plt.axis('equal')
#    plt.show()
#
#    mx=np.zeros(en)
#    my=np.zeros(en)
#
#    for i in range(en):
#        mx[i]=xm+R*np.cos(i*h)
#        my[i]=ym+R*np.sin(i*h)