2016里约奥运圣火模拟

Table of Contents

奥运圣火

今年的奥运圣火非常有创意, 没有做复杂的点火装置, 全靠后部的圣火盆衬托. 而这个装置的精妙设计以及恰到好处的反光却带来了惊艳的视觉效果:

c7959f5d683a4e34882bc17cedd35849_th.gif

很容易看出它的几何结构, 中间的圆不变, 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 这样的函数.

做好的效果是这样:

flame.gif

由于没有处理前后关系,所以会有一种缺少层次感的感觉(暂时也没想到什么好办法解决).

导出 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()

参考文献

Author: expye(Zihao Ye)

Email: expye@outlook.com

Date: 2016-08-06

Last modified: 2021-10-25 Mon 21:58

Licensed under CC BY-NC 4.0