/*
* lgshex.c
* (c) Neil Gershenfeld  03/09/03
* simple hexagonal 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 64
#define DEPTH 24

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

/*
* 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, */
	 5, 6, 7, 8,36,		/*  5, 6, 7, 8, 9, */
	10,11,12,13,14,		/* 10,11,12,13,14, */
	15,16,17, 9,19,		/* 15,16,17,18,19, */
	20,42,22,23,24,		/* 20,21,22,23,24, */
	25,26,27,28,29,		/* 25,26,27,28,29, */
	30,31,32,33,34,		/* 30,31,32,33,34, */
	35,18,37,38,39,		/* 35,36,37,38,39, */
	40,41,21,43,44,		/* 40,41,42,43,44, */
	45,46,47,48,49,		/* 45,46,47,48,49, */
	50,51,52,53,54,		/* 50,51,52,53,54, */
	55,56,57,58,59,		/* 55,56,57,58,59, */
	60,61,62,63};		/* 60,61,62,63     */
unsigned char Table1[MAX_BITS] = {
	 0, 1, 2, 3, 4,		/*  0, 1, 2, 3, 4, */
	 5, 6, 7, 8,18,		/*  5, 6, 7, 8, 9, */
	10,11,12,13,14,		/* 10,11,12,13,14, */
	15,16,17,36,19,		/* 15,16,17,18,19, */
	20,42,22,23,24,		/* 20,21,22,23,24, */
	25,26,27,28,29,		/* 25,26,27,28,29, */
	30,31,32,33,34,		/* 30,31,32,33,34, */
	35, 9,37,38,39,		/* 35,36,37,38,39, */
	40,41,21,43,44,		/* 40,41,42,43,44, */
	45,46,47,48,49,		/* 45,46,47,48,49, */
	50,51,52,53,54,		/* 50,51,52,53,54, */
	55,56,57,58,59,		/* 55,56,57,58,59, */
	60,61,62,63};		/* 60,61,62,63     */

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,"Hex 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] = 63;
	while (1) {
	/*****************
	* collision step *
	*****************/
	for (i = 0; i < SIZE; ++i)
		for (j = 0; j < SIZE; j += 2) {
			L2[i+SIZE*j] = Table[L1[i+SIZE*j]];
			L2[i+SIZE*(j+1)] = Table1[L1[i+SIZE*(j+1)]];
		}
	/*******************
	* propagation step *
	*******************/
	for (i = 1; i < (SIZE-1); i += 2)
		for (j = 1; j < (SIZE-1); ++j) {
			L1[i+SIZE*j] = 
				(L2[i+SIZE*(j+1)] & 1)  +
				(L2[(i-1)+SIZE*(j+1)] & 2) +
				(L2[(i-1)+SIZE*j] & 4) +
				(L2[i+SIZE*(j-1)] & 8) +
				(L2[(i+1)+SIZE*j] & 16) +
				(L2[(i+1)+SIZE*(j+1)] & 32);
			L1[(i+1)+SIZE*j] =
				(L2[(i+1)+SIZE*(j+1)] & 1) +
				(L2[i+SIZE*j] & 2) +
				(L2[i+SIZE*(j-1)] & 4) +
				(L2[(i+1)+SIZE*(j-1)] & 8) +
				(L2[(i+2)+SIZE*(j-1)] & 16) +
				(L2[(i+2)+SIZE*j] & 32);
		}
	/*******************************
	* periodic boundary conditions *
	*******************************/
	for (i = 1; i < (SIZE-1); i += 2) {
		L1[i+SIZE*0] = 
			(L2[i+SIZE*(0+1)] & 1)  +
			(L2[(i-1)+SIZE*(0+1)] & 2) +
			(L2[(i-1)+SIZE*0] & 4) +
			(L2[i+SIZE*(SIZE-1)] & 8) +
			(L2[(i+1)+SIZE*0] & 16) +
			(L2[(i+1)+SIZE*(0+1)] & 32);
		L1[(i+1)+SIZE*0] =
			(L2[(i+1)+SIZE*(0+1)] & 1) +
			(L2[i+SIZE*0] & 2) +
			(L2[i+SIZE*(SIZE-1)] & 4) +
			(L2[(i+1)+SIZE*(SIZE-1)] & 8) +
			(L2[(i+2)+SIZE*(SIZE-1)] & 16) +
			(L2[(i+2)+SIZE*0] & 32);
		L1[i+SIZE*(SIZE-1)] = 
			(L2[i+SIZE*0] & 1)  +
			(L2[(i-1)+SIZE*0] & 2) +
			(L2[(i-1)+SIZE*(SIZE-1)] & 4) +
			(L2[i+SIZE*(SIZE-1-1)] & 8) +
			(L2[(i+1)+SIZE*(SIZE-1)] & 16) +
			(L2[(i+1)+SIZE*0] & 32);
		L1[(i+1)+SIZE*(SIZE-1)] =
			(L2[(i+1)+SIZE*0] & 1) +
			(L2[i+SIZE*(SIZE-1)] & 2) +
			(L2[i+SIZE*(SIZE-1-1)] & 4) +
			(L2[(i+1)+SIZE*(SIZE-1-1)] & 8) +
			(L2[(i+2)+SIZE*(SIZE-1-1)] & 16) +
			(L2[(i+2)+SIZE*(SIZE-1)] & 32);
	}
	for (j = 1; j < (SIZE-1); ++j) {
		L1[0+SIZE*j] = 
			(L2[0+SIZE*(j+1)] & 1)  +
			(L2[(SIZE-1)+SIZE*j] & 2) +
			(L2[(SIZE-1)+SIZE*(j-1)] & 4) +
			(L2[0+SIZE*(j-1)] & 8) +
			(L2[(0+1)+SIZE*(j-1)] & 16) +
			(L2[(0+1)+SIZE*j] & 32);
		L1[SIZE-1+SIZE*j] = 
			(L2[SIZE-1+SIZE*(j+1)] & 1)  +
			(L2[(SIZE-1-1)+SIZE*(j+1)] & 2) +
			(L2[(SIZE-1-1)+SIZE*j] & 4) +
			(L2[SIZE-1+SIZE*(j-1)] & 8) +
			(L2[0+SIZE*j] & 16) +
			(L2[0+SIZE*(j+1)] & 32);
	}
	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) | (L1[i+SIZE*j] & 32)) != 0);
		Data[i][j][1] = 255*(((L1[i+SIZE*j] & 2) | (L1[i+SIZE*j] & 4) | (L1[i+SIZE*j] & 8)) != 0);
		Data[i][j][2] = 255*(((L1[i+SIZE*j] & 8) | (L1[i+SIZE*j] & 16) | (L1[i+SIZE*j] & 32)) != 0);		
	}
	XPutImage(D,W,Gc,I,0,0,0,0,SIZE,SIZE);
	}
}