Simulate the Belousov Zhabotinsky

belousov.c

#include "raylib.h"
#include <math.h>

#define SCREEN_WIDTH 500
#define SCREEN_HEIGHT 500
#define CELL_SIZE 4    
#define GRID_WIDTH (SCREEN_WIDTH / CELL_SIZE)
#define GRID_HEIGHT (SCREEN_HEIGHT / CELL_SIZE)

typedef struct {
	float a;
	float b;
} Cell;

Cell grid[GRID_WIDTH][GRID_HEIGHT];
Cell next[GRID_WIDTH][GRID_HEIGHT];


const float Da = 0.082f;
const float Db = 0.041f;
const float f = 0.014f;
const float k = 0.040f;
const float dt = 0.05f;


int spawnTimer = 0;
const int spawnInterval = 120; 

float laplace(float v[GRID_WIDTH][GRID_HEIGHT], int x, int y) {
	float sum = 0;
	int dx[] = {0, 1, 0, -1, 1, 1, -1, -1};  
	int dy[] = {-1, 0, 1, 0, -1, 1, 1, -1};

	for(int i = 0; i < 8; i++) {
		int nx = (x + dx[i] + GRID_WIDTH) % GRID_WIDTH;
		int ny = (y + dy[i] + GRID_HEIGHT) % GRID_HEIGHT;
		sum += v[nx][ny];
	}
	return sum - 8 * v[x][y];  
}

void updateGrid(void) {
	float a[GRID_WIDTH][GRID_HEIGHT];
	float b[GRID_WIDTH][GRID_HEIGHT];

	for(int x = 0; x < GRID_WIDTH; x++) {
		for(int y = 0; y < GRID_HEIGHT; y++) {
			a[x][y] = grid[x][y].a;
			b[x][y] = grid[x][y].b;
		}
	}

	for(int x = 0; x < GRID_WIDTH; x++) {
		for(int y = 0; y < GRID_HEIGHT; y++) {
			float A = grid[x][y].a;
			float B = grid[x][y].b;

			float dA = Da * laplace(a, x, y) + A * (1.0f - A) - (k + f) * A - B * A * A;
			float dB = Db * laplace(b, x, y) + (A * A * B - B * f) * 0.5f;

			next[x][y].a = A + dt * dA;
			next[x][y].b = B + dt * dB;

			next[x][y].a = fmaxf(0.0f, fminf(1.0f, next[x][y].a));
			next[x][y].b = fmaxf(0.0f, fminf(1.0f, next[x][y].b));
		}
	}

	for(int x = 0; x < GRID_WIDTH; x++) {
		for(int y = 0; y < GRID_HEIGHT; y++) {
			grid[x][y] = next[x][y];
		}
	}
}

void addNewSeeds(void) {
	int numSeeds = GetRandomValue(2, 4);
	for(int i = 0; i < numSeeds; i++) {
		int px = GetRandomValue(0, GRID_WIDTH-1);
		int py = GetRandomValue(0, GRID_HEIGHT-1);

		float intensity = GetRandomValue(70, 100) / 100.0f;

		for(int x = px-2; x < px+2; x++) {
			for(int y = py-2; y < py+2; y++) {
				int xx = (x + GRID_WIDTH) % GRID_WIDTH;
				int yy = (y + GRID_HEIGHT) % GRID_HEIGHT;
				grid[xx][yy].a = intensity;
				grid[xx][yy].b = 1.0f - intensity;
			}
		}
	}
}

Color getVibrantColor(float a, float b) {

	float value = a / (a + b + 0.1f);


	if (value < 0.2f) {
		return (Color){0, 0, 255, 255};          
	} else if (value < 0.4f) {
		return (Color){0, 255, 255, 255};        
	} else if (value < 0.6f) {
		return (Color){255, 0, 255, 255};        
	} else if (value < 0.8f) {
		return (Color){255, 255, 0, 255};        
	} else {
		return (Color){255, 0, 0, 255};          
	}
}

int main(void) {
	InitWindow(SCREEN_WIDTH, SCREEN_HEIGHT, "BZ Reaction Simulation");
	SetTargetFPS(60);


	for(int x = 0; x < GRID_WIDTH; x++) {
		for(int y = 0; y < GRID_HEIGHT; y++) {
			grid[x][y].a = GetRandomValue(0, 10) / 100.0f;
			grid[x][y].b = GetRandomValue(0, 10) / 100.0f;
		}
	}


	addNewSeeds();

	while (!WindowShouldClose()) {

		spawnTimer++;
		if(spawnTimer >= spawnInterval) {
			addNewSeeds();
			spawnTimer = 0;
		}

		for(int i = 0; i < 2; i++) {
			updateGrid();
		}

		BeginDrawing();
		ClearBackground(BLACK);

		for(int x = 0; x < GRID_WIDTH; x++) {
			for(int y = 0; y < GRID_HEIGHT; y++) {
				float a = grid[x][y].a;
				float b = grid[x][y].b;

				Color color = getVibrantColor(a, b);


				Vector2 center = {
					x * CELL_SIZE + CELL_SIZE/2.0f,
					y * CELL_SIZE + CELL_SIZE/2.0f
				};
				DrawCircleV(center, CELL_SIZE/2.0f, color);
			}
		}

		DrawFPS(10, 10);
		EndDrawing();
	}

	CloseWindow();
	return 0;
}