Sage (terminal velocity of different objects)

sage_terminal_velocity.py

## Credits: 
##
## Collaborative effort:
## - Michel Dusseault wrote the code. 
## - I worked out some of the differential equations.
##

var('t')
g = 9.81  # m/s^2
rho = 1.2   # Rho: Density of air at sea level @ 25C
r = 0.11  # Radius
A = 2*pi*r^2  # Half the surface area of a sphere

labels = ['Cannon ball', 'Soccer ball', 'Car', 'Car (Edge)', 'Pillow', 'Anvil', 'Piano', 'Coyote']
colors = ['red', 'blue', 'orange', 'purple', 'cyan', 'magenta', 'black', 'brown']
C_ds = [0.43, 0.43, 1.28, 0.43, 1.28, 0.8, 1.28, 1.2]
Ms = [43.9, 0.43, 2300.0, 2300.0, 0.5, 90.0, 245.0, 21.0]
As = [2*pi*r^2, 2*pi*r^2, 3.0*1.5, 1.5*1.3, 0.66*0.55, 0.5*0.3, 1.55*1.4, 0.35*0.2]

k_func(C_d, A, p, m) = (C_d*A*p)/(2*m)
Ks = [k_func(C_ds[i], As[i], rho, Ms[i]) for i in xrange(len(labels))]

# True as long as g and k are positive.
v(t, k) = sqrt(g/k)*tanh(sqrt(k*g)*t)
x(t, k) = log(cosh(sqrt(g*k)*t))/k

# Free fall
v_ff(t) = g*t
x_ff(t) = (g*t^2)/2
themax = 30
p = plot(x_ff(t), 0, themax, color='green', label='Free fall')
q = plot(v_ff(t), 0, themax, color='green', label='Free fall')
for i in xrange(len(labels)):
    p += plot(x(t, k=Ks[i]), 0, themax, axes_labels=['Time (sec)', 'Distance fallen (m)'], label=labels[i], rgbcolor=colors[i])
    q += plot(v(t, k=Ks[i]), 0, themax, axes_labels=['Time (sec)', 'Velocity (m/s)'], label=labels[i], rgbcolor=colors[i])

p.show()
q.show()