2016里约奥运圣火模拟
奥运圣火
今年的奥运圣火非常有创意, 没有做复杂的点火装置, 全靠后部的圣火盆衬托. 而这个装置的精妙设计以及恰到好处的反光却带来了惊艳的视觉效果:
很容易看出它的几何结构, 中间的圆不变, 30 个十字架型以相同的角速度旋转, 但是初相并不一样, 由此带来炫目的螺旋效果.
公式很容易写出(当然我只是目测得到的公式, 不一定对): \[ x = r\times \sin\left(\frac{\theta}{4} + \frac{\pi i}{2} + \omega t\right)\] \[ y = r\times \cos\left(\theta\right) + r\times\cos\left(\frac{\theta}{4} + \frac{\pi i}{2} + \omega t\right)\cos\left(\theta\right)\] \[ z = r\times \sin\left(\theta\right) + r\times\cos\left(\frac{\theta}{4} + \frac{\pi i}{2} + \omega t\right)\sin\left(\theta\right)\]
由于对这个几何体特别有感觉, 我打算用 python 模拟一下.
实现
之前用 matplotlib 这个库也只是做物理实验报告之类的, 并没有写过动画效果, 在看了这位仁兄的Lorenz System之后, 感觉用 matplotlib.animation.FuncAnimation 这个函数还是不难写的, 只需要写一个 animate 每帧更新一下就可以了.
更新时有些细节需要注意, 对于所有 ax.plot, 更新它的坐标需要两步:
line.set_data(x, y) line.set_3d_properties(z)
这是因为 set_data 只能接收两个参数.
对于所有 ax.scatter, 更新它的坐标只需要一步:
pts._offsets3d = (x, y, z)
因为 offsets 只接收两个参数, 而且没有类似于 set_3d_properties 这样的函数.
做好的效果是这样:
由于没有处理前后关系，所以会有一种缺少层次感的感觉(暂时也没想到什么好办法解决).
导出 gif
可以参考这里.
代码
"""sketch map disk disk | | sail: o - o - o - o o | o -- sail | windmill: o - o - o - o - o | \ o sail | o .... . .--windmill ring: . . . .--windmill .... """ import matplotlib as mpl from mpl_toolkits.mplot3d import Axes3D import numpy as np import matplotlib.pyplot as plt from matplotlib import animation # parameters num_frames = 150 # number of animation frames num_disks = 5 # number of disks per sail num_sails = 4 # number of sails per windmill num_windmills = 30 # number of windmills r = 10 # radius # init matplotlib plt.clf() fig = plt.figure() ax = fig.gca(projection='3d') ax.set_xlim3d(-20, 20) ax.set_ylim3d(-20, 20) ax.set_zlim3d(-15, 15) ax.axis('off') # init ring, disks and sails # ring: A single plot featuring the ring ring = ax.plot([], [], [], c='y') # disks: disks scattered on the sails, grouped by distance to windmill axis disks = [ax.scatter([], [], [], c='y', s=10*_) for _ in list(range(num_disks-1, 0, -1)) + [num_disks]] # sails: sails of all windmills sails = sum([ax.plot([], [], [], c='y') for _ in range(num_sails*num_windmills)], []) def init(): ring[0].set_data([], []) ring[0].set_3d_properties([]) for sail in sails: sail.set_data([], []) sail.set_3d_properties([]) for disk in disks: disk._offsets3d = ([], [], []) return ring + sails + disks def animate(frame): theta = np.linspace(-np.pi, np.pi, num_windmills) phi = frame * 0.1 # x, y, z: 1d vector w/ shape(num_windmills, ) # 3d coordinate of windmill axis, used as base position to compute coordinate of disks z = r * np.sin(theta) x = np.zeros(num_windmills) y = r * np.cos(theta) # set coordinate of ring ring[0].set_data(x, y) ring[0].set_3d_properties(z) # X, Y, Z stores the 3d coordinates of all disks X = np.zeros((num_disks, num_sails, num_windmills)) Y = np.zeros((num_disks, num_sails, num_windmills)) Z = np.zeros((num_disks, num_sails, num_windmills)) # the following function could be accelerated with pure tensor operations. # however for readability we still use loops. sails_iter = iter(sails) for i in range(num_sails): for j in range(num_disks): ratio = (j + 1.) / num_disks X[j, i] = x + ratio * r * np.sin(theta / 4 + i * np.pi / 2 + phi) Y[j, i] = y + ratio * r * \ np.cos(theta / 4 + i * np.pi / 2 + phi) * np.cos(theta) Z[j, i] = z + ratio * r * \ np.cos(theta / 4 + i * np.pi / 2 + phi) * np.sin(theta) # set coordinate of sails for k in range(num_windmills): sail = next(sails_iter) sail.set_data([x[k], X[-1, i, k]], [y[k], Y[-1, i, k]]) sail.set_3d_properties([z[k], Z[-1, i, k]]) # set coordinate of disks for j, disk in enumerate(disks): disk._offsets3d = (X[j].flatten(), Y[j].flatten(), Z[j].flatten()) ax.view_init(0, frame) fig.canvas.draw() return ring + sails + disks anim = animation.FuncAnimation(fig, animate, init_func=init, frames=num_frames, interval=30, blit=True) anim.save('flame.gif', fps=30, writer='imagemagick') # uncomment the following line for prompt result # plt.show()