heat.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

#define SDL_MAIN_HANDLED
#include <SDL2/SDL.h>

void* xmalloc(size_t size) {
	void* ptr = malloc(size);
	if (ptr == NULL) {
		fprintf(stderr, "Failed to allocate %lld bytes\n", size);
		exit(EXIT_FAILURE);
	}
	return ptr;
}

void memset_float(float* dst, float val, size_t n) {
	for (size_t i = 0; i < n; i++) {
		dst[i] = val;
	}
}

int clampi(int x, int min, int max) {
	if (x < min) {
		return min;
	}
	if (x > max) {
		return max;
	}
	return x;
}

int width = 640;
int height = 480;

SDL_Window* window;
SDL_Renderer* renderer;

bool running = true;
bool paused = false;

int num_points;
float dx;
float dt;
float alpha;
float r;
float* heat_current;
float* heat_old;

float* a;
float* b;
float* c;
float* d;

void init_sim() {
	dx = 1.0f / num_points;
	r = alpha * dt / (dx * dx);

	heat_current = xmalloc(num_points * sizeof(float));
	heat_old = xmalloc(num_points * sizeof(float));

	for (int i = 0; i < num_points; i++) {
		float x = (float)i / (num_points - 1);
		float a = 6 * (x - 0.5);
		heat_old[i] = expf(-a * a);
	}

	a = xmalloc((num_points - 2) * sizeof(float));
	b = xmalloc((num_points - 3) * sizeof(float));
	c = xmalloc((num_points - 3) * sizeof(float));
	d = xmalloc((num_points - 2) * sizeof(float));
}

void solve_tridiagonal(int n, float* a, float* b, float* c, float* d) {
	for (int i = 0; i < n - 1; i++) {
		float k = -c[i] / a[i];
		a[i + 1] += k * b[i];
		d[i + 1] += k * d[i];
	}

	d[n - 1] /= a[n - 1];
	for (int i = n - 2; i >= 0; i--) {
		d[i] = (d[i] - b[i] * d[i + 1]) / a[i];
	}
}

void step() {
	heat_current[0] = 0;
	heat_current[num_points - 1] = 0;

	memset_float(a, 1 + 2 * r, num_points - 2);
	memset_float(b, -r, num_points - 3);
	memset_float(c, -r, num_points - 3);
	for (int i = 0; i < num_points - 2; i++) {
		d[i] = heat_old[i + 1];
	}
	d[0] += r * heat_old[0];
	d[num_points - 3] += r * heat_old[num_points - 1];

	solve_tridiagonal(num_points - 2, a, b, c, d);

	memcpy(heat_current + 1, d, (num_points - 2) * sizeof(float));

	float* tmp = heat_current;
	heat_current = heat_old;
	heat_old = tmp;
}

int main() {
	SDL_Init(SDL_INIT_EVERYTHING);

	window = SDL_CreateWindow("Heat", SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED, width, height, 0);
	renderer = SDL_CreateRenderer(window, -1, SDL_RENDERER_PRESENTVSYNC);

	num_points = width;
	dt = 0.01;
	alpha = 0.01f;
	init_sim();

	while (running) {
		SDL_Event e;
		while (SDL_PollEvent(&e)) {
			switch (e.type) {
			case SDL_QUIT:
				running = false;
				break;
			case SDL_KEYDOWN:
				if (e.key.keysym.sym == SDLK_SPACE) {
					paused = !paused;
				}
			default:
				break;
			}
		}

		int edit_radius = 15;

		int mousex;
		int mousey;
		Uint32 buttons = SDL_GetMouseState(&mousex, &mousey);
		if (buttons & SDL_BUTTON(SDL_BUTTON_LEFT)) {
			float val = 1 - 2.0f / height * mousey;

			int leftx = clampi(mousex - edit_radius, 0, num_points);
			int rightx = clampi(mousex + edit_radius, 0, num_points);

			for (int x = leftx; x < rightx; x++) {
				heat_old[x] = val;
				heat_current[x] = val;
			}
		}

		if (!paused) {
			step();
		}

		SDL_SetRenderDrawColor(renderer, 0, 0, 0, 255);
		SDL_RenderClear(renderer);

		SDL_SetRenderDrawColor(renderer, 255, 255, 255, 255);
		for (int i = 0; i < num_points; i++) {
			float val = heat_current[i];
			float y = (1 - val) / 2 * height;
			SDL_RenderDrawPoint(renderer, i, y);
		}

		SDL_RenderPresent(renderer);
	}

	return 0;
}