波の回折と干渉。Blenderのpythonプログラム

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

--- ここから ---
# memo
# 1. select plane
# 2. select object data property and delete all shape keys
# 3. run this program

# scratch
# 1. delete a plane mesh
# 2. add a plane mesh
# 3. select object property and change scale x,y to 5.0,10.0
# 4. select plane and change the mode to edit mode
# 5. ctrl+R to devide plane. mouse scroll up to numbers of lines
# 6. add a material

from operator import itemgetter, countOf
import bpy
import numpy as np

TIMESTEP = 3
T_STEP = 20
SLIT_POS = 10
WALL_POS = 10
SRC_POW = 0.4

me = bpy.context.object

a = []

for v in me.data.vertices:
a.append([v.co.y, v.co.x, v.co.z, v.index])

a.sort()
X_SU = countOf(map(itemgetter(0),a), -1.0)
Y_SU = int(len(a)/X_SU)
print('X_SU %d' % X_SU, 'Y_SU %d' % Y_SU)
b = []
for i in range(0,len(a),X_SU):
b.extend(sorted(a[i:i+X_SU], key=itemgetter(1)))
print('i=%d '%i, 'len(b)=%d'%len(b))

aa = np.array(b)

z = np.zeros((Y_SU+2, X_SU+2))
z1 = z.copy()
k = z.copy()

edge = np.full((Y_SU+2,X_SU+2),False)
wall = int(X_SU/WALL_POS)
slit1 = int(Y_SU/2+Y_SU/SLIT_POS+0.5)-1
slit2 = int(Y_SU/2-Y_SU/SLIT_POS+0.5)-1

edge[:,0] = True
edge[:,-1] = True
edge[:,wall] = True
edge[slit1,wall] = False
edge[slit2,wall] = False

for i in range(Y_SU+2):
for j in range(X_SU+2):
if edge[i][j]:
print('o',end='')
else:
print('.',end='')
print(':')

for i in range(0,bpy.context.scene.frame_end,TIMESTEP):

if (i%100) == 0:
print(int(i/TIMESTEP))
else:
print(int(i/TIMESTEP), end=' ')

# set both y-end to absorbed end
z[ 0,:] = z[ 1,:] - k[ 1,:] * 2.32475437
z[-1,:] = z[-2,:] - k[-2,:] * 2.32475437

z[int(Y_SU/2+0.5)-1,:wall] = SRC_POW*np.sin(2*np.pi*i/(TIMESTEP*T_STEP))

for ix in range(1,X_SU+1):
for iy in range(1,Y_SU+1):
if not edge[iy][ix]:
z1[iy,ix] = np.sum(z[iy-1:iy+2,ix]) + z[iy,ix-1] + z[iy,ix+1]
if edge[iy-1,ix]:
z1[iy,ix] += z[iy,ix]
if edge[iy+1,ix]:
z1[iy,ix] += z[iy,ix]
if edge[iy,ix-1]:
z1[iy,ix] += z[iy,ix]
if edge[iy,ix+1]:
z1[iy,ix] += z[iy,ix]
z1[iy,ix] = z1[iy,ix]/5 + k[iy,ix]
k = (z1 - z)

iy = 1
ix = 1

for j in aa:
me.data.vertices[int(j[3])].co.z = z[iy,ix]
ix += 1
if ix GT X_SU: # replace GT to angled blacket (youtube doesn't allow it)
ix = 1
iy += 1
bpy.context.scene.frame_set(i)

bpy.ops.object.shape_key_add(from_mix=False)
bpy.ops.object.active_shape_key_index
me.active_shape_key.value = 1
me.active_shape_key.keyframe_insert(data_path="value", frame=i)
me.active_shape_key.value = 0
me.active_shape_key.keyframe_insert(data_path="value", frame=i-TIMESTEP)
me.active_shape_key.keyframe_insert(data_path="value", frame=i+TIMESTEP)

z = z1.copy()

print('finished.')
--- ここまで ---