# Reaction-Diffusion Simulation Using Gray-Scott Model
# https://en.wikipedia.org/wiki/Reaction-diffusion_system
# FB - 20151016
import math
import random
from PIL import Image, ImageDraw
imgx = 256; imgy = 256 # image size
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)
pixels = image.load()
steps = 1000
dA = 1.0 # Diffusion Rate A
dB = 0.5 # Diffusion Rate B
f = 0.0545
k = 0.062
h = 0.01 # time step size

# 3x3 Convolution Matrix
weights = [[0.05, 0.2, 0.05], [0.2, -1.0, 0.2], [0.05, 0.2, 0.05]]

##arA = [[[1.0 for x in range(imgx)] for y in range(imgy)] for z in range(2)]
##arB = [[[0.0 for x in range(imgx)] for y in range(imgy)] for z in range(2)]
##arB[0][(imgy - 1) / 2][(imgx - 1) / 2] = 1.0 # seed
arA = [[[random.random() for x in range(imgx)] for y in range(imgy)] for z in range(2)]
arB = [[[random.random() for x in range(imgx)] for y in range(imgy)] for z in range(2)]

# simulation
p = 0; print "%" + str(p).zfill(2)
z = 1
for i in range(steps):
    z = 1 - z
    for iy in range(imgy):
        for ix in range(imgx):
            lapA = 0.0
            lapB = 0.0
            for jy in range(3):
                for jx in range(3):
                    kx = ix + (jx - 1)
                    ky = iy + (jy - 1)
                    if kx < 0: kx = imgx - 1
                    if kx > imgx - 1: kx = 0
                    if ky < 0: ky = imgy - 1
                    if ky > imgy - 1: ky = 0
                    lapA += arA[z][ky][kx] * weights[jy][jx]
                    lapB += arB[z][ky][kx] * weights[jy][jx]
            a = arA[z][iy][ix]
            b = arB[z][iy][ix]
            abb = a * b * b
            an = a + (dA * lapA - abb + f * (1.0 - a)) * h
            bn = b + (dB * lapB + abb - (k + f) * b) * h
            # prevent negative chemical amounts
            if an < 0: an = 0.0
            if bn < 0: bn = 0.0
            arA[1 - z][iy][ix] = an
            arB[1 - z][iy][ix] = bn

    pn = 100 * (i + 1) / steps # percent completed
    if pn != p:
        p = pn
        print "%" + str(p).zfill(2)

# paint the final state
z = 1 - z
aMin = arA[z][0][0]; aMax = aMin
bMin = arB[z][0][0]; bMax = bMin
for iy in range(imgy):
    for ix in range(imgx):
        a = arA[z][iy][ix]
        b = arB[z][iy][ix]
        if a < aMin: aMin = a
        if a > aMax: aMax = a
        if b < bMin: bMin = b
        if b > bMax: bMax = b

if aMin != aMax:
    for iy in range(imgy):
        for ix in range(imgx):
            a = arA[z][iy][ix]
            cA = int(255 * (a - aMin) / (aMax - aMin))
            pixels[ix, iy] = (cA, cA, cA)
    # label = "f = " + str(f) + " k = " + str(k)
    # draw.text((0, 0), label, (0, 255, 0))
    image.save("ReactionDiffusionSim_A.png", "PNG")

if bMin != bMax:
    for iy in range(imgy):
        for ix in range(imgx):
            b = arB[z][iy][ix]
            cB = int(255 * (b - bMin) / (bMax - bMin))
            pixels[ix, iy] = (cB, cB, cB)
    # label = "f = " + str(f) + " k = " + str(k)
    # draw.text((0, 0), label, (0, 255, 0))
    image.save("ReactionDiffusionSim_B.png", "PNG")