Reaction-diffusion model on cellular automata
Alasdair Turner (2009) A Simple Model of the Belousov-Zhabotinsky Reaction from First Principles
implementation note: discovery.ucl.ac.uk/id/eprint/17241/1/17241.pdf
en.wikipedia.org/wik…zhabotinsky_reaction
Log in to post a comment.
const W = 300; const N1 = 0; const N2 = 10; const N3 = 20; const N4 = 40; const N5 = 60; const N6 = 80; const N7 = 100; const N8 = 200; const N9 = 400; const N = [N1, N2, N3, N4, N5, N6, N7, N8, N9]; let done = 0; let D = 60; Canvas.setpenopacity(-0.75); const t = new Turtle(); for (let x = 0; x < 3; x++) { for (let y = 0; y < 3; y++) { let cx = -100 + 200 / 3 * (x + 0.5); let cy = -100 + 200 / 3 * (y + 0.5); for (let r = D/2; r <= D/2 + 1; r += 0.1) { t.jump(cx - r, cy); t.seth(-90); t.circle(r); } } } let a = [], a_buf = []; let b = [], b_buf = []; let c = [], c_buf = []; for (let x = 0; x < W; x++) { a[x] = []; a_buf[x] = []; b[x] = []; b_buf[x] = []; c[x] = []; c_buf[x] = []; for (let y = 0; y < W; y++) { a[x][y] = Math.random(); a_buf[x][y] = 0.0; b[x][y] = Math.random(); b_buf[x][y] = 0.0; c[x][y] = Math.random(); c_buf[x][y] = 0.0; } } function walk(frame) { if (frame == N[done]) { let cx = -100 + 200 / 3 * (done % 3 + 0.5); let cy = -100 + 200 / 3 * (Math.floor(done / 3) + 0.5);; for (let x = 0; x < W; x++) { for (let y = 0; y < W; y++) { if ((x - W/2)**2 + (y - W/2)**2 < (W/2)**2) { let ac = a[x][y]; t.jump(cx + (x - W/2) * D / W, cy + (y - W/2) * D / W); t.seth(360 * Math.random()); t.forward(ac * 0.2); } } } done++; } for (let x = 0; x < W; x++) { for (let y = 0; y < W; y++) { let an = 0.0, bn = 0.0, cn = 0.0; for (let dx = -1; dx <= +1; dx++) { for (let dy = -1; dy <= +1; dy++) { an += a[(x + dx + W) % W][(y + dy + W) % W]; bn += b[(x + dx + W) % W][(y + dy + W) % W]; cn += c[(x + dx + W) % W][(y + dy + W) % W]; } } an /= 9; bn /= 9; cn /= 9; a_buf[x][y] = an + an * (bn - cn); b_buf[x][y] = bn + bn * (cn - an); c_buf[x][y] = cn + cn * (an - bn); } } for (let x = 0; x < W; x++) { for (let y = 0; y < W; y++) { a[x][y] = Math.min(Math.max(a_buf[x][y], 0.0), 1.0); b[x][y] = Math.min(Math.max(b_buf[x][y], 0.0), 1.0); c[x][y] = Math.min(Math.max(c_buf[x][y], 0.0), 1.0); } } return done < N.length; }