Magimenal Mirai: ベロウソフ・ジャボチンスキー反応のモデルで遊びたかった

　※注意
　　この記事には強い点滅を伴うgifが貼ってあります。

　毎度ながら「遊ぶ」まで至らずに終わりました。ただしつい先日学んだ関数内関数を使えたというだけで個人的には満足なのです。
　なお関数内関数を使う必要があったかどうかについては考えてはいけない。
　これを参考にした。
---------------------------------------------------------------------------------
#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def defs(k1, k2, g, n, nei):
    def BZreaction(plate):
        padded=zeropad(plate)
        output=np.zeros(padded.shape)
        i=nei
        while i<padded.shape[0]-nei:
            j=nei
            while j<padded.shape[0]-nei:
                if padded[i, j]==0:
                    b=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>=n)
                    a=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>0)-b
                    output[i, j]=a/k1+b/k2
                elif padded[i, j]>=n:
                    output[i, j]=0
                else:
                    s=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1])
                    b=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>=n)
                    a=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>0)-b
                    output[i, j]=s/a+b+1+g
                j+=1
            i+=1

        return output[nei:-1*nei, nei:-1*nei]

    def zeropad(plate):
        padded=np.zeros((plate.shape[0]+2*nei, plate.shape[0]+2*nei))
        padded[nei:plate.shape[0]+nei, nei:plate.shape[1]+nei]=plate
        return padded

    return BZreaction, zeropad

#setting
g=10
k1=2
k2=1
n=20
l=50
plate=np.random.rand(l*l).reshape(l, l)*n
nei=1

BZreaction, zeropad=defs(k1, k2, g, n, nei)

fig=plt.figure()

ims=[]
for a in range(500):
plate=BZreaction(plate)
im=plt.imshow(plate, animated=True)
ims.append([im])

ani=animation.ArtistAnimation(fig, ims, interval=50, blit=True, repeat_delay=1000)
ani.save('sample.gif', writer='imagemagick')
---------------------------------------------------------------------------------