/*
* sor.c
* (c) Neil Gershenfeld  3/2/03
* Successive Overrelaxation of elliptic boundary value problem
*/

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <stdio.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 alpha,u[SIZE][SIZE];
int i,j,nsteps,step;

main() {
	printf("Number of time steps?  "); scanf("%d",&nsteps);
	printf("Alpha?  "); scanf("%f",&alpha);
	
	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, "SOR 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 = 0; i < SIZE; ++i) {
		u[i][0] = 0.0;
		u[i][SIZE-1] = -1.0;
		u[0][i] = 0;
		u[SIZE-1][i] = 1;
	}

	for (step = 1; step <= nsteps; ++step) {
		for (i = 1; i < (SIZE-1); ++i)
			for (j = 1; j < (SIZE-1); ++j) {
				u[i][j] = (1.0 - alpha) * u[i][j] +
					alpha * 0.25 * (u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1]);
				Data[i][j][0] = (int) (255*(1+u[i][j])/2.0);
				Data[i][j][1] = (int) (128+255*(u[i][j])/2.0);
				Data[i][j][2] = (int) (255*(1-u[i][j])/2.0);
			}
		XPutImage(D, W, Gc, I, 0, 0, 0, 0, SIZE, SIZE);		
	}
}