# this version of the Euler-Cromer code should look quite like the 1-D version
# except that acceleration, velocity, and position are all defined to be 3-D
# vectors.  You'll notice in the line "funct1.plot(pos=(r.x,r.y))" that when you
# just want the x or y component of a vector, you call for it by appending
# ".x" or ".y" to the vector name.

from __future__ import division
from visual import *

from visual.graph import *

# set up for the graph
graphA = gdisplay(title="y vs. x",xtitle="x", ytitle="y")
funct0 = gcurve(gdisplay=graphA, color=color.cyan) #will be the no-drag trajectory plot
#funct1 = gcurve(gdisplay = graphA, color=color.red) #un-comment this line

# constant
D = 0.001               #in m, projectile diameter

gamma = 0.25            #in Ns^2/m^4, quadratic drag factor for sphere in air at STP
betta = 0.00016         #in Ns/m^2, linear drag factor for sphere in air at STP

rho = 1000              #in kg/m^3, density of the projectile, that of water
m = rho*(pi*D**3)/6     #in kg, projectile mass
gmag=9.8                #the magnitude of the acceleration due to gravity
g = vector(0,-gmag,0)   #the vector for acceleration due to gravity

# set initial conditions & time step
t=0
dt=0.002
ri = vector(0,0,0)   #initial position vector
vi = vector(15,15,0) #initial velocity vector

r0 = ri    #initializing the no-drag position vector            
v0 = vi    #initializing the no-drag velocity vector

#r1=ri #un-comment this line
#v1=vi #un-comment this line

bm = betta*6/(pi*rho*D**2) #b/m in terms of density and diameter
cm = gamma*6/(pi*rho*D)  #c/m in terms of density and diameter

# plot initial position
funct0.plot(pos=(r0.x,r0.y))
#fucnt1.plot(pos=(r1.x, r1.y)) #un-comment this line


# apply the Euler-Cromer method to each component
while not (r0.y < 0): #continue the calculation until the projectile lands
    a0 = g         #update no-drag acceleration vector
    #a1 = ...       # complete this line of code
    v0 = v0 + a0*dt #update no-drag velocity vector
    #v1 = ...       # complete this line of code
    r0 = r0 + v0*dt #update no-drag position vector
    #r1 = ...       # complete this line of code
    t=t+dt
    funct0.plot(pos=(r0.x,r0.y))  # add new points to the no-drag trajectory plot
    #funct1.plot(pos = (r1.x, r1.y)) #un-comment this line