I used Python along with some common libraries to run a Belousov-Zhabotinsky simulation I found online. After playing around with it, I decided I wanted a symmetric pattern (actually to make a drink coaster) and my favorite symmetry is the kind you get packing rhombuses with 120° rotation. I already worked on embedding these in a 2D array for hex CAs, so it wasn't hard to apply to this case.
Code: Select all
import math
import numpy as np
from scipy.ndimage import convolve
from sys import argv
from PIL import Image
# Simulation method modified from https://scipython.com/blog/simulating-the-belousov-zhabotinsky-reaction/
# Side length of base rhombus (stored as an nXn grid)
n = 200 if len(argv) < 2 else int(argv[1])
# Value for all Belousov-Zhabotinsky reaction coefficients.
param = 1 if len(argv) < 3 else float(argv[2])
# Number of iterations to run.
niter = 500 if len(argv) < 4 else int(argv[3])
# File to write image to.
imgfile = "./bzrhombus.png" if len(argv) < 5 else argv[4]
# Reaction parameters.
alpha, beta, gamma = param, param, param
# Shape of hexagonal neighborhood out to radius 3.
FILTER_STR = '''
1111000
1111100
1111110
1111111
0111111
0011111
0001111
'''
FILTER_INT = [[int(c) for c in row] for row in FILTER_STR.split()]
FILTER = np.array([FILTER_INT])/sum(sum(row) for row in FILTER_INT)
def update(p,arr):
"""Update arr[p] to arr[q] by evolving in time."""
# Count the average amount of each species in the neighborhood around each cell
# by convolution with FILTER
s = convolve(arr[p], FILTER, mode='wrap')
# Apply the reaction equations
q = (p+1) % 2
arr[q,0] = s[0] + s[0]*(alpha*s[1] - gamma*s[2])
arr[q,1] = s[1] + s[1]*(beta*s[2] - alpha*s[0])
arr[q,2] = s[2] + s[2]*(gamma*s[0] - beta*s[1])
# Ensure the species concentrations are kept within [0,1].
np.clip(arr[q], 0, 1, arr[q])
return arr
def equivalence(i, j, n):
"""Find the equivalent coordinates in the base rhombus by rotating through symmetries."""
rot = 0
while True:
if ((i // n) + (j //n)) % 3 == 0:
return i % n, j % n
i, j = -j - 1, i - j
rot += 1
# Initialize the base rhombus with random amounts of A, B and C.
data = np.random.random(size=(2, 3, n, n))
# Build a 3nX3n patch based on 120-degree rhombic packing.
# This can be tiled with simple toroidal wrapping at the boundaries.
arr = np.zeros((2, 3, n * 3, n * 3))
for i in range(arr.shape[2]):
for j in range(arr.shape[3]):
ti, tj = equivalence(i, j, n)
arr[:, :, i, j] = data[:, :, ti, tj]
# run the iteration
for i in range(0, niter):
arr = update(i % 2, arr)
# Create an image from just the top nXn square, equivalent to the base rhombus
disp = np.concatenate((arr[(i + 1) % 2].transpose(2, 1, 0)[:n, :n, :], np.full((n, n, 1), 1, dtype=np.uint8)), axis=2)
img = Image.fromarray((disp * 255).astype(np.uint8), mode='RGBA')
w, h = img.size
img = img.transform((int(w * 3), int(h * math.sqrt(3))),
Image.Transform.AFFINE, (0.5, 0.5 / math.sqrt(3), -0.5 * w, 0, 1 / math.sqrt(3), 0),
resample=Image.Resampling.BICUBIC)
img.save(imgfile)
Since the built-in convolution does not have the rotational symmetry I need, I duplicated the initial nXn data into a 3nX3n array, at which point, the patch can be tiled without any rotation. This means doing about 9x as much computation as needed. I prefer this to adding code to pad around boundaries each pass, though that would also work.
Here's a rhombus I produced with the command "python bz_rhombus.py 200 1 500 ./bzrhombus.png"
- bzrhombus.png (212.89 KiB) Viewed 654 times
- bzhex.png (581.94 KiB) Viewed 654 times
- bzhexemboss.png (795.29 KiB) Viewed 654 times