Central field motion -motion of satellite around earth

from pylab import*

gm=10.0
x=20.0
y=0.0
#Satellite
cir=Circle((x,y),radius=2,  fc='g')
gca().add_patch(cir)
r=sqrt(x*x+y*y)
#velocity=sqrt(gm/r), we get circular path
# for less or more velocities we get elliptical or hyperbolic path
f=0.8 #choose different f values, for example f=0.8, 1.0, 1.2, 1.4 ,see difference
vx=0
vy=f*sqrt(gm/r)
h=0.05
xpos=[x]
ypos=[y]
#Earth
cir=Circle((0,0),radius=5,  fc='y')
gca().add_patch(cir)
 
for i in range(30000):
    r=(x*x+y*y)**0.5
    vx=vx-(gm*x*h/(r*r*r))
    vy=vy-(gm*y*h/(r*r*r))
    x=x+vx*h
    y=y+vy*h
    xpos.append(x)
    ypos.append(y)
plot(xpos,ypos,'.b')

axis('scaled')
xlabel('x')
ylabel('Y')
show()