/*
* lgssq.c
* (c) Neil Gershenfeld  03/09/03
* simple square lattice gas
*/

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

#define SIZE 500
#define MAX_BITS 16
#define DEPTH 24

#define NumBits(i) ((i&1)+((i&2)>>1)+((i&4)>>2)+((i&8)>>3))

/*
* X variables *
*/
Display *D;
Window  W;
XVisualInfo Info;
Visual  *V;
XSetWindowAttributes Attr;
GC      Gc;
XImage  *I;
char     Data[SIZE][SIZE][4];
unsigned char L1[SIZE*SIZE],L2[SIZE*SIZE];
int     S;

/*
* math variables
*/
int i,j,m,a,c,ran;
unsigned char Table[MAX_BITS] = {
	 0, 1, 2, 3, 4,		/*  0, 1, 2, 3, 4, */
	10, 6, 7, 8, 9,		/*  5, 6, 7, 8, 9, */
	 5,11,12,13,14,		/* 10,11,12,13,14, */
	15};			/* 15              */

main() {
	/***************
	* set up for X *
	***************/
	D = XOpenDisplay("");
	S = DefaultScreen(D);
	if (XMatchVisualInfo(D,S,DEPTH,TrueColor,&Info) == 0) {
		printf("Display can not handle %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,"Square Lattice");
   	XMapRaised(D,W);
        Gc = XCreateGC (D, W, 0L, (XGCValues *) 0);
        I = XCreateImage(D, V, DEPTH, ZPixmap, 0, (char *) Data, SIZE, SIZE, 8, 0);
	/*********************
	* initialize lattice *
	*********************/
	m = 134456;
	a = 8121;
	c = 28411;
	ran = 1;
	for (i = 0; i < SIZE; ++i)
		for (j = 0; j < SIZE; ++j) {
			/* random initial conditions
			ran = (ran*a+c) % m;
			L1[i+SIZE*j] = 2 << ((unsigned char) (6 *
				ran/((double) m)));
			*/
			L1[i+SIZE*j] = 0;
		}
	L1[0] = L1[SIZE-1] = L1[SIZE*(SIZE-1)] = L1[SIZE+SIZE*(SIZE-1)] = 0;
	for (i = 4*SIZE/10; i < 6*SIZE/10; ++i)
		for (j = 4*SIZE/10; j < 6*SIZE/10; ++j)
			L1[i+SIZE*j] = 15;
	while (1) {
	/*****************
	* collision step *
	*****************/
	for (i = 0; i < SIZE; ++i)
		for (j = 0; j < SIZE; ++j) {
			L2[i+SIZE*j] = Table[L1[i+SIZE*j]];
		}
	/*******************
	* propagation step *
	*******************/
	for (i = 1; i < (SIZE-1); ++i)
		for (j = 1; j < (SIZE-1); ++j) {
			L1[i+SIZE*j] = 
				(L2[i+SIZE*(j+1)] & 1)  +
				(L2[(i-1)+SIZE*j] & 2) +
				(L2[i+SIZE*(j-1)] & 4) +
				(L2[(i+1)+SIZE*j] & 8);
		}
	/*******************************
	* periodic boundary conditions *
	*******************************/
	for (i = 1; i < (SIZE-1); ++i) {
		L1[i+SIZE*0] = 
			(L2[i+SIZE*(0+1)] & 1)  +
			(L2[(i-1)+SIZE*0] & 2) +
			(L2[i+SIZE*(SIZE-1)] & 4) +
			(L2[(i+1)+SIZE*0] & 8);
		L1[i+SIZE*(SIZE-1)] = 
			(L2[i+SIZE*0] & 1)  +
			(L2[(i-1)+SIZE*(SIZE-1)] & 2) +
			(L2[i+SIZE*(SIZE-1-1)] & 4) +
			(L2[(i+1)+SIZE*(SIZE-1)] & 8);
	}
	for (j = 1; j < (SIZE-1); ++j) {
		L1[0+SIZE*j] = 
			(L2[0+SIZE*(j+1)] & 1)  +
			(L2[(SIZE-1)+SIZE*j] & 2) +
			(L2[0+SIZE*(j-1)] & 4) +
			(L2[1+SIZE*j] & 8);
		L1[SIZE-1+SIZE*j] = 
			(L2[SIZE-1+SIZE*(j+1)] & 1)  +
			(L2[(SIZE-1-1)+SIZE*j] & 2) +
			(L2[SIZE-1+SIZE*(j-1)] & 4) +
			(L2[0+SIZE*j] & 8);
	}
	for (i = 0; i < SIZE; ++i)
		for (j = 0; j < SIZE; ++j) {
		Data[i][j][0] = 255*(((L1[i+SIZE*j] & 1) | (L1[i+SIZE*j] & 2)) != 0);
		Data[i][j][1] = 255*(((L1[i+SIZE*j] & 2) | (L1[i+SIZE*j] & 4)) != 0);
		Data[i][j][2] = 255*(((L1[i+SIZE*j] & 4) | (L1[i+SIZE*j] & 8) | (L1[i+SIZE*j] & 1)) != 0);
		}
	XPutImage(D,W,Gc,I,0,0,0,0,SIZE,SIZE);
	}
}