/*
* diffadi.c
* (c) Neil Gershenfeld  3/2/03
* 2D diffusion by ADI
*/

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define SIZE 500
#define DEPTH 24

Display  *D;
Visual   *V;
int      S;
Window   W;
GC       Gc;
XImage   *I;
XVisualInfo Info;
char Data[SIZE][SIZE][4];

float d,dx,dt,alpha,u[SIZE][SIZE],a[SIZE],b[SIZE],c[SIZE],
	cprime[SIZE], uprime[SIZE];
int j,k,i,n,nsteps,step;

main() {
	printf("Diffusion constant?  "); scanf("%f",&d);
	printf("Space step?  "); scanf("%f",&dx);
	printf("Time step?  "); scanf("%f",&dt);
	printf("Number of time steps?  "); scanf("%d",&nsteps);

	alpha = d * dt / (2.0 * dx * dx);

	D = XOpenDisplay("");
	S = DefaultScreen(D);
	if (XMatchVisualInfo(D, S, DEPTH, TrueColor, &Info) == 0) {
		printf("Display does not support %d bit TrueColor Visual\n",DEPTH);
		return;
	}
	V = Info.visual;
	W = XCreateSimpleWindow(D, DefaultRootWindow(D), 0, 0,
		SIZE, SIZE, 0, 0, 0);       
	XStoreName(D, W, "diffadi output");
	XMapRaised(D, W);
	Gc = XCreateGC (D, W, 0L, (XGCValues *) 0);
	I = XCreateImage(D, V, DEPTH, ZPixmap, 0, (char *) Data,
		SIZE, SIZE, 8, 0);

	for (i = 1; i < (SIZE-1); ++i) {
		a[i] = -alpha;
		b[i] = 1+2.0*alpha;
		c[i] = -alpha;
	}
	a[0] = c[0] = a[SIZE-1] = c[SIZE-1] = 0;
	b[0] = b[SIZE-1] = 1;
	
	for (j = 0; j < SIZE; ++j)
		for (k = 0; k < SIZE; ++k)
			u[j][k] = rand()/((float) RAND_MAX);
			
	for (step = 1; step <= nsteps; ++step) {
		for (k = 1; k < (SIZE-1); ++k) {      	
			cprime[0] = c[0] / b[0];
			uprime[0] = ( (1.0 - 2.0*alpha) * u[j][k] +
			alpha * (u[j][k+1] + u[j][k-1]) ) / b[0];
			for (j = 1; j < SIZE; ++j) {
				cprime[j] = c[j] / (b[j] - a[j] * cprime[j-1]);
				uprime[j] = ( ( (1.0 - 2.0*alpha) * u[j][k] +
				alpha * (u[j][k+1] + u[j][k-1]) ) - a[j] *
				uprime[j-1]) / (b[j] - a[j] * cprime[j-1]);
			}
			for (j = (SIZE-2); j > 0; --j)
				u[j][k] = uprime[j] - cprime[j] * u[j+1][k];
		}
		for (j = 1; j < (SIZE-1); ++j) {      	
			cprime[0] = c[0] / b[0];
			uprime[0] = ( (1.0 - 2.0*alpha) * u[j][k] +
			alpha * (u[j+1][k] + u[j-1][k]) ) / b[0];
			for (k = 1; k < SIZE; ++k) {
				cprime[k] = c[k] / (b[k] - a[k] * cprime[k-1]);
				uprime[k] = ( ( (1.0 - 2.0*alpha) * u[j][k] +
				alpha * (u[j+1][k] + u[j-1][k]) ) - a[k] *
				uprime[k-1]) / (b[k] - a[k] * cprime[k-1]);
			}
			for (k = (SIZE-2); k > 0; --k) {
				u[j][k] = uprime[k] - cprime[k] * u[j][k+1];
				Data[j][k][0] = (int) 255*u[j][k];
				Data[j][k][1] = (int) 255*u[j][k];
				Data[j][k][2] = (int) 255*u[j][k];
			}
		}
		XPutImage(D, W, Gc, I, 0, 0, 0, 0, SIZE, SIZE);
	}
}