from visual import *

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

# SET PARAMETERS=================
mass = 0.020 # of rolling object
a = 40.0   # surface inclination angle (degrees)
r1 = 0.022  # outer radius of rolling object (meters)
mom_inertia=(2.0/5.0)*mass*r1*r1 # moment of inertia
x = 1.0   # distance roll (meters)
h= 1.14 # height of end of inclined plane(meters)
air_f = 0.5 # air friction constant
#=================================
# Vpython x-y-z axes
#    y
# xxxxxxx
#    y
# Max=Nsin(a)-fcos(a)
# May=Ncos(a)-Mg+fsin(a)
# Rf=Icom(a/R)
i_fraction=mom_inertia/(mass*r1*r1)
grav = 9.81
pi=3.14159
air_density=1.29
dt = 0.01
r2d=180.0/3.14159
time = 0.0
t = 0.05 # surface thickness (meters)
scene.background=(1,1,1)
autoscale = 1
kvalue=0.5*air_f*air_density*pi*r1**2/mass
a = a*(3.14159/180.0) # convert to radians
t2 = t/2.0
xx=x*cos(a)
yy=x*sin(a)
incline = vector(xx,yy,0.0)
ramp = box(pos=(xx*0.5,yy*0.5,0.0), axis=incline, length=x, height=t, width=0.5, color=color.blue)
floor = box(pos=(-x*0.5,-h-(t/2.0),0), length=1.5*x, height=t, width=0.5, color=color.blue)
# starting position
z0 = 0.0
tr = 0.0
dx = (r1+t2)*sin(a)
x0 = x*cos(a)-dx
y0 = x*sin(a)+(r1+t2+tr)*cos(a)
pos0 = vector(x0,y0,z0)
roller = sphere(pos=pos0, radius=r1, color=color.red)
roller.velocity = vector(0,0,0)
# acceleration
con=1.0/(1.0+i_fraction)
acc1 = vector(-con*grav*sin(a)*cos(a),-con*grav*sin(a)*sin(a),0.0)
acc2 = vector(0.0,-grav,0.0)
scene.autoscale=0
scene.uniform=1
done=0
free=0
set=0
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
    if (free==0):
        roller.acceleration = acc1+vector(axf,ayf,0)
    else:
        roller.acceleration = acc2+vector(axf,ayf,0)        
    if (done==0):
        roller.velocity=roller.velocity+roller.acceleration*0.5*dt
        done=1
    else:
        roller.velocity=roller.velocity+roller.acceleration*dt
    roller.pos = roller.pos + roller.velocity*dt
    time=time+dt
    if roller.x < -dx:
        if (set==0):
            v0x=roller.velocity.x
            v0y=roller.velocity.y
            v0=(v0x*v0x+v0y*v0y)**0.5
            a0=atan(v0y/v0x)*r2d
            set=1
        free=1
        if roller.y < -h+r1:
            p1=-round((roller.pos.x+dx)*1000.0)/1000.0
            p2=-round(roller.velocity.x*1000.0)/1000.0
            print "x impact = ",p1, " meters  x-velocity = ",p2
            p1=round(v0*1000.0)/1000.0
            p2=round(a0*1000.0)/1000.0
            print "launch velocity = ",p1," m/s launch angle = ",-p2," degrees"
            roller.pos=pos0
            roller.velocity = vector(0.0,0.0,0.0)
            time = 0.0
            free=0
            done=0