Curling in Vpython code

from __future__ import division

from visual import *

from visual.graph import*

#Parameters of surface

icebase = box(pos=(0,0,0), size = (40,260,.00125), color=color.white)

ring1b = cylinder(pos=(0,95,.49), axis=(0,0,.00125), radius=17, color=color.blue)

ring2b = cylinder(pos=(0,-95,.49), axis=(0,0,.00125), radius=17, color=color.blue)

ring3w = cylinder(pos=(0,95,.51), axis=(0,0,.00125), radius=12 , color=color.white)

ring4w = cylinder(pos=(0,-95,.51), axis=(0,0,.00125), radius=12, color=color.white)

ring5r = cylinder(pos=(0,95,.53), axis=(0,0,.00125), radius=7, color=color.red)

ring6r = cylinder(pos=(0,-95,.53), axis=(0,0,.00125), radius=7,color=color.red)

ring7w = cylinder(pos=(0,95,.55), axis=(0,0,.00125), radius=2, color=color.white)

ring8w = cylinder(pos=(0,-95,.55), axis=(0,0,.00125), radius=2, color=color.white)

hogline1 = box(pos=(0,-55,.55), size= (40,1.3,.00125), color=color.black)

hogline2 = box(pos=(0,55,.55), size= (40,1.3,.00125), color=color.black)

tline1= box(pos=(0,0,.61), size=(.8,250,.00125), color=color.black)

tline2= box(pos=(-5,0,.61), size=(.8,157.59,.00125), color=color.black)

tline3= box(pos=(5,0,.61), size=(.8,157.59,.00125), color=color.black)

tline4= box(pos=(0,95,.63), size=(40,.8,.00125), color=color.black)

tline4= box(pos=(0,-95,.63), size=(40,.8,.00125), color=color.black)

backline1=box(pos=(0,112.5,.65), size=(40,.8,.00125),color=color.black)

backline2=box(pos=(0,-112.5,.65), size=(40,.8,.00125),color=color.black)

#Rocks

rock1=sphere(pos=(0,-1,.67), radius=2.5, color=color.yellow)

rock2=sphere(pos=(0,85,.67),radius=2.5, color=color.cyan)

rock3=sphere( pos=(0,30,.67), radius=2.5, color=color.cyan)

rock4=sphere( pos=(0,70,.67), radius=2.5, color=color.yellow)

rock1.velocity=vector(0,25.5,0)

rock2.velocity=vector(0,0,0)

rock3.velocity=vector(0,0,0)

rock4.velocity=vector(0,0,0)

#constants of simulation

coefficient_friction=.4

m=1

g=-9.8

time=0

dt=.001

#friction

k = .5*m*(pow(mag(rock1.velocity),2)+pow(mag(rock2.velocity),2)+pow(mag(rock3.velocity),2)+pow(mag(rock4.velocity),2))

#kinetic energy

p = m*(rock1.velocity.y+rock2.velocity.y+rock3.velocity.y+rock4.velocity.y)

graph1 = gdisplay(x=650, y=0, width=600, height=400, title='Momentum vs. Time',

xtitle='time', ytitle='Momentum', xmax=20, xmin=0,

ymax=1*p, ymin=-1*p, foreground=color.black, background=color.white)

pcurve = gcurve(gdisplay = graph1, color=color.black)

#momentum graph

graph2 = gdisplay(x=650, y=400, width=600, height=400, title='Kinetic Energy vs. Time',

xtitle='time', ytitle='Kinetic Energy', xmax=20, xmin=0,

ymax=1.2*k, ymin=0, foreground=color.black, background=color.white)

kcurve = gcurve(gdisplay = graph2, color=color.black)

#kinetic energy graph

while true:

rate(1000)

f_friction=m*g*(coefficient_friction)

a=g*(coefficient_friction)

time += dt

if rock1.velocity.y>0:

rock1.velocity.y=rock1.velocity.y+a*dt

if rock2.velocity.y>0:

rock2.velocity.y=rock2.velocity.y+a*dt

if rock3.velocity.y>0:

rock3.velocity.y=rock3.velocity.y+a*dt

if rock4.velocity.y>0:

rock4.velocity.y=rock4.velocity.y+a*dt

if rock1.velocity.y<0:

rock1.velocity.y=rock1.velocity.y-a*dt

if rock2.velocity.y<0:

rock2.velocity.y=rock2.velocity.y-a*dt

if rock3.velocity.y<0:

rock3.velocity.y=rock3.velocity.y-a*dt

if rock4.velocity.y<0:

rock4.velocity.y=rock4.velocity.y-a*dt

if rock3.pos.y>=rock4.pos.y-2*rock1.radius:

rock4.velocity.y= rock3.velocity.y

rock3.velocity.y=0

if rock4.pos.y>=rock2.pos.y-2*rock1.radius:

rock2.velocity.y= rock4.velocity.y

rock4.velocity.y=0

if rock1.pos.y>=rock3.pos.y-2*rock1.radius:

rock3.velocity.y=rock1.velocity.y

rock1.velocity.y=0

if rock2.pos.y<=rock4.pos.y-2*rock1.radius:

rock4.velocity.y=rock2.velocity.y

rock2.velocity.y=0

if rock2.pos.y>icebase.pos.y+.5*icebase.size.y-rock2.radius:

rock2.velocity = -rock2.velocity

rock1.pos=rock1.pos+rock1.velocity*dt

rock2.pos=rock2.pos+rock2.velocity*dt

rock3.pos=rock3.pos+rock3.velocity*dt

rock4.pos=rock4.pos+rock4.velocity*dt

p = m*(rock1.velocity.y+rock2.velocity.y+rock3.velocity.y+rock4.velocity.y)

k = k = .5*m*(pow(mag(rock1.velocity),2)+pow(mag(rock2.velocity),2)+pow(mag(rock3.velocity),2)+pow(mag(rock4.velocity),2))

if true:

pcurve.plot(pos = (time,p))

kcurve.plot(pos = (time,k))