from pylab import *
from numpy import *

dt=0.001
N=600

x=empty(N+1)
z=empty(N+1)
vx=empty(N+1)
vz=empty(N+1)
v=empty(N+1)
ax=empty(N+1)
az=empty(N+1)

m=0.057
g=-9.81

x[0]=0.
z[0]=1.5
vx[0] = 35.
vz[0] = 0.

for n in range(0,N):
    v[n]= sqrt( vx[n]*vx[n] + vz[n]*vz[n] )
    ax[n]= -0.001 * v[n] * vx[n] / m
    az[n]= (m*g - 0.001 * v[n] * vz[n])/m
    x[n+1]=x[n]+vx[n]*dt
    z[n+1]=z[n]+vz[n]*dt
    vx[n+1]=vx[n]+ax[n]*dt
    vz[n+1]=vz[n]+az[n]*dt

plot(x,z)

for n in range(0,N):
    if(x[n] <= 11.89 and x[n+1] > 11.89):
        print(n)
        print(x[n])
        print(z[n])
        
for n in range(0,N):
    if(z[n] >= 0.034 and z[n+1] < 0.034):
        print(n)
        print(z[n])
        print(x[n])