from pylab import *

dt=1
N=100000

g=9.81
Rz=6371000.

x=empty(N+1)
y=empty(N+1)
vx=empty(N+1)
vy=empty(N+1)
r=empty(N+1)
ax=empty(N+1)
ay=empty(N+1)

x[0]=0.
y[0]=10. * Rz
vx[0] = 1900. 
vy[0] = 0.

for n in range(0,N):
    r[n]=sqrt(x[n]*x[n] + y[n]*y[n])
    ax[n]=-g*Rz*Rz * x[n]/(r[n]*r[n]*r[n])
    ay[n]=-g*Rz*Rz * y[n]/(r[n]*r[n]*r[n])
    
    x[n+1]=x[n]+vx[n]*dt
    y[n+1]=y[n]+vy[n]*dt
    vx[n+1]=vx[n]+ax[n]*dt
    vy[n+1]=vy[n]+ay[n]*dt

plot(x,y)
axis('square')
xlim([-3*y[0],3*y[0]])
ylim([-3*y[0],3*y[0]])
gca().add_patch(plt.Circle((0, 0), Rz)) #zemegula

for n in range(0,N-1):
    if(r[n] >= Rz and r[n+1] < Rz): 
        print("čas dopadu =", n*dt, "s")
        
for n in range(0,N-1):
    if(x[n] < 0 and x[n+1] >= 0): 
        print("čas obehu =", n*dt, "s")