3次元版の波。自由端反射のpythonプログラム смотреть онлайн

これは立方体の中を伝わる波です。右のグラフを倒したものが左のグラフです。3次元にすることは簡単だが、どこがどこへ伝達したかが分かるように表現するのが難しい。
試行錯誤の末、この形にしてみたが、これでも分からないと思います。。
これは慣性と伝達の単純処理を繰り返す、自由端反射のpythonプログラムです。
仕組みは1次元版と2次元版に記載しています。

ソースリストは以下の通り。

--- ここから ---
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

X_SU = 140
Y_SU = 40
H_SU = 44 #3+8*BLK # 11+8xn 01.2.3.4.5.6.7.8.90

fig = plt.figure(figsize = (12.8, 7.2))
ax = fig.add_subplot(121, projection="3d")
ax2 = fig.add_subplot(122, projection="3d")

plt.subplots_adjust(left=0, bottom=0, right=0.99, top=1, wspace=0.01)

x, y = np.meshgrid(range(0,X_SU),range(0,Y_SU))
hx, hxh = np.meshgrid(range(0,X_SU),range(0,H_SU))
hy, hyh = np.meshgrid(range(0,Y_SU),range(0,H_SU))
z = np.full((H_SU,Y_SU,X_SU),0.0)

z1 = z.copy()
k = z.copy()
z[2,-2,2] = 20.0 * Y_SU*2 # 左下端

def update(i):
global z,z1,k
print(i)
z[:,0,:] = z[:,1,:]
z[:,-1,:] = z[:,-2,:]
z[:,:,0] = z[:,:,1]
z[:,:,-1] = z[:,:,-2]
z[0,:,:] = z[1,:,:]
z[-1,:,:] = z[-2,:,:]
for ih in range(1,H_SU-1):
for ix in range(1,X_SU-1):
for iy in range(1,Y_SU-1):
z1[ih,iy,ix] = (np.sum(z[ih,iy-1:iy+2,ix]) + z[ih,iy,ix-1] + z[ih,iy,ix+1] + z[ih-1,iy,ix] + z[ih+1,iy,ix])/7 + k[ih,iy,ix]
if i GT 6: # GTは不等号に置き換えること (youtubeは不等号を許可していないため)
k = z1 - z
ax.clear()
vminmax = np.max(z)

ax.plot_surface(x, y, z[1], cmap="viridis", vmin=-vminmax, vmax=vminmax)

X5 = int(X_SU/5)
H2 = int(H_SU/2)+2
ax.plot_wireframe(z[::-1,::-1,0 ], hyh, hy+H2, linewidth=0.5)
ax.plot_wireframe(z[::-1,::-1,X5 ]+X5, hyh, hy+H2, linewidth=0.5)
ax.plot_wireframe(z[::-1,::-1,X5*2]+X5*2, hyh, hy+H2, linewidth=0.5)
ax.plot_wireframe(z[::-1,::-1,X5*3]+X5*3, hyh, hy+H2, linewidth=0.5)
ax.plot_wireframe(z[::-1,::-1,X5*4]+X5*4, hyh, hy+H2, linewidth=0.5)
ax.plot_wireframe(z[::-1,::-1,-1]+X5*5-1, hyh, hy+H2, linewidth=0.5)

H2 = H_SU*2
ax.plot_surface(x, y, z[-2]+H2, cmap="viridis", vmin=-vminmax+H2, vmax=vminmax+H2)

X5 = int(X_SU/10)
ax.set_zlim(0,H2)
ax.set_xlim(X5,X_SU-X5)
ax.set_ylim(Y_SU*2,-Y_SU)

ax2.clear()

ax2.plot_wireframe(z[1]-int(H_SU*0.5), y, x, linewidth=0.5)

X5 = int(X_SU/5)
ax2.plot_surface(hyh, hy, z[:,:,0 ], cmap="viridis", vmin=-vminmax, vmax=vminmax)
ax2.plot_surface(hyh, hy, z[:,:,X5 ]+X5, cmap="viridis", vmin=-vminmax+X5, vmax=vminmax+X5)
ax2.plot_surface(hyh, hy, z[:,:,X5*2]+X5*2, cmap="viridis", vmin=-vminmax+X5*2, vmax=vminmax+X5*2)
ax2.plot_surface(hyh, hy, z[:,:,X5*3]+X5*3, cmap="viridis", vmin=-vminmax+X5*3, vmax=vminmax+X5*3)
ax2.plot_surface(hyh, hy, z[:,:,X5*4]+X5*4, cmap="viridis", vmin=-vminmax+X5*4, vmax=vminmax+X5*4)
ax2.plot_surface(hyh, hy, z[:,:,-1]+X5*5-1, cmap="viridis", vmin=-vminmax+X5*5-1, vmax=vminmax+X5*5-1)

ax2.plot_wireframe(z[-2]+int(H_SU*1.5), y, x, linewidth=0.5)

H2 = int(H_SU*0.5)
ax2.set_xlim(H2+H_SU,-H2)
ax2.set_zlim(X5,X_SU-X5)
ax2.set_ylim(Y_SU*2,-Y_SU)

z = z1.copy()

ani = animation.FuncAnimation(fig, update, interval=120, blit=False, save_count=400)
#ani.save("nami3d_fe_2scr.mp4")
plt.show()
--- ここまで ---