charged particles in a box

charges.c

#include <SDL.h>
#include <SDL_image.h>
#include <stdio.h>
#include <math.h>

#define WIDTH 200
#define HEIGHT 200

#define MAX_PARS 128

#define TIME_SCALE 300

/*#define VSYNC*/

Uint8 lin2srgb(double i)
{
  return (Uint8)(0.5 + 255.0 * ((i <= 0.0031308) ? (12.92 * i) : ((1.055 * pow(i, 1/2.4)) - 0.055)));
}

void drawpnt(SDL_Renderer *ren, double x, double y, double r, double g, double b)
{
	int cx = (int)x;
	int cy = (int)y;
	double d, wx1, wx2, wx3, wy1, wy2, wy3, w;
	d = x - cx;
	wx1 = 0.5 * (d - 1.0) * (d - 1.0);
	wx2 = 0.75 - (d - 0.5) * (d - 0.5);
	wx3 = 0.5 * d * d;
	d = y - cy;
	wy1 = 0.5 * (d - 1.0) * (d - 1.0);
	wy2 = 0.75 - (d - 0.5) * (d - 0.5);
	wy3 = 0.5 * d * d;

	w = wx1 * wy1;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx - 1, cy - 1);
	w = wx2 * wy1;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx    , cy - 1);
	w = wx3 * wy1;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx + 1, cy - 1);
	w = wx1 * wy2;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx - 1, cy    );
	w = wx2 * wy2;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx    , cy    );
	w = wx3 * wy2;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx + 1, cy    );
	w = wx1 * wy3;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx - 1, cy + 1);
	w = wx2 * wy3;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx    , cy + 1);
	w = wx3 * wy3;
	SDL_SetRenderDrawColor(ren, lin2srgb(r * w), lin2srgb(g * w), lin2srgb(b * w), 255);
	SDL_RenderDrawPoint(ren, cx + 1, cy + 1);
}

struct particle {
	double x, y;
	double m, q;
	double fx, fy;
	double vx, vy;	
};

void calcforces(struct particle *c, int n)
{
	int i,j;
	double dx, dy, fx, fy, qqrrr;
	for (i = 0; i < n; i++) {
		c[i].fx = 0;
		c[i].fy = 0;
	}
	for (i = 0; i < n; i++) {
		for (j = i + 1; j < n; j++) {
			dx = c[j].x - c[i].x;
			dy = c[j].y - c[i].y;
			qqrrr = (c[i].q * c[j].q) / pow(dx * dx + dy * dy, 1.5);
			fx = qqrrr * dx;
			fy = qqrrr * dy;
			c[i].fx -= fx;
			c[i].fy -= fy;
			c[j].fx += fx;
			c[j].fy += fy;
		}
	}
}

void updatepos(struct particle *c, int n, double t)
{
	double ovx, ovy;
	int i;
	for (i = 0; i < n; i++) {
		ovx = c[i].vx;
		ovy = c[i].vy;
		c[i].vx += (c[i].fx / c[i].m) * t;
		c[i].vy += (c[i].fy / c[i].m) * t;
		c[i].x += 0.5 * (ovx + c[i].vx) * t;
		c[i].y += 0.5 * (ovy + c[i].vy) * t;
		if (c[i].x < 0) {
			c[i].x = -c[i].x;
			c[i].vx = -c[i].vx;
		}
		if (c[i].x > WIDTH) {
			c[i].x = 2 * WIDTH - c[i].x;
			c[i].vx = -c[i].vx;
		}

		if (c[i].y < 0) {
			c[i].y = -c[i].y;
			c[i].vy = -c[i].vy;
		}
		if (c[i].y > HEIGHT) {
			c[i].y = 2 * HEIGHT - c[i].y;
			c[i].vy = -c[i].vy;
		}
			
	}
}

void drawpars(SDL_Renderer *ren, struct particle *c, int n)
{
	int i;
	for (i = 0; i < n; i++)
		drawpnt(ren, c[i].x, c[i].y, 16.0/9.0, 16.0/9.0, 16.0/9.0);
}

int main(int argc, char **argv)
{
	if (SDL_Init(SDL_INIT_VIDEO | SDL_INIT_EVENTS | SDL_INIT_TIMER) != 0) {
		fprintf(stderr, "SDL_Init() failed:  %s\n", SDL_GetError());
		return 1;
	}
	SDL_Window *win	= SDL_CreateWindow(
		"points",
		SDL_WINDOWPOS_UNDEFINED,
		SDL_WINDOWPOS_UNDEFINED,
		WIDTH,
		HEIGHT,
		0
	);
	if (win == NULL) {
		fprintf(stderr, "Window creation failed:  %s\n", SDL_GetError());
		return 1;
	}
#ifdef VSYNC
	SDL_Renderer *ren = SDL_CreateRenderer(win, -1, SDL_RENDERER_ACCELERATED | SDL_RENDERER_PRESENTVSYNC);
#else
	SDL_Renderer *ren = SDL_CreateRenderer(win, -1, SDL_RENDERER_ACCELERATED);
#endif
	if (ren == NULL) {
		fprintf(stderr, "Renderer creation failed:  %s\n", SDL_GetError());
		return 1;
	}
	
	struct particle chrgs[MAX_PARS];
	int n = 0;

	Uint64 t = SDL_GetPerformanceCounter(), lt;
	Uint64 freq = SDL_GetPerformanceFrequency();
	int done = 0;
	while (!done) {
		lt = t;
		t = SDL_GetPerformanceCounter();
		SDL_Event e;
		while (SDL_PollEvent(&e)) {
			switch (e.type) {
			case SDL_MOUSEBUTTONUP:
				if (e.button.button == SDL_BUTTON_LEFT) {
					if (n < MAX_PARS) {
						chrgs[n].x = e.button.x;
						chrgs[n].y = e.button.y;
						chrgs[n].m = 1;
						chrgs[n].q = 1;
						chrgs[n].fx = 0;
						chrgs[n].fy = 0;
						chrgs[n].vx = 0;
						chrgs[n].vy = 0;
						n++;
					}
				} else {
					n = 0;
				}
				break;
			case SDL_QUIT:
				done = 1;
			}
		}
		calcforces(chrgs, n);
		updatepos(chrgs, n, TIME_SCALE * ((double)(t-lt))/((double)freq));
		SDL_SetRenderDrawColor(ren, 0, 0, 0, 255);
		SDL_RenderClear(ren);
		SDL_SetRenderDrawBlendMode(ren, SDL_BLENDMODE_ADD);
		drawpars(ren, chrgs, n);
		SDL_RenderPresent(ren);
#ifndef VSYNC
		SDL_Delay(1);
#endif
	}

	SDL_DestroyRenderer(ren);
	SDL_DestroyWindow(win);
	SDL_Quit();
	return 0;
}

Makefile

all:
  gcc -Wall -O2 -o charges charges.c `sdl2-config --cflags --libs`