Someone

from visual import*
from visual.graph import* 

#In Xanadu did Kublai Khan a stately pleasure-dome decree#

# Kinematic graphs
fg_color = color.white
bg_color = color.black
gc_color = color.red

#velocity graph
vel_graph = gdisplay(x=610, y=100, width=300, height=200,title= 'Car velocity vs. time',
                     xtitle='t (s)', ytitle='v (m/s)', xmax= 6., xmin= 0., ymax=3., ymin=-3.,
                     foreground=fg_color, background=bg_color)
vel_Plot = gcurve(color=gc_color)

#position graph
pos_graph = gdisplay(x=610, y=300, width=300, height=200, title= 'Car position',
                     xtitle='x position (m)', ytitle='y position (m)', xmax= 2.5, xmin= -1, ymax=2.5,
                     ymin= 0., foreground=fg_color, background=bg_color)
pos_Plot = gcurve (color=gc_color)

#Work graph
work_graph = gdisplay(x=610, y=600, width=300, height=200, title='Energy Loss to Friction vs. Time', xtitle='t(s)',
                      ytitle='Work/mass', xmax= 2., xmin=0., ymax=10., ymin=0., foreground=fg_color,
                      background=bg_color)
work_plot = gcurve(color=gc_color)

#input_parameters
sign=1.
mu = .025         #coeffecient of friction
R = .38   #radius of track
r = .01          #center of mass distance
R=R-r

g = 9.8     #gravity
h = .37      #starting height
k = 0       #traversal number

car = sphere(pos=(0,h,0), radius=0.01, color=color.green)
car.trail = curve(color=car.color)

#initial conditions
#theta = acos(.26/.34)
theta = asin(.26/.34)
omega = 0.
alfa = -g*sin(theta)/R - sign*mu*(omega*omega+g*cos(theta)/R)
car.vel = vector(0.,0.,0.)
max_time = 2
time = 0.
dt = .00005

#initial angular velocity and velocity around loop
omega= -sqrt(2*g*(h-R*(1-cos(theta)))*(1-mu/tan(theta)))/R
v=omega*R

#Initial angular acceleration alfa
alfa = -g*sin(theta)/R + mu*(g*cos(theta)/R + omega*omega)

scale = .07
bv = arrow(pos=car.pos, axis = scale*car.vel,color = color.blue)

#THE LOOP
while time < max_time:
    rate(1500 )
    #make one step
    #new angle
    theta += omega*dt + .5*alfa*dt*dt

    #new cartesian coordinates
    car.pos = vector(R*sin(theta), R*(1-cos(theta)), 0.)

    #new Velocity
    v = omega*R
    car.vel = vector(v*cos(theta),v*sin(theta),0)

    #plot points
    vel_Plot.plot(pos=(time, v))
    pos_Plot.plot(pos=(car.pos.x, car.pos.y))

    #new ball arrow
    bv.pos = car.pos
    bv.axis = scale*car.vel

    #new car trail
    car.trail.append(pos = car.pos)

    #New time
    time += dt

    #new angular acceleration
    omega_old=omega
    if omega < 0:
        sign = -1
    else:
        sign = 1

    alfa = -g*sin(theta)/R - sign*mu*(omega*omega + g*cos(theta)/R)

    #calculate new angular velocity

    omega += alfa*dt

    #Print out the traversal number and maximum height
    #if omega changes sign
    if(omega_old<0. and omega>0.) or (omega_old>0. and omega<0.) :
        k=k+1
        print"pass number =",k
        print "maximum car height for pass", k, "= ", car.pos.y