from visual import *
from math import *
scene.autoscale=0

scene.background=color.black


k=2.0
l=4.0
r=10.0
v=0.0
ball=sphere(color=color.yellow,radius=1.0)
ball.track=curve(color=color.red)
ball.pos.x=r
ball.pos.y=0.0
ball.pos.z=-k*(1-2/r)**(-0.5)

def deltatheta(r):
        return(l/(r*r))

def force(r):
        dum=r*r
        return(-1.0/dum+l*l/(dum*r)-3l*l/(dum*dum))

l1=curve(color=color.green)
l2=curve(color=color.green)
l3=curve(color=color.green)
l4=curve(color=color.green)
l5=curve(color=color.green)
l6=curve(color=color.green)
l7=curve(color=color.green)
l8=curve(color=color.green)
c1=curve(color=color.green)
c2=curve(color=color.green)
c3=curve(color=color.green)
c4=curve(color=color.green)
c5=curve(color=color.green)
c6=curve(color=color.green)
c7=curve(color=color.green)
c8=curve(color=color.green)



for r in arange(2.05,20.0,0.1):
        z=-k*(1-2/r)**(-0.5)
        theta=0.0
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l1.append(pos=temp)
        theta=2*math.pi/8
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l2.append(pos=temp)
        theta=4*math.pi/8
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l3.append(pos=temp)
        theta=6*math.pi/8
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l4.append(pos=temp)
        theta=8*math.pi/8
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l5.append(pos=temp)
        theta=10*math.pi/8
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l6.append(pos=temp)
        theta=12*math.pi/8
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l7.append(pos=temp)
        theta=14*math.pi/8
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        l8.append(pos=temp)

for theta in arange(0,2*math.pi,0.02):
        r=2.05
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c1.append(pos=temp)
        r=2.1
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c2.append(pos=temp)
        r=2.25
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c3.append(pos=temp)
        r=2.5
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c4.append(pos=temp)
        r=5.0
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c5.append(pos=temp)
        r=10.0
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c6.append(pos=temp)
        r=15.0
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c7.append(pos=temp)
        r=20.0
        z=-k*(1-2/r)**(-0.5)
        temp=vector(r*math.cos(theta),r*math.sin(theta),z)
        c8.append(pos=temp)
        
        
dt=0.01
t=0


while(1):
        #rate(500)
    r=r+dt*v
    v=v+dt*force(r)
    theta=theta+dt*deltatheta(r)
    ball.pos.x=r*math.cos(theta)
    ball.pos.y=r*math.sin(theta)
    ball.pos.z=-k*(1-2/r)**(-0.5)
    ball.track.append(pos=ball.pos)