# Coin sliding off a turntable simulation.
#
# Usage:
# python <script_name> [mu_value [print_option]]
#
# Arguments:
# mu_value - The fractional part of the ratio of kinetic to static friction.
# print_option - 0 to print summary only, 1 to print data for plotting.
#
# Example:
# The following simulates a coin sliding off a turntable with mu=0.8.
# Data suitable for use with gnuplot are saved to the files r8.dat (position),
# v8.dat (velocity vectors), and a8.dat (acceleration vectors).
#
# python <script_name> 8 1
from __future__ import print_function
import copy
import math
import sys # Compute the inertial acceleration vector due to kinetic friction.
# The magnitude of the acceleration is constant but the direction is
# against the relative velocity between the coin and the point on the
# turntable directly beneath the coin.
def accel (x, v, omega, amag) :
rsq = x[0]*x[0] + x[1]*x[1]
rmag = math.sqrt(rsq)
rhat = [x[0]/rmag, x[1]/rmag]
rw = rmag*omega
vp = [-rw*rhat[1], rw*rhat[0]]
vrel = [v[0]-vp[0], v[1]-vp[1]]
vrelsq = vrel[0]*vrel[0] + vrel[1]*vrel[1]
if (vrelsq < 1e-15) :
a = [-amag*rhat[0], -amag*rhat[1]]
else :
vrelmag = math.sqrt(vrelsq)
vrelhat = [vrel[0]/vrelmag, vrel[1]/vrelmag]
a = [-amag*vrelhat[0], -amag*vrelhat[1]]
return a# Advance state by timestep=dt using RK4.
def step (x, v, omega, amag, dt) :
x0 = copy.copy(x)
v0 = copy.copy(v)
a0 = accel (x, v, omega, amag)
(x, v) = step_pos_vel (x0, v0, v0, a0, 0.5*dt)
v1 = copy.copy(v)
a1 = accel (x, v, omega, amag)
(x, v) = step_pos_vel (x0, v1, v0, a1, 0.5*dt)
v2 = copy.copy(v)
a2 = accel (x, v, omega, amag)
(x, v) = step_pos_vel (x0, v2, v0, a2, dt)
a = accel (x, v, omega, amag)
v3 = [(v0[0] + 2.0*(v1[0]+v2[0]) + v[0]) / 6.0,
(v0[1] + 2.0*(v1[1]+v2[1]) + v[1]) / 6.0]
a3 = [(a0[0] + 2.0*(a1[0]+a2[0]) + a[0]) / 6.0,
(a0[1] + 2.0*(a1[1]+a2[1]) + a[1]) / 6.0]
(x, v) = step_pos_vel (x0, v3, v0, a3, dt)
return (x, v)# Advance position and velocity.
def step_pos_vel (x, xdot, v, vdot, dt) :
return ([x[0]+xdot[0]*dt, x[1]+xdot[1]*dt],
[v[0]+vdot[0]*dt, v[1]+vdot[1]*dt])# Integrate state and record data until the coin flies off the turntable.
def integ (x, v, omega, amag, dt, rfinal, mode, smu) :
t = 0
if (mode > 0) :
fd_r = open ("r" + smu + ".dat", 'w')
fd_v = open ("v" + smu + ".dat", 'w')
fd_a = open ("a" + smu + ".dat", 'w')
while (vector_mag(x) < rfinal) :
if (mode > 0) :
print_xva (t, x, v, omega, amag, fd_r, fd_v, fd_a)
(x, v) = step (x, v, omega, amag, dt)
t = t + dt
if (mode > 0) :
print_xva (t, x, v, omega, amag, fd_r, fd_v, fd_a)
print ("t=", t)
print ("x=", x[0], ', ', x[1], ', magnitude = ', vector_mag(x))
print ("v=", v[0], ', ', v[1], ', magnitude = ', vector_mag(v))def print_xva (t, x, v, omega, amag, fd_r, fd_v, fd_a) :
print (x[0], x[1], file=fd_r)
if (abs(int(t*5+0.5) - t*5) < 1e-12) :
a = accel (x, v, omega, amag)
print (x[0], x[1], v[0], v[1], file=fd_v)
print (x[0], x[1], a[0], a[1], file=fd_a)def vector_mag (vec) :
return math.sqrt (vec[0]*vec[0] + vec[1]*vec[1])def main () :
# Get the fractional mu value and output option from the argument list.
# mu is the ratio of kinetic friction to static friction;
# the value must be in the range [0.0, 1.0). This is enforced
# by prefixing "0." to the input value.
if (len(sys.argv) > 1) :
smu = sys.argv[1]
mu = float("0." + smu)
else :
smu = "8"
mu = 0.8
# mode = 0 means print just summary data to stdout.
# A non-zero value prints position, velocity, and acceleration
# data to files in a format suitable for use with gnuplot.
if (len(sys.argv) > 2) :
mode = int(sys.argv[2])
else :
mode = 0
# Turntable angular velocity, in radians per time unit.
omega = 1
# The integration time step.
# This is close to optimal for RK4 for this value of omega.
dt = 0.005
# The initial radius of the coin and and the radius of the turntable.
rinit = 1
rfinal = 4
# The magnitude of the acceleration due to kinetic friction.
amag = mu*rinit*omega*omega
# Initial state; coin just started slipping (so non-slip velocity).
x = [rinit, 0.0]
v = [0.0, rinit*omega]
# Integrate state and record data.
integ (x, v, omega, amag, dt, rfinal, mode, smu)main()
