Bouncing ball simulation in VPython, using three different types of dissipative force to get the decay.

bouncing.py

from visual import *
from math import *

g=9.80
size=0.20
mass=0.2
drop_height = 2.0
k_ball = 1000
drag=0.40
drag2=0.32
fraction = 0.0425

floor = box(pos=(0,0,0), length=5, width=1, height=0.05, color=color.cyan)
ball = sphere(pos=(-floor.length/2, drop_height,0), radius=size, color=color.green)
ball.p=mass*vector(0.2,0,0)
ball.m=mass
ball.trail=curve(color=ball.color)

F_grav=ball.m*vector(0,-g,0)

t=0
tmax=50
dt=0.0001
nbounce=0
go=1

oldvel=ball.p.y/ball.m
oldheight=ball.pos.y
lastbounce=ball.pos.y


while go>0:
    #rate(1000)
    if ball.p.y<0:
        vhat = vector(0,-1,0)
    else:
        vhat = vector(0,+1,0)
        
    if ball.pos.y<ball.radius:
        stretch=ball.radius-ball.pos.y
        F_spring = k_ball*stretch*vector(0,1,0)
        #F_drag = -drag*(ball.p.mag/(ball.m*1))*vhat
        #F_drag = -drag2*(ball.p.mag2/(ball.m*ball.m))*vhat
        F_drag = -fraction*F_spring.mag*vhat
    else:
        F_spring = vector(0,0,0)
        F_drag=vector(0,0,0)
    #F_drag.x=0
    #F_drag.z=0
    F_total = F_grav + F_spring + F_drag
    ball.p = ball.p + F_total*dt
    ball.pos=ball.pos+dt*ball.p/ball.m
    ball.trail.append(pos=ball.pos)
    if (ball.p.y<0) and (oldvel>0):
        #print("Bounce: ", nbounce, "  Height: ", ball.pos.y)
        lastbounce=ball.pos.y
    if (ball.p.y>0) and (oldvel<0):
        nbounce=nbounce+1
        print("  Bounce: ",nbounce, "Time: ", t, "Height: ", lastbounce)
    oldvel=ball.p.y/ball.m
    oldheight=ball.pos.y
    if lastbounce<1.05*ball.radius:
        go=0
    

    t = t+dt
    
print(ball.pos.x+floor.length/2)