/*
* diffimp.c
* (c) Neil Gershenfeld  3/2/03
* 1D diffusion implicit differences
*/

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

#define NPTS 500

float d=1.0, dx=1.0 ,dt, alpha, u[NPTS], a[NPTS], b[NPTS], c[NPTS],
	cprime[NPTS], uprime[NPTS];
int i, nsteps, step;
XPoint ubuf[NPTS], unewbuf[NPTS];
int S;
Display *D;
Window W;
GC Gc,GcRev;

main() {
	printf("Time step?  "); scanf("%f",&dt);
	printf("Number of time steps?  "); scanf("%d",&nsteps);
	D = XOpenDisplay("");
	S = DefaultScreen(D);
	W = XCreateSimpleWindow(D,DefaultRootWindow(D),0,0,
		NPTS,NPTS,1,WhitePixel(D,S),WhitePixel(D,S));
	XStoreName(D,W,"Implicit Differences");
	XMapRaised(D,W);
	Gc = XCreateGC (D,W,0L,(XGCValues *) 0);
	XSetForeground(D,Gc,BlackPixel(D,S));
	XSetBackground(D,Gc,WhitePixel(D,S));
	XSetLineAttributes(D,Gc,0,LineSolid,CapButt,JoinMiter);
	GcRev = XCreateGC(D,W,0L,(XGCValues *) 0);
	XSetForeground(D,GcRev,WhitePixel(D,S));
	XSetBackground(D,GcRev,BlackPixel(D,S));
	XSetLineAttributes(D,GcRev,0,LineSolid,CapButt,JoinMiter);
	alpha = d * dt / (dx * dx);
	for (i = 0; i < NPTS; ++i) {
		u[i] = 0;
		unewbuf[i].x = ubuf[i].x = i;
		unewbuf[i].y = (int) ((NPTS-1) * (1.0 - (0.1/1.1)));
		a[i] = -alpha;
		b[i] = 1+2.0*alpha;
		c[i] = -alpha;
	}
	u[NPTS/2] = 1.0;
	for (step = 1; step <= nsteps; ++step) {
		cprime[0] = c[0] / b[0];
		uprime[0] = u[0] / b[0];
		for (i = 1; i < NPTS; ++i) {
			cprime[i] = c[i] / (b[i] - a[i] * cprime[i-1]);
			uprime[i] = (u[i] - a[i] * uprime[i-1]) / (b[i] - a[i] * cprime[i-1]);
		}
		u[NPTS-1] = uprime[NPTS-1];
		for (i = (NPTS-2); i >= 0; --i)
			u[i] = uprime[i] - cprime[i] * u[i+1];
		for (i = 0; i < NPTS; ++i) {
			ubuf[i].y = unewbuf[i].y;
			unewbuf[i].y = (int) ((NPTS-1) * (1.0 - ((u[i] + 0.1)/1.1)));
		}
		XDrawLines(D,W,GcRev,ubuf,NPTS,CoordModeOrigin);
		XDrawLines(D,W,Gc,unewbuf,NPTS,CoordModeOrigin);
		XFlush(D);
		usleep(100);
	}
}