from pylab import *
from numpy import *

dt=0.001
N=3500

x=empty(N+1)
y=empty(N+1)
z=empty(N+1)
vx=empty(N+1)
vy=empty(N+1)
vz=empty(N+1)
ax=empty(N+1)
ay=empty(N+1)
az=empty(N+1)
ux=empty(N+1)
uy=empty(N+1)
uz=empty(N+1)
u=empty(N+1)

m=0.43
g=-9.81

wx = 0.
wy = 3.
wz = 0.

x[0]=5.
y[0]=0.
z[0]=0.

v0 = 30.
uhol = 40. # v stupňoch

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

for n in range(0,N):
    ux[n] = vx[n] - wx
    uy[n] = vy[n] - wy
    uz[n] = vz[n] - wz
    u[n]= sqrt( ux[n]*ux[n] + uy[n]*uy[n] + uz[n]*uz[n] )
    ax[n]= -0.006 * u[n] * ux[n] / m
    ay[n]= -0.006 * u[n] * uy[n] / m
    az[n]= (m*g - 0.006 * u[n] * uz[n])/m
    x[n+1]=x[n]+vx[n]*dt
    y[n+1]=y[n]+vy[n]*dt
    z[n+1]=z[n]+vz[n]*dt
    vx[n+1]=vx[n]+ax[n]*dt
    vy[n+1]=vy[n]+ay[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.11 and z[n+1] < 0.11): 
        print(z[n])
        print(x[n])
        print(y[n])