from visual import *

# vpython model of an object rolling on an
# inclined surface without slipping

# SET PARAMETERS=================
mass=0.020 # mass of rolling object in kg
shape = 'sphere' #solid sphere
#shape = 'cylinder' #solid cylinder
#shape = 'ring'
air_friction = 0.5 # Air friction constant
a = 2.0   # surface inclination angle (degrees)
r1 = 0.12  # outer radius (meters)
r2 = 0.0  # inner radius (meters)
length = 0.0 # Cylinder length or ring width
x = 1.0   # roll distance (meters)
#=================================

air_density=1.29
pi=3.14159
scene.background=(1,1,1)
autoscale = 1
a = a*(3.14159/180.0) # convert to radians
t = 0.05 # surface thickness (meters)
t2 = t/2.0
incline = vector(x*cos(a),x*sin(a),0.0)
floor = box(axis=incline, length=x, height=t, width=0.5, color=color.blue)
# starting position
z0 = 0.0
tr = 0.0
if shape=='cylinder':
    z0 = -0.15
if shape=='ring':
    tr = t2
x0 = (x/2.0)*cos(a)-(r1+t2)*sin(a)
y0 = (x/2.0)*sin(a)+(r1+t2+tr)*cos(a)
pos0 = vector(x0,y0,z0)
if shape=='sphere':
  roller = sphere(pos=pos0, radius=r1, color=color.red)
  con=5.0/7.0
  area=pi*r1*r1
elif shape=='cylinder':
  roller = cylinder(pos=pos0, axis=(0,0,0.3),radius=r1, color=color.red)
  con=2.0/3.0
  area=2.0*r1*length
elif shape=='ring':
  roller = ring(pos=pos0, axis=(0,0,0.3),radius=r1, thickness=(r1-r2), color=color.red)
  con = 1.0/(1.5+0.5*(r2/r1)**2)
  area=2.0*r1*length
else:
  print "Invalid Choice ",shape
kvalue=air_friction*air_density*area/mass
roller.velocity = vector(0,0,0)
# acceleration
acc = vector(-con*9.8*sin(a)*cos(a),-con*9.8*sin(a)*sin(a),0.0)
dt = 0.01
time = 0.0
scene.autoscale=0
scene.uniform=1
while 1:
    rate(100)
    vmag2=(roller.velocity.x**2+roller.velocity.y**2)
    vmag=vmag2**0.5
    axf=-kvalue*vmag*roller.velocity.x
    ayf=-kvalue*vmag*roller.velocity.y
    roller.acceleration = acc+vector(axf,ayf,0)
    time = time + dt
    roller.pos = roller.pos + roller.velocity*dt
    if roller.pos.x < -((x/2.0)*cos(a)+(r1+t2)*sin(a)):
        roller.pos = pos0
        roller.velocity = vector(0.0,0.0,0.0)
        print "Time of impact = ",time, " seconds"
        time = 0.0
    else:
        roller.velocity = roller.velocity + roller.acceleration*dt