Description
Python 2.7.5, NumPy 1.7.1, SciPy 0.13.0, Matplotlib 1.3.1
This is a script to simulate the rolling of a coin( a circle) on a surface without slipping in 3 dimensions. I tried to paste on a shorter version of my code but it is not easy because of the physics part. Please bare with me. My question relates to the animation part. If one uncomments the last two lines related to function animation just before the pyplot.show(), we will get a colored circle. But if we keep the function animation, the circle in motion will only have one color. The circle is a 3D line object and, it is being updated by update_line function. It seems the color information is left behind during the update_line or animation period. I tried many methods like, passing in the color argument into the update_line function and use the set_c method, or passing in a colormap and a scale. None of them worked for me. Does anyone know how to give the rolling circle a color gradient/cycle, to make it looks more like a rolling without slipping motion? Thanks in advance! Tony :)
import numpy as np
from matplotlib import pyplot
import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation
def rotation_matrix(axis,theta):
axis = axis/np.sqrt(np.dot(axis,axis))
a = np.cos(theta/2)
b,c,d = axis*np.sin(theta/2)
return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
[2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
[2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])
xx = 5
R = 5*np.sqrt(2)
ncp = 10 # number of 'new' CP, not including origin
onetheta = 30 # the angle of triangles which approximate a circle
angles = range(0,360+onetheta,onetheta) #[0,2,...,360]
theta = np.array(angles) * np.pi / 180 #in radians
C = np.zeros((ncp+1,3,3))
# Set the orientation of z w.r.t lab_z
leananl=np.pi/5
C[0,:,:] = np.dot(rotation_matrix(np.array([1,0,0]),leananl),np.eye(3))
Cprime = np.zeros((ncp+1,3,3))
CM = np.zeros((ncp+1,3)) #CM position vector array
CM[0,:]= R*C[0,:,1] #CM(t0)
CMprime = np.zeros((ncp+1,3))
z = np.array([0,0,1])
cirCP = np.zeros((len(angles),3,ncp)) #circle of the coin
rotheta = onetheta*np.pi/180 # theta
CP = np.zeros((ncp+1,3))
CPprime = np.zeros((ncp+1,3))
lowerangle = -np.radians(np.linspace(0,ncp,ncp)/ncp)
omega_z = np.zeros(ncp+1)
w = np.zeros((ncp+1,3))
w[:,2]=w[:,2]+0.0005
w[:,1]=w[:,1]+0.001
w[:,0]=w[:,0]-0.0005
for j in range(ncp):
betaj=np.pi/2 - np.arccos(np.dot(C[j,:,1],z))
alphaangl = np.arcsin(np.sin(rotheta/2)/np.cos(betaj))
kappaj = 2*alphaangl
phi = np.arccos(np.cos(alphaangl)/np.cos(np.radians(onetheta/2)))
twodt=2*phi/w[j+1,1] #phi/wy
Avecinspace = np.dot(rotation_matrix(z,kappaj/2),
C[j,:,0]*R*np.sin(rotheta/2)*2)
Cprime[j,:,:] = np.dot(rotation_matrix(z,kappaj),C[j,:,:])
C[j+1,:,:] = np.dot(rotation_matrix(Cprime[j,:,0],w[j+1,0]*twodt),Cprime[j,:,:])
CPprime[j,:] = CP[j,:] +Avecinspace
CP[j+1,:] = CPprime[j,:] + Cprime[j,:,0]*w[j+1,2]*twodt*R
CMprime[j,:]=CM[j,:] + Cprime[j,:,0]*w[j+1,2]*twodt*R
CM[j+1,:] = CP[j+1,:] + C[j+1,:,1]*R
for i in range(len(angles)):
cirCP[i,:,j] = np.dot(rotation_matrix(C[j,:,2],theta[i]),CP[j,:]-CM[j,:]) + CM[j,:]
omega_z[j+1] = omega_z[j] + w[j+1,2]*twodt
### plot figure
fig1 = pyplot.figure(1)
ax = p3.Axes3D(fig1)
ax.view_init(elev=35, azim=-45)
def update_line(x, line, cirCP, CP):
for indl, linein in enumerate(line):
for li, angx1 in zip(linein,range(len(angles))):
li.set_data(cirCP[0:2+indl,0,x],cirCP[0:2+indl,1,x])
li.set_3d_properties(cirCP[0:2+indl,2,x])
return line
line = [ax.plot(*[cirCP[j:j+2,i,0] for i in range(3)]) for j in range(len(angles))]
ax.set_xbound(-10,10)
ax.set_ybound(0,20)
ax.set_zbound(0,10)
line_ani = animation.FuncAnimation(fig1, update_line, ncp, fargs=(line,cirCP,CP),
interval=500, blit=False,repeat=True)
pyplot.show()