from visual import *
from visual.graph import *
#adapted by Jon Marr, 9-March-2015

launchvel=vector(0,0,0)
# velocity of Rocket needed to circle Moon and return to Earth

#set fixed field of view by creating two, essentially invisible, dots at corners
dot = sphere(pos=(-9e8,-9e8,0),radius=0.01,color=color.black)
dot2 = sphere(pos=(9e8,9e8,0),radius=0.01,color=color.black)
scaler = 5
Earthr = 6.37e6
Moonr = 1.75e6

Earth = sphere(pos=(0,0,0), radius=Earthr*scaler, color=color.blue)
Moon = sphere(pos=(4e8,0,0), radius=Moonr*scaler, color=color.yellow)
Rocket = sphere(pos=(6.3701e6,0,0), radius=1.5e6*scaler, color=color.white)
Earth.velocity = vector(0,0,0)
Moon.velocity = vector(0,1000,0)
Rocket.velocity = launchvel
#initial positions and velocities

Earth.m=5.97e24 # mass of Earth in kg
Moon.m=7e22 # mass of Moon
Rocket.m=173 # mass of Ranger satellite
Earth.p=Earth.m*Earth.velocity
Moon.p=Moon.m*Moon.velocity
Rocket.p=Rocket.m*Rocket.velocity
# intial momenta of all bodies

G=6.67e-11 # Newton's Grav. constant in mks units

dt=1 #(time steps of 10 seconds)
t=0 # t is time in seconds
thr=0  # time in hours (to use for ease of calculating end of while loop)
while (thr < 24*120) and mag(Rocket.pos)>=Earthr: #watch for 4 months
    rate(6000)
    t=t+dt
    thr=t/3600 # time in hours = time in seconds/(60*60)

    #body1 = Earth; body2 = Moon; body3 = Ranger
    # first calculate force between bodies 1 and 2
    deltarEM=Earth.pos-Moon.pos # deltarEM = rE - rM
    #print mag(deltarEM)
    FEM = -(G*Earth.m*Moon.m/mag(deltarEM)**3)*deltarEM
    FME = -FEM

    # now calculate forces on Ranger
    deltarRE=Rocket.pos-Earth.pos 
    #print deltarRE.x,deltarRE.y,deltarRE.z
    FRE = -(G*Earth.m*Rocket.m/mag(deltarRE)**3)*deltarRE
    deltarRM=Rocket.pos-Moon.pos
    #print deltarRM.x,deltarRM.y,deltarRM.z
    FRM = -(G*Moon.m*Rocket.m/mag(deltarRM)**3)*deltarRM

    if mag(deltarRM) < 5e7 or mag(deltarRE) < 5e7:
        dt = 1
    else:
        dt = 10

    # now do p update using Net force on each body
    Earth.p = Earth.p + FEM*dt
    Moon.p = Moon.p + FME*dt
    Rocket.p = Rocket.p + (FRE+FRM)*dt

    # pos. update for each body
    Earth.pos = Earth.pos + (Earth.p/Earth.m)*dt
    Moon.pos = Moon.pos + (Moon.p/Moon.m)*dt
    Rocket.pos = Rocket.pos + (Rocket.p/Rocket.m)*dt

    if mag(Rocket.pos)<Earthr:
        print "CRASH to EARTH"

print "trip time (days)=",t/(24*3600)