from pylab import *
from numpy import *

dt=0.001
N=15000

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=4.
g=-1.62

x[0]=0.
z[0]=1.8

v0 = 14.4
uhol = 44. # v stupňoch

theta = uhol * pi/180.
vx[0] = v0 * cos(theta)
vz[0] = v0 * sin(theta)

for n in range(0,N):
    v[n]= sqrt( vx[n]*vx[n] + vz[n]*vz[n] )
    ax[n]= -0.002 * v[n] * vx[n] / m
    az[n]= (m*g - 0.002 * 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)
axis('square')

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