from __future__ import division, print_function
from vpython import *
scene.width = scene.height = 600

sff=2e5
sf=0.5
G = 6.7e-11
mEarth = 6e24
mcraft = 15e3
mMoon = 7e22
deltat = 10

Earth = sphere(pos=vector(0,0,0), radius=6.4e6, color=color.cyan)
craft = sphere(pos=vector(-10*Earth.radius, 0,0), radius=1e6,
               color=color.yellow, make_trail=True)
Moon = sphere(pos=vector(4.0e8,0,0), radius=1.75e6, color=color.white)
# elliptical orbit (0, 3.3e3,0)
# crashing on the moon (0,3.272e3,0)
# figure 8, almost (0,3.2735e3,0)
vcraft = vector(0,3.2735e3,0)
pcraft = mcraft*vcraft
print('p=', pcraft)
t = 0
#Fgrav = 1
parr = arrow(pos=craft.pos, axis=sf*pcraft, color=color.green, shaftwidth = craft.radius*2)#
#farr = arrow(pos=craft.pos, axis=sff*Fgrav, color=color.blue, shaftwidth = craft.radius/2)
farr=arrow(pos=craft.pos, color=color.red, shaftwidth=craft.radius*2)
#farrm=arrow(pos=craft.pos, color=color.red, shaftwidth=craft.radius/2)

#scene.autoscale = False ##turn off automatic camera zoom

scene.center = (Earth.pos + Moon.pos)/2

while t < 100*24*60*60:
    rate(5e3)     
    rvec=craft.pos-Earth.pos
    rvecm=craft.pos-Moon.pos
    if mag(rvecm) < Moon.radius:
        break
    rhat=norm(rvec)
    rhatm=norm(rvecm)
    fmag=G*mEarth*mcraft/mag(rvec)**2
    fmagm=G*mMoon*mcraft/mag(rvecm)**2
#    print('Fmag = ',fmag)
    Fgrav=fmag*rhat+fmagm*rhatm
#    Fgrav=fmag*rhat
    pcraft=pcraft-Fgrav*deltat
    craft.pos = craft.pos + (pcraft/mcraft)*deltat
    parr.pos = craft.pos
    parr.axis=sf*pcraft*10
    farr.pos = craft.pos
    farr.axis = -sff*Fgrav    
    t = t+deltat