Skip to content

Instantly share code, notes, and snippets.

@colejhudson

colejhudson/schrodinger.c

Created November 8, 2018 00:03
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save colejhudson/d13666ea937b8bf8727c132eb7d3ec41 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save colejhudson/d13666ea937b8bf8727c132eb7d3ec41 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
schrodinger.c
/* Schrodinger equation explorer */
/* Michael Creutz */
/* [email protected] */
/* to compile: */
/* cc -O -L/usr/X11R6/lib schrodinger.c -lX11 -lm */
/* this solves the Schrodinger equation in a potential */
/* displays amplitude squared and the potential superposed */
/* the left mouse button perturbs the wave function */
/* the right mouse button perturbs the potential */
/* the damp button evolves system towards the ground state */
/* version of August 2006 */
# include <X11/Xlib.h>
# include <X11/Xutil.h>
# include <X11/Xos.h>
# include <X11/Xatom.h>
# include <X11/cursorfont.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <malloc.h>
# include <math.h>
# include <unistd.h>
/* NSITES=number of lattice sites
fft (as implemented here) REQUIRES a power of two */
# define NSITES 64
# define NPOINTS 128
double repsi[NSITES],impsi[NSITES]; /* real and imaginary parts of psi */
double rekpsi[NSITES],imkpsi[NSITES]; /* field in momentum space */
double potential[NSITES]; /* potential in Schrodinger equation */
double repfactor[NSITES],impfactor[NSITES]; /* for updating potential term */
double rekfactor[NSITES],imkfactor[NSITES]; /* for updating kinetic term */
double cn[NSITES],sn[NSITES]; /* for cosine and sin of angles */
double pi; /* 3.14... */
int iamp[NPOINTS+1],ipot[NPOINTS+1];
XPoint xpoints[2*NPOINTS]; /* for drawing waves */
struct timeval timeout; /* timer for speed control */
/* dimensions for buttons */
# define BUTTONWIDTH 84
# define BUTTONHEIGHT 18
int paused=0; /* is the system paused? */
int damped=0; /* is the system damped? */
int speed=1; /* updating delay proportional to speed; between 0 and 10 */
int hbar=5; /* between 1 and 10 */
int delay=0; /* counter for implementing speed control */
struct timeval timeout; /* timer for speed control */
void drawbutton(),openwindow(),makebuttons(),update(),repaint(),
cleanup(),drawit(),renorm(),oldupdate(),newupdate(),
makefactor(),fastfourier();
# define line(x1,y1,x2,y2) XDrawLine(display,window,gc,x1,y1,x2,y2)
# define color(i) XSetForeground(display,gc,translate[i])
/* various window stuff */
Display *display;
int screen;
static char *progname;
Window window,quitbutton,pausebutton,vconstbutton,
dampbutton,speedbutton,hbarbutton,makebutton();
XColor xcolor,colorcell;
Colormap cmap;
GC gc;
int windowwidth,windowheight,playwidth,playheight;
XFontStruct *font=NULL;
int font_height,font_width;
XSizeHints size_hints;
int darkcolor,lightcolor,black,white;
long translate[256]; /* for converting colors */
int main(argc,argv)
int argc;
char **argv;
{int i,j,k,p;
double psipeg,w;
unsigned int width, height;
XEvent report;
progname=argv[0];
openwindow(argc,argv);
makebuttons();
/* set up trig stuff */
pi=4*atan(1.0);
for (i=0;i<NSITES;i++)
{cn[i]=cos(2*pi*i/(1.*NSITES));
sn[i]=sin(2*pi*i/(1.*NSITES));
}
/* initial configuration */
for (i=0;i<NSITES;i++)
{
repsi[i]=1+cn[i]+cn[(2*i)%NSITES];
impsi[i]=0;
potential[i]=0.5;
}
renorm();
makefactor();
/* loop forever, looking for events */
while(1)
{if ((0==XPending(display))&&(0==paused))
{if (delay) /* don't update yet */
{delay--;
/* this use of select() seems a kludge to me;
why can't usleep() be more standard? */
timeout.tv_sec=0;
timeout.tv_usec=50000; /* .05 sec per delay unit */
select(0,NULL,NULL,NULL,&timeout);
}
else
{delay=speed;
update();
}
}
else
{XNextEvent(display,&report);
switch (report.type)
{case Expose: /* the window is exposed and needs to be redrawn */
if ((report.xexpose.window)!=window) break; /* cuts down flashing,
but you might remove this line if things aren't being redrawn */
if (report.xexpose.count!=0) break; /* more in queue, wait for them */
repaint();
break;
case ConfigureNotify: /* the window has been resized */
width=report.xconfigure.width;
height=report.xconfigure.height;
if ((width<size_hints.min_width)||(height<size_hints.min_height))
{fprintf(stderr,"%s: window too small to proceed.\n",progname);
cleanup();
}
else if ((width!=windowwidth)||(height!=windowheight))
{windowwidth=width;
windowheight=height;
playwidth=windowwidth;
playheight=(windowheight-3*BUTTONHEIGHT);
makebuttons();
}
break;
case ButtonPress: /* which button was pressed? */
if (report.xbutton.window==quitbutton)
cleanup();
else if (report.xbutton.window==pausebutton)
{paused=1-paused;
drawbutton(pausebutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"pause",2-4*paused);
}
else if (report.xbutton.window==dampbutton)
{damped=1-damped;
drawbutton(dampbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"damp",2-4*damped);
makefactor();
}
else if (report.xbutton.window==vconstbutton)
{drawbutton(vconstbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"V=const",-2);
for (i=0;i<NSITES;i++)
{potential[i]=0.5;
}
drawit();
drawbutton(vconstbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"V=const", 2);
makefactor();
}
else if (report.xbutton.window==speedbutton)
{ /* reset speed */
speed=10-(11*report.xbutton.x)/BUTTONWIDTH;
if (speed<0) speed=0;
if (speed>10) speed=10;
delay=speed;
drawbutton(speedbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"speed",-2);
drawbutton(speedbutton,1+((10-speed)*(BUTTONWIDTH-2))/11,1,
(BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
}
else if (report.xbutton.window==hbarbutton)
{ /* reset hbar */
hbar=1+(10*report.xbutton.x)/BUTTONWIDTH;
if (hbar<1) hbar=1;
if (hbar>10) hbar=10;
makefactor();
drawbutton(hbarbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"hbar",-2);
drawbutton(hbarbutton,1+((hbar-1)*(BUTTONWIDTH-2))/10,1,
(BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
}
else
/* button one changes wave function, two updates, and three
modifies potential */
{i=report.xbutton.x;
j=report.xbutton.y;
if (j<playheight) /* add in Lorentzian where clicked */
{p=(i*NSITES)/playwidth; /* site where clicked */
psipeg=5*(j-playheight)/(1.0*playheight);
if (report.xbutton.button==Button1)
{for (k=0;k<NSITES;k++)
{
w=.5*(1+cos(2*pi*(p-k)/(1.*NSITES))); /* peaks at p */
w*=w; /* sharpen the peak */
w*=w;
w*=w;
w*=w;
repsi[k]=psipeg*w+repsi[k]*(1-w);
}
renorm();
}
else if (report.xbutton.button==Button3)
{for (k=0;k<NSITES;k++)
{
w=.5*(1+cos(2*pi*(p-k)/(1.*NSITES))); /* peaks at p */
w*=w; /* sharpen the peak */
w*=w;
w*=w;
w*=w;
potential[k]=-.2*psipeg*w+potential[k]*(1-w);
}
makefactor();
}
else if (report.xbutton.button==Button2) update();
drawit();
}
else update(); /* do a sweep when mouse clicked outside play area*/
}
break;
default:
break;
} /* end of switch */
} /* end of if XPending */
} /* end of while(1) */
} /* end of main */
void update()
{int i,iteration;
double temp;
for (iteration=0;iteration<5;iteration++)
{ /* update the potential term */
for (i=0;i<NSITES;i++)
{
temp =impsi[i]*repfactor[i]+repsi[i]*impfactor[i];
repsi[i]=repsi[i]*repfactor[i]-impsi[i]*impfactor[i];
impsi[i]=temp;
}
/* go to momentum space and update the kinetic term */
fastfourier(repsi,impsi,rekpsi,imkpsi,NSITES,1,1,1);
for (i=0;i<NSITES;i++)
{
temp =imkpsi[i]*rekfactor[i]+rekpsi[i]*imkfactor[i];
rekpsi[i]=rekpsi[i]*rekfactor[i]-imkpsi[i]*imkfactor[i];
imkpsi[i]=temp;
}
/* go back to position space */
fastfourier(rekpsi,imkpsi,repsi,impsi,NSITES,1,1,-1);
}
if (damped) renorm();
drawit();
return;
}
# define size 10
# define dt .05
/* pfactor and kfactor normalize potential and kinetic terms */
/* their only purpose here is to make it easy to change things */
# define pfactor (.6*5./hbar)
# define kfactor (1.0*hbar/(5.*size))
# define dx (size/(1.*NSITES))
# define dampfactor .5
void makefactor() /* make the diagonal matrices used in updating */
{int i,j;
double temp,nsitesinv;
nsitesinv=1.0/NSITES;
for (i=0;i<NSITES;i++)
{temp=pfactor*dt*potential[i];
repfactor[i]=cos(temp);
impfactor[i]=sin(temp);
if (damped)
{temp=exp(-dampfactor*temp);
repfactor[i]*=temp;
impfactor[i]*=temp;
}
j=i;
if (j>=NSITES/2) j-=NSITES;
temp=kfactor*dt*j*j;
rekfactor[i]=nsitesinv*cos(temp);
imkfactor[i]=nsitesinv*sin(temp);
if (damped)
{temp=exp(-dampfactor*temp);
rekfactor[i]*=temp;
imkfactor[i]*=temp;
}
}
return;
}
void renorm() /* normalizes the wave function */
{int i;
double norm;
norm=0.0;
for(i=0;i<NSITES;i++)
norm+=(repsi[i]*repsi[i]+impsi[i]*impsi[i]);
if (0.00001<norm)
{norm=1/sqrt(norm*dx);
for (i=0;i<NSITES;i++)
{repsi[i]*=norm;
impsi[i]*=norm;
}
}
return;
}
#define min(a,b) ((a>b)? b : a)
#define max(a,b) ((a>b)? a : b)
void drawit() /* draws the wave function */
{int i,k1,k2,j,jmax;
double mag;
/* calculate potential and amplitude to display as integers */
j=0;
for (i=0;i<NSITES;i++)
{mag=sqrt(repsi[i]*repsi[i]+impsi[i]*impsi[i]);
k1=playheight*(1-mag);
k2=playheight*(1-potential[i]);
iamp[j]=max(0,min(playheight,k2));
ipot[j]=max(0,min(playheight,k1));
j+=(NPOINTS/NSITES);
}
iamp[NPOINTS]=iamp[0];
ipot[NPOINTS]=ipot[0];
/* interpolate between NSITES to NPOINTS */
jmax=NPOINTS/NSITES;
for (j=1;j<jmax;j++)
for (i=0;i<NSITES;i++)
{
iamp[i*jmax+j]=((jmax-j)*iamp[i*jmax]+j*iamp[(i+1)*jmax])/jmax;
ipot[i*jmax+j]=((jmax-j)*ipot[i*jmax]+j*ipot[(i+1)*jmax])/jmax;
}
/* now draw it */
/* first the red part */
for (i=0;i<NPOINTS;i++)
{
k1=iamp[i];
k2=ipot[i];
xpoints[i].x=xpoints[2*NPOINTS-1-i].x=(i*playwidth)/(NPOINTS-1);
xpoints[i].y=max(k1,k2);
xpoints[2*NPOINTS-1-i].y=max(k2,min(k1,k2)); /* for later use */
}
xpoints[NPOINTS ].x=playwidth;
xpoints[NPOINTS ].y=playheight;
xpoints[NPOINTS+1].x=0;
xpoints[NPOINTS+1].y=playheight;
color(0);
XFillPolygon(display,window,gc,xpoints,NPOINTS+2,Nonconvex,CoordModeOrigin);
/* fix the two points for the magenta part */
for (i=NPOINTS-2;i<NPOINTS;i++)
{
k1=iamp[i];
k2=ipot[i];
xpoints[2*NPOINTS-1-i].x=(i*playwidth)/(NPOINTS-1);
xpoints[2*NPOINTS-1-i].y=max(k2,min(k1,k2));
}
color(12);
XFillPolygon(display,window,gc,xpoints,NPOINTS*2,Nonconvex,CoordModeOrigin);
/* now the green */
for (i=0;i<NPOINTS;i++)
{
k1=iamp[i];
k2=ipot[i];
xpoints[i].y=min(k1,k2);
}
color(2);
XFillPolygon(display,window,gc,xpoints,NPOINTS*2,Nonconvex,CoordModeOrigin);
/* and the cyan on top */
xpoints[NPOINTS ].x=playwidth;
xpoints[NPOINTS ].y=0;
xpoints[NPOINTS+1].x=0;
xpoints[NPOINTS+1].y=0;
color(3);
XFillPolygon(display,window,gc,xpoints,NPOINTS+2,Nonconvex,CoordModeOrigin);
return;
}
void repaint() /* this fixes the window up whenever it is uncovered */
{color(9);
drawbutton(pausebutton, 0,0,BUTTONWIDTH,BUTTONHEIGHT,"pause", 2-4*paused);
drawbutton(dampbutton, 0,0,BUTTONWIDTH,BUTTONHEIGHT,"damp", 2-4*damped);
drawbutton(quitbutton, 0,0,BUTTONWIDTH,BUTTONHEIGHT,"quit", 2);
drawbutton(vconstbutton, 0,0,BUTTONWIDTH,BUTTONHEIGHT,"V=const", 2);
drawbutton(speedbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"speed",-2);
drawbutton(speedbutton,1+((10-speed)*(BUTTONWIDTH-2))/11,1,
(BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
drawbutton(hbarbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"hbar",-2);
drawbutton(hbarbutton,1+((hbar-1)*(BUTTONWIDTH-2))/10,1,
(BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
drawbutton(window,0,playheight+BUTTONHEIGHT,playwidth,BUTTONHEIGHT,
"left mouse button perturbs wave function (psi^2=magenta)",0);
drawbutton(window,0,playheight+2*BUTTONHEIGHT,playwidth,BUTTONHEIGHT,
"right mouse button perturbs potential (green)",0);
color(4);
XFillRectangle(display,window,gc,0,playheight,windowwidth,BUTTONHEIGHT);
drawit();
return;
}
void openwindow(argc,argv)
/* a lot of this is taken from basicwin in the Xlib Programming Manual */
int argc;
char **argv;
{char *window_name="Schrodinger";
char *icon_name="waves";
long event_mask;
Pixmap icon_pixmap;
char *display_name=NULL;
int i;
# define icon_bitmap_width 16
# define icon_bitmap_height 16
static unsigned char icon_bitmap_bits[] = {
0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
/* open up the display */
if ((display=XOpenDisplay(display_name))==NULL)
{fprintf(stderr,"%s: cannot connect to X server %s\n",
progname,XDisplayName(display_name));
exit(-1);
}
screen=DefaultScreen(display);
cmap=DefaultColormap(display,screen);
darkcolor=black=BlackPixel(display,screen);
lightcolor=white=WhitePixel(display,screen);
translate[0]=white;
translate[1]=black;
translate[2]=lightcolor;
translate[3]=darkcolor;
for(i=4;i<256;i++)
translate[i]=translate[i%4];
/* a bunch of colors to play with */
if (XAllocNamedColor(display,cmap,"salmon",&colorcell,&xcolor))
darkcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"wheat",&colorcell,&xcolor))
lightcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"red",&colorcell,&xcolor))
translate[0]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"blue",&colorcell,&xcolor))
translate[1]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"green",&colorcell,&xcolor))
translate[2]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"cyan",&colorcell,&xcolor))
translate[3]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"orange",&colorcell,&xcolor))
translate[4]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"purple",&colorcell,&xcolor))
translate[5]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"yellow",&colorcell,&xcolor))
translate[6]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"pink",&colorcell,&xcolor))
translate[7]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"brown",&colorcell,&xcolor))
translate[8]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"grey",&colorcell,&xcolor))
translate[9]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"turquoise",&colorcell,&xcolor))
translate[10]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"gold",&colorcell,&xcolor))
translate[11]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"magenta",&colorcell,&xcolor))
translate[12]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"navy",&colorcell,&xcolor))
translate[13]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"tan",&colorcell,&xcolor))
translate[14]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"violet",&colorcell,&xcolor))
translate[15]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"white",&colorcell,&xcolor))
translate[16]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"black",&colorcell,&xcolor))
translate[17]=colorcell.pixel;
/* feel free to type in more colors, I got bored */
for(i=18;i<256;i++)
translate[i]=translate[i%18];
/* make the main window */
windowwidth=6*BUTTONWIDTH;
windowheight=300;
playwidth=windowwidth;
playheight=(windowheight-3*BUTTONHEIGHT);
window=XCreateSimpleWindow(display,RootWindow(display,screen),
0,0,windowwidth,windowheight,4,translate[14],lightcolor);
/* make the icon */
icon_pixmap=XCreateBitmapFromData(display,window,
icon_bitmap_bits,icon_bitmap_width,icon_bitmap_height);
size_hints.flags=PPosition | PSize | PMinSize;
size_hints.min_width=windowwidth;
size_hints.min_height=300;
/* this is probably silly, and not really tested, but here goes */
#ifdef X11R3
size_hints.x=x;
size_hints.y=y;
size_hints.width=windowwidth;
size_hints.height=windowheight;
XSetStandardProperties(display,window,window_name,icon_name,
icon_pixmap,argv,argc,&size_hints);
#else
{XWMHints wm_hints;
XClassHint class_hints;
XTextProperty windowName, iconName;
if (XStringListToTextProperty(&window_name,1,&windowName)==0)
{fprintf(stderr,"%s: structure allocation for windowName failed.\n"
,progname);
exit(-1);
}
if (XStringListToTextProperty(&icon_name,1,&iconName)==0)
{fprintf(stderr,"%s: structure allocation for iconName failed.\n"
,progname);
exit(-1);
}
wm_hints.initial_state=NormalState;
wm_hints.input=True;
wm_hints.icon_pixmap=icon_pixmap;
wm_hints.flags=StateHint|IconPixmapHint|InputHint;
class_hints.res_name=progname;
class_hints.res_class="Basicwin";
XSetWMProperties(display,window,&windowName,&iconName,
argv,argc,&size_hints,&wm_hints,&class_hints);
}
#endif
/* pick the events to look for */
event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
XSelectInput(display,window,event_mask);
/* pick font: 9x15 is supposed to almost always be there */
if ((font=XLoadQueryFont(display,"9x15"))==NULL)
if ((font=XLoadQueryFont(display,"fixed"))==NULL)
{fprintf(stderr,"%s: Sorry, having font problems.\n",progname);
exit(-1);
}
font_height=font->ascent+font->descent;
font_width=font->max_bounds.width;
/* make graphics context: */
gc=XCreateGC(display,window,0,NULL);
XSetFont(display,gc,font->fid);
XSetForeground(display,gc,black);
XSetBackground(display,gc,lightcolor);
/* show the window */
XMapWindow(display,window);
return;
}
void makebuttons()
{ /* first destroy any old buttons */
XDestroySubwindows(display,window);
/* now make the new buttons */
quitbutton =makebutton(0,windowheight-3*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);
vconstbutton=makebutton(BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
pausebutton=makebutton(2*BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
dampbutton=makebutton(3*BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
speedbutton=makebutton(windowwidth-BUTTONWIDTH,
windowheight-3*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);
hbarbutton=makebutton(windowwidth-2*BUTTONWIDTH,
windowheight-3*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);
return;
}
void cleanup()
{XUnloadFont(display,font->fid);
XFreeGC(display,gc);
XCloseDisplay(display);
exit(0);
}
Window makebutton(xoffset,yoffset,xsize,ysize)
int xoffset,yoffset,xsize,ysize;
/* Puts button of specified dimensions on window. Nothing is drawn. */
{Window buttonwindow;
long event_mask;
buttonwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
xsize,ysize,0,black,lightcolor);
event_mask=ButtonPressMask|ExposureMask; /* look for mouse-button presses */
XSelectInput(display,buttonwindow,event_mask);
XMapWindow(display,buttonwindow);
return buttonwindow;
}
void drawbutton(buttonwindow,xoffset,yoffset,xsize,ysize,text,state)
Window buttonwindow;
int xoffset,yoffset,xsize,ysize,state;
char * text;
/* Draw a button in buttonwindow of specified dimensions with text
centered. Color of button determined by sign of "state",
size of border determined by magnitude.
I don't need no stinking tookit. */
{int textlength,i,j;
int cdark,clight,cup,cdown;
int cleft,cright,cbutton,ctext;
cup=lightcolor;
cdown=darkcolor;
cdark=black;
clight=white;
if (state<0)
{cbutton=cdown;
ctext=clight;
cleft=cdark;
cright=clight;
}
else
{cbutton=cup;
ctext=cdark;
cleft=clight;
cright=cdark;
}
j=abs(state); /* width for button borders */
XSetForeground(display,gc,cbutton);
XFillRectangle(display,buttonwindow,gc,xoffset+j,yoffset+j,
xsize-2*j,ysize-2*j);
XSetForeground(display,gc,cleft);
XFillRectangle(display,buttonwindow,gc,xoffset,yoffset,xsize,j);
XFillRectangle(display,buttonwindow,gc,xoffset,yoffset,j,ysize);
XSetForeground(display,gc,cright);
for (i=0;i<j;i++)
{XDrawLine(display,buttonwindow,gc,
xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
XDrawLine(display,buttonwindow,gc,
xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
}
if (NULL!=text)
{textlength=strlen(text);
XSetForeground(display,gc,ctext);
XDrawString(display,buttonwindow,gc,
xoffset+(xsize-textlength*font_width)/2,
yoffset+(ysize+font_height)/2-2,
text,textlength);
}
XSetForeground(display,gc,black);
return;
}
void fastfourier(double *rex,double *imx, double *ref,double *imf,
int n,int d1,int d2,int direction)
/* recursively calls itself on two halves of the vector */
/* (rex[],imx[]) is the input vector
(ref[],imf[]) is the output
n=vector length
d1=spacing between elements of input vector
d2= same for output
direction= +/-1 for forward and reverse fourier transforms
vector is not normalized */
{int j,Ndn,nd2,i1,i2;
double c,s,temp1,temp2,temp3;
if (n<=2) /* do two elements case by hand */
{*ref=*rex+*(rex+d1);
*imf=*imx+*(imx+d1);
*(ref+d2)=*rex-*(rex+d1);
*(imf+d2)=*imx-*(imx+d1);
return;
}
nd2=n>>1; /* half the vector length */
/* recursively call fastfourier */
fastfourier(rex ,imx ,ref ,imf ,nd2,2*d1,d2,direction);
fastfourier(rex+d1,imx+d1,ref+nd2*d2,imf+nd2*d2,nd2,2*d1,d2,direction);
/* combine the results */
Ndn=NSITES/n;
for (j=0;j<nd2;j++)
{c=cn[j*Ndn];
s=direction*sn[j*Ndn];
i1=d2*j;
i2=d2*(j+nd2);
temp1 =ref[i1]+c*ref[i2]-s*imf[i2];
temp2 =imf[i1]+c*imf[i2]+s*ref[i2];
temp3 =ref[i1]-c*ref[i2]+s*imf[i2];
imf[i2]=imf[i1]-c*imf[i2]-s*ref[i2];
ref[i1]=temp1;
imf[i1]=temp2;
ref[i2]=temp3;
}
return;
}
Raw
xautomalab.c
/* X-automalab, a totalistic cellular automaton simulator for X windows.
Michael Creutz
[email protected]
compiling: cc -O -o xautomalab -L/usr/X11R6/lib xautomalab.c -lX11
version 2.1
August 2006
minor modifications to eliminate gcc -Wall warnings, March 1999
left-right error depending on bit order fixed, June 1999
This program is still evolving; the latest release can be found via
http://thy.phy.bnl.gov/www/xtoys/xtoys.html
version history
1.1 adds save/restore as a gif file
1.2 allows window resizing
1.3 allows block selection while running and improves picture loading
1.4 supports monochrome
1.5 uses my own gif routines
1.6 handles unusual screen depths and comments things better
1.7 prettier buttons and colors
1.8 fixed a color table error in the gif save routine
2.0 adds speed gadget (it seemed too fast on an alpha)
2.1 minor fixes; lzw restored
*/
# include <X11/Xlib.h>
# include <X11/Xutil.h>
# include <X11/Xos.h>
# include <X11/Xatom.h>
# include <X11/cursorfont.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <malloc.h>
# define BPW (sizeof(long))
/* size and position parameters for buttons, etc. */
# define PLAYTOP 130
# define PLAYLEFT 10
# define BARHEIGHT 45
# define BARWIDTH 120
# define BARTOP 2
# define BARLEFT 2
# define BUTTONWIDTH 112
# define BUTTONHEIGHT 18
# define RIGHTBUTTONS 11
# define RIGHTBOXHEIGHT (BUTTONHEIGHT*RIGHTBUTTONS)
# define XHEIGHT (2*BUTTONHEIGHT)
static char *progname;
int nrows,ncols,volume; /* lattice dimensions (plus two for boundaries) */
/* ncols will be truncated to
a multiple of bytes per long word (BPW),
and then two units are lost for boundaries */
int block, /* size of cells displayed on the screen */
blockshift; /* log_2(block) */
unsigned char *field[2]={NULL,NULL}; /* pointers to old and new system */
int old=0,new=1; /* which field is current? */
int table[2][256]; /* for the rule table */
int nbrhd[3][3]; /* to tag active neighbors */
int birth[9],survivor[9]; /* when are there births or survivors */
/* things for regulating the updating speed */
int speed=0; /* updating delay proportional to speed */
int delay=0; /* counter for implementing speed control */
struct timeval timeout; /* timer for speed control */
int paused=0, /* is the system paused? */
boundary, /* 0,1,2,3 for open, live, periodic, flow */
pencolor=1, /* for sketching */
trace, /* is the tracer on? */
xpast, /* are we to xor with the past? */
neighbors, /* how many live neighbors */
colshift; /* how to shift to get the next column */
long mask0x13,mask0x01,mask0x03,mask0x10; /* some useful masks */
char stringbuffer[256]; /* generally useful */
unsigned char *work=NULL; /* points to a work array */
/* various window stuff */
Display *display;
int screen;
Window window,quitbutton,pausebutton,playground,rulebar,
fillbutton,nhood,rightbox,
srbutton,speedbutton,makebutton();
XColor xcolor,colorcell;
Colormap cmap;
GC gc,gcpen;
XFontStruct *font=NULL;
XSizeHints size_hints;
XImage *spinimage=NULL;
int windowwidth,windowheight;
int font_height,font_width;
int darkcolor,lightcolor,black,white;
long translate[256]; /* for converting colors */
/* text for buttons */
char * righttext[RIGHTBUTTONS]={"open","live","periodic","flow",
"boundary","tracer","clear trace","xor past","","","cell size"};
void settable(),openwindow(),makebuttons(),update(),repaint(),cleanup(),
fixboundary(),showpic(),drawbutton(),setboundary(),
settrace(),setrule(),sketch(),setnbrhd(),savepic(),loadpic(),
setxpast(),putneighbors();
int main(argc,argv)
int argc;
char **argv;
{unsigned int width, height;
int i,j,x,y;
long bitcheck;
XEvent report;
progname=argv[0];
/* initial lattice size */
block=1;
nrows=200/block;
ncols=(~(BPW-1))&(200/block);
volume=nrows*ncols;
/* initial rule */
for (i=0;i<3;i++)
for (j=0;j<3;j++)
nbrhd[i][j]=1;
nbrhd[0][0]=nbrhd[2][1]=0;
trace=1;
xpast=0;
boundary=1;
for (i=0;i<9;i++)
birth[i]=survivor[i]=0;
birth[2]=birth[5]=birth[6]=1;
survivor[4]=survivor[5]=1;
settable();
openwindow(argc,argv);
makebuttons();
/* set the masks for parallel updating */
for (i=0;i<sizeof(long);i++)
{mask0x13=0x13|(mask0x13<<8);
mask0x01=0x01|(mask0x01<<8);
mask0x03=0x03|(mask0x03<<8);
mask0x10=0x10|(mask0x10<<8);
}
/* check for byte order in a word */
bitcheck=1;
if (* (char*) (&bitcheck) )
colshift=(-1); /* LSBfirst */
else
colshift=1; /* MSBfirst */
/* loop forever, looking for events */
while(1)
{if ((0==XPending(display))&&(0==paused))
{if (delay) /* don't update yet */
{delay--;
/* this use of select() seems a kludge to me;
why can't usleep() be more standard? */
timeout.tv_sec=0;
timeout.tv_usec=100000; /* .1 sec per delay unit */
select(0,NULL,NULL,NULL,&timeout);
}
else
{delay=speed;
update();
}
}
else
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if ((report.xexpose.window)!=window) break;
/* cuts down flashing, but you might remove this
if things aren't being redrawn */
if (report.xexpose.count!=0) break; /* more in queue, wait */
repaint();
break;
case ConfigureNotify:
width=report.xconfigure.width;
height=report.xconfigure.height;
if ((width<size_hints.min_width)||(height<size_hints.min_height))
{fprintf(stderr,"%s: window too small to proceed.\n",progname);
cleanup();
}
else if ((width!=windowwidth)||(height!=windowheight))
{windowwidth=width;
windowheight=height;
nrows=(height-PLAYTOP-14)/block;
ncols=(~(BPW-1))&((width-PLAYLEFT-BUTTONWIDTH-20)/block);
volume=nrows*ncols;
XClearWindow(display,window);
makebuttons();
fixboundary();
repaint();
}
break;
case ButtonPress:
if (report.xbutton.window==quitbutton)
cleanup();
else if (report.xbutton.window==pausebutton)
{paused=1-paused;
drawbutton(pausebutton,0,0,68,18,"pause",1-2*paused);
}
else if (report.xbutton.window==fillbutton)
{drawbutton(fillbutton,0,0, 68,18,"clear",-1);
for (i=0;i<volume;i++)
field[old][i]=0;
fixboundary();
showpic();
drawbutton(fillbutton,0,0, 68,18,"clear",1);
}
else if (report.xbutton.window==rightbox){
y=report.xbutton.y/BUTTONHEIGHT;
if (y<4)
{setboundary(y);
showpic();
}
else if (y==5)
settrace(0);
else if (y==6)
settrace(1);
else if (y==7)
setxpast(0);
else if (y==8)
setxpast(1);
else if (y==9)
{i=(1<<(4*report.xbutton.x/BUTTONWIDTH));
if (i!=block) /* reset for new block size */
{block=i;
nrows=(windowheight-PLAYTOP-10)/block;
ncols=(~(BPW-1))&((windowwidth-PLAYLEFT-BUTTONWIDTH-20)/block);
volume=nrows*ncols;
/* clear old window */
XSetForeground(display,gcpen,lightcolor);
XFillRectangle(display,window,gcpen,0,0,
windowwidth,windowheight);
makebuttons();
fixboundary();
showpic();
}
}
}
else if (report.xbutton.window==rulebar)
setrule((report.xbutton.x*(neighbors+1))/BARWIDTH+
((report.xbutton.y*3)/BARHEIGHT)*10);
else if (report.xbutton.window==playground) /* sketch */
{x=report.xbutton.x/block;
y=report.xbutton.y/block;
sketch(x,y);
}
else if (report.xbutton.window==nhood)
{x=(3*report.xbutton.x)/BARHEIGHT;
y=(3*report.xbutton.y)/BARHEIGHT;
setnbrhd(x,y);
}
else if (report.xbutton.window==srbutton)
{if (2*report.xbutton.y<XHEIGHT)
{drawbutton(srbutton,0,0,68,18,"save",-1);
savepic(field[old],ncols,nrows);
drawbutton(srbutton,0,0,68,18,"save",1);
}
else
{drawbutton(srbutton,0,XHEIGHT/2,68,18,"restore",-1);
loadpic(field[old],ncols,nrows);
drawbutton(srbutton,0,XHEIGHT/2,68,18,"restore",1);
}
}
else if (report.xbutton.window==speedbutton)
{ /* reset speed */
speed=10-(11*report.xbutton.x)/BUTTONWIDTH;
if (speed<0) speed=0;
if (speed>10) speed=10;
delay=speed;
drawbutton(speedbutton,0,0,BUTTONWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(BUTTONWIDTH-2))/11,1,
(BUTTONWIDTH-2)/11,16,"",2);
}
else
update(); /* do a sweep when mouse clicked */
break;
default:
break;
} /* end of switch */
} /* end of if XPending */
} /* end of while(1) */
} /* end of main */
void settable() /* set up the rule table from birth and survivor arrays */
/* The use of a table in allows much more general rules
with no loss in speed over these totalistic ones.
Any ideas as to how to add easily selected further rules? */
{int i,j,bits,mask;
/* mask out the unused neighbors */
if (colshift>0) /* correct for byte order */
mask=nbrhd[0][1]
|nbrhd[2][1]<<1
|nbrhd[1][2]<<2
|nbrhd[1][0]<<3
|nbrhd[0][2]<<4
|nbrhd[2][2]<<5
|nbrhd[0][0]<<6
|nbrhd[2][0]<<7;
else
mask=nbrhd[2][1]
|nbrhd[0][1]<<1
|nbrhd[1][2]<<2
|nbrhd[1][0]<<3
|nbrhd[2][2]<<4
|nbrhd[0][2]<<5
|nbrhd[2][0]<<6
|nbrhd[0][0]<<7;
for (i=0;i<256;i++)
{bits=0; /* count the set bits in i */
j=i&mask;
while (j)
{bits+=(1&j);
j>>=1;
}
table[0][i]=birth[bits];
table[1][i]=survivor[bits];
}
return;
}
void sketch(int x, int y) /* sketch in playground starting at point x,y */
{int sketching,i,oldx,oldy,deltax,deltay,newx,newy;
Window root,child;
XEvent sketchreport;
unsigned int keys_buttons;
int window_x,window_y;
newx=x;
newy=y;
field[old][x+ncols*y]^=pencolor;
sketching=1;
while (sketching)
{showpic();
XNextEvent(display,&sketchreport);
switch (sketchreport.type)
{case MotionNotify:
oldx=newx;
oldy=newy;
if (!XQueryPointer(display,playground,
&root,&child,&window_x,&window_y,
&newx,&newy,&keys_buttons))
{sketching=0;
break;
}
if (blockshift)
{newx>>=blockshift;
newy>>=blockshift;
}
if ((newx>ncols)|(newy>nrows)|(newx<0)|(newy<0))
{sketching=0;
break;
}
deltax=abs(newx-oldx);
deltay=abs(newy-oldy);
if (deltax>deltay)
for (i=1;i<=deltax;i++)
field [old][oldx+i*(newx-oldx)/deltax
+ncols*(oldy+i*(newy-oldy)/deltax)]^=pencolor;
else
for (i=1;i<=deltay;i++)
field [old][oldx+i*(newx-oldx)/deltay
+ncols*(oldy+i*(newy-oldy)/deltay)]^=pencolor;
break;
case ButtonRelease:
default:
sketching=0;
break;
}
}
fixboundary();
showpic();
return;
}
void update()
{int i,topi=volume-ncols;
putneighbors((long *) field[old],(long *) field[new]);
for (i=ncols;i<topi;i++) /* bottom and top will be fixed by boundary */
field[new][i]^=table[field[old][i]][work[i]];
old=new;
new=1-new;
fixboundary();
showpic();
return;
}
void putneighbors(long *source, long *destination)
/* initializes the new field and puts neighbors in work array */
/* this updating is done BPW spins at a time using multi-spin coding */
{int i,rowshift=ncols/BPW,topi,cprshift,cmrshift;
unsigned long * temp;
temp=(unsigned long*) work;
topi=(volume/BPW);
/* depending on xpast and trace, set up the destination */
if (xpast&trace)
for (i=0;i<topi;i++)
{destination[i]=(source[i] &mask0x10) /* keep old trace bit */
|((source[i]>>1)&mask0x01) /* previous spin for xpast */
|((source[i]<<1)&mask0x03) /* old state to second bit */
|((source[i]<<3)&mask0x10); /* new trace bit */
source[i]&=mask0x01; /* clean out the source */
}
else if (xpast) /* include previous state into destination */
for (i=0;i<topi;i++)
{destination[i]=(source[i] &mask0x10) /* keep trace bit */
|((source[i]>>1)&mask0x01) /* previous spin for xpast */
|((source[i]<<1)&mask0x03); /* old state to second bit */
source[i]&=mask0x01; /* clean out the source */
}
else if (trace)
for (i=0;i<topi;i++)
{destination[i]=(source[i] &mask0x10) /* keep old trace bit */
|((source[i]<<1)&mask0x03) /* old state to second bit */
|((source[i]<<3)&mask0x10); /* new trace bit */
source[i]&=mask0x01; /* clean out the source */
}
else /* neither trace nor xpast */
for (i=0;i<topi;i++)
{destination[i]=(source[i] &mask0x10) /* keep old trace bit */
|((source[i]<<1)&mask0x03); /* old state to second bit */
source[i]&=mask0x01; /* clean out the source */
}
/* now put the neighbors into the work array */
topi=(volume/BPW)-rowshift;
cprshift=colshift+rowshift;
cmrshift=colshift-rowshift;
for (i=rowshift;i<topi;i++)
{temp[i]=(source[i ]>>8)
|(source[i-colshift]<<(8*(BPW-1))) /* W neighbor */
|(source[i ]<<9)
|(source[i+colshift]>>(8*(BPW-1)-1)) /* E neighbor */
|(source[i+rowshift]<<2 ) /* S neighbor */
|(source[i-rowshift]<<3 ) /* N neighbor */
|(source[i+rowshift]>>4)
|(source[i-cmrshift]<<(8*(BPW-1)+4)) /* SW neighbor */
|(source[i+rowshift]<<13)
|(source[i+cprshift]>>(8*(BPW-1)-5)) /* SE neighbor */
|(source[i-rowshift]>>2)
|(source[i-cprshift]<<(8*(BPW-1)+6)) /* NW neighbor */
|(source[i-rowshift]<<15)
|(source[i+cmrshift]>>(8*(BPW-1)-7));/* NE neighbor */
}
}
void showpic() /* display the field */
{int row,col,i1,i2,color,j,j1,j2,blocktop=block;
char *picture=(*spinimage).data;
if (block>4) blocktop=block-1;
if (8==(*spinimage).depth)
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[old][row*ncols+col]];
j=block*(col+block*ncols*row);
if (color!=picture[j])
for (i1=0;i1<blocktop;i1++)
{j1=i1*block*ncols+j;
for (i2=0;i2<blocktop;i2++)
picture[j1+i2]=color;
}
}
else
{for (j=0;j<volume;j++)
picture[j]=translate[field[old][j]];
}
}
else /* depth is not 8, use xputpixel (this is really ugly) */
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[old][row*ncols+col]];
if (color!=XGetPixel(spinimage,j1=block*col,j2=block*row))
for (i2=0;i2<blocktop;i2++)
for (i1=0;i1<blocktop;i1++)
XPutPixel(spinimage,j1+i1,j2+i2,color);
}
else
for (row=0;row<nrows;row++)
{j1=row*ncols;
for (col=0;col<ncols;col++)
XPutPixel(spinimage,col,row,translate[field[old][j1+col]]);
}
}
XPutImage(display,playground,gc,spinimage,0,0,0,0,block*ncols,block*nrows);
XSync(display,0);
return;
}
void fixboundary()
{int i,bv;
if (boundary<2) /* dead or alive */
{bv=pencolor*boundary;
for (i=0;i<ncols;i++)
{field[old][i]=bv;
field[old][volume-1-i]=bv;
}
for (i=0;i<volume;i+=ncols)
{field[old][i]=bv;
field[old][i+ncols-1]=bv;
}
}
else if (2==boundary) /* periodic */
{for (i=0;i<ncols;i++)
{field[old][i]=field[old][volume-2*ncols+i];
field[old][volume-ncols+i]=field[old][ncols+i];
}
for (i=0;i<volume;i+=ncols)
{field[old][i]=field[old][i+ncols-2];
field[old][i+ncols-1]=field[old][i+1];
}
}
else /* flow: alive on top, dead on sides and bottom */
{for (i=0;i<ncols;i++)
{field[old][i]=pencolor;
field[old][volume-1-i]=0;
}
for (i=0;i<volume;i+=ncols)
{field[old][i]=0;
field[old][i+ncols-1]=0;
}
}
return;
}
void repaint()
/* this fixes the window up whenever it is uncovered */
{int i;
drawbutton(quitbutton,0,0, 68,18,"quit",1);
drawbutton(fillbutton,0,0, 68,18,"clear",1);
drawbutton(srbutton, 0,0, 68,18,"save",1);
drawbutton(srbutton, 0,18,68,18,"restore",1);
drawbutton(pausebutton,0,0,68,18,"pause",1-2*paused);
drawbutton(speedbutton,0,0,BUTTONWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(BUTTONWIDTH-2))/11,1,
(BUTTONWIDTH-2)/11,16,"",2);
drawbutton(window,PLAYLEFT-4,PLAYTOP-4,block*ncols+8,block*nrows+8,NULL,4);
setnbrhd(-1,-1);
/* draw the boundary condition buttons etc */
for (i=0;i<RIGHTBUTTONS;i++)
drawbutton(rightbox,0,i*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT,righttext[i],
((i!=4)&&(i!=10)&&(i!=7)) /* buttons without borders */
*(1-2*((i==boundary)||(i==5)))); /* press boundary and tracer */
setxpast(-1);
/* fix the block buttons */
for (i=1;i<=4;i++)
{
sprintf(stringbuffer,"%d",(1<<i)/2);
drawbutton(rightbox,(i-1)*BUTTONWIDTH/4,9*BUTTONHEIGHT,
BUTTONWIDTH/4,BUTTONHEIGHT,
stringbuffer,1-2*(i==(blockshift+1)));
}
/* write various strings */
XDrawString(display,window,gc,BARWIDTH+BARLEFT+2,BARTOP+font_height-2,
"births",6);
XDrawString(display,window,gc,BARWIDTH+BARLEFT+2,BARTOP+BARHEIGHT-2,
"survivors",9);
XDrawString(display,window,gc,70,BARHEIGHT+BARTOP+2*font_height,
"neighborhood",12);
sprintf(stringbuffer,"%d by %d lattice ",ncols-2,nrows-2);
XDrawImageString(display,window,gc,5,BARHEIGHT+BARTOP+5*font_height,
stringbuffer,strlen(stringbuffer));
if ((*spinimage).depth>1)
{XSetForeground(display,gcpen,translate[1]);
XDrawString(display,window,gcpen,5,BARTOP+BARHEIGHT+4*font_height,
"young",5);
XSetForeground(display,gcpen,translate[3]);
XDrawString(display,window,gcpen,60,BARTOP+BARHEIGHT+4*font_height,
"old",3);
XSetForeground(display,gcpen,translate[2]);
XDrawString(display,window,gcpen,98,BARTOP+BARHEIGHT+4*font_height,
"moribund",8);
}
showpic();
return;
}
void setboundary(int value)
/* take action if boundary buttons hit */
{drawbutton(rightbox,0,(BUTTONHEIGHT*boundary),BUTTONWIDTH,BUTTONHEIGHT,
righttext[boundary],1);
boundary=value;
drawbutton(rightbox,0,(BUTTONHEIGHT*boundary),BUTTONWIDTH,BUTTONHEIGHT,
righttext[boundary],-1);
fixboundary();
return;
}
void settrace(int value)
/* take action if trace buttons hit */
{int i;
if (1==value) /* clear the trace bit */
{drawbutton(rightbox,0,6*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT,
"clear trace",-1);
for (i=0;i<volume;i++)
field[old][i]&=0x3;
showpic();
drawbutton(rightbox,0,6*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT,
"clear trace",1);
}
else
{trace=1-trace;
drawbutton(rightbox,0,5*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT,
"tracer",1-2*trace);
}
return;
}
void setrule(int value)
/* take action if rule buttons hit */
/* call with negative value just to redraw but not change rule */
{int i;
if (value>=0)
{if (value<10)
birth[value]=1-birth[value];
else if (value>=20)
survivor[value-20]=1-survivor[value-20];
}
/* fix the rulebar */
XSetForeground(display,gcpen,lightcolor);
XFillRectangle(display,rulebar,gcpen,0,0,BARWIDTH,BARHEIGHT);
for (i=0;i<(neighbors+1);i++)
{drawbutton(rulebar,i*(BARWIDTH/(neighbors+1)),0,
BARWIDTH/(neighbors+1),BARHEIGHT/3,NULL,1-2*birth[i]);
drawbutton(rulebar,i*(BARWIDTH/(neighbors+1)),2*BARHEIGHT/3,
BARWIDTH/(neighbors+1),BARHEIGHT/3,NULL,1-2*survivor[i]);
sprintf(stringbuffer,"%d",i);
XDrawString(display,rulebar,gc,
18-2*neighbors+BARWIDTH*i/(neighbors+1),2*BARHEIGHT/3-2,
stringbuffer,1);
}
settable();
return;
}
void setnbrhd(int vx, int vy) /* fix the neighborhood buttons */
{int i,j;
neighbors=0;
if ((vx>=0)&(vy>=0)) /* toggle this neighbor */
nbrhd[vx][vy]=1-nbrhd[vx][vy];
nbrhd[1][1]=0;
for (i=0;i<3;i++)
for (j=0;j<3;j++)
{drawbutton(nhood,i*BARHEIGHT/3,j*BARHEIGHT/3,
BARHEIGHT/3,BARHEIGHT/3,NULL,1-2*nbrhd[i][j]);
neighbors+=nbrhd[i][j];
}
setrule(-1); /* redraw the rule bar */
return;
}
void setxpast(int value) /* fix up the xpast button */
{int i;
if ((1==value)&xpast) /* swap past and present */
{drawbutton(rightbox,0,8*BUTTONHEIGHT,BUTTONWIDTH,XHEIGHT/2,"reverse",-1);
for (i=0;i<volume;i++)
field[old][i]=((field[old][i]>>1)&1)|((field[old][i]&1)<<1);
showpic();
drawbutton(rightbox,0,8*BUTTONHEIGHT,BUTTONWIDTH,XHEIGHT/2,"reverse",1);
}
if (0==value) xpast=1-xpast;
/* repaint the button */
if (xpast)
{drawbutton(rightbox,0,7*BUTTONHEIGHT,BUTTONWIDTH,XHEIGHT/2,"xor past",-1);
drawbutton(rightbox,0,8*BUTTONHEIGHT,BUTTONWIDTH,XHEIGHT/2,"reverse",1);
}
else
{drawbutton(rightbox,0,7*BUTTONHEIGHT,BUTTONWIDTH,XHEIGHT/2,"xor past",1);
drawbutton(rightbox,0,8*BUTTONHEIGHT,BUTTONWIDTH,XHEIGHT/2,NULL,0);
}
return;
}
void openwindow(int argc, char **argv)
/* a lot of this is taken from basicwin in the Xlib Programming Manual */
{char *window_name="Automalab";
char *icon_name="Automalab";
long event_mask;
Pixmap icon_pixmap;
char *display_name=NULL;
int i;
# define icon_bitmap_width 16
# define icon_bitmap_height 16
static char icon_bitmap_bits[] = {
0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
/* open up the display */
if ((display=XOpenDisplay(display_name))==NULL)
{fprintf(stderr,"%s: cannot connect to X server %s\n",
progname,XDisplayName(display_name));
exit(-1);
}
screen=DefaultScreen(display);
cmap=DefaultColormap(display,screen);
darkcolor=black=BlackPixel(display,screen);
lightcolor=white=WhitePixel(display,screen);
translate[0]=white;
translate[1]=black;
translate[2]=white;
translate[3]=black;
for(i=4;i<256;i++)
translate[i]=translate[i%4];
/* for monochrome, invert black and white when traced */
translate[8]=black;
translate[17]=white;
translate[18]=black;
translate[19]=white;
if (DefaultDepth(display,screen)>1)
{ /* set the various colors */
if (XAllocNamedColor(display,cmap,"steel blue",&colorcell,&xcolor))
darkcolor=colorcell.pixel; /* pressed button color */
if (XAllocNamedColor(display,cmap,"pale turquoise",&colorcell,&xcolor))
lightcolor=colorcell.pixel; /* background color */
if (XAllocNamedColor(display,cmap,"light gray",&colorcell,&xcolor))
translate[0]=colorcell.pixel; /* dead cells */
if (XAllocNamedColor(display,cmap,"red",&colorcell,&xcolor))
translate[1]=colorcell.pixel; /* newborn */
if (XAllocNamedColor(display,cmap,"sea green",&colorcell,&xcolor))
translate[2]=colorcell.pixel; /* moribund */
if (XAllocNamedColor(display,cmap,"blue",&colorcell,&xcolor))
translate[3]=colorcell.pixel; /* living */
if (XAllocNamedColor(display,cmap,"orange",&colorcell,&xcolor))
translate[4]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"gold",&colorcell,&xcolor))
translate[5]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"yellow",&colorcell,&xcolor))
translate[6]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"white",&colorcell,&xcolor))
translate[7]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"grey",&colorcell,&xcolor))
translate[8]=colorcell.pixel; /* tracing color */
/* fill out the color table for future uses */
for(i=9;i<256;i++)
translate[i]=translate[i%8];
}
translate[16]=translate[8]; /* for tracing */
/* make the main window */
windowwidth=(block*ncols+PLAYLEFT+BUTTONWIDTH+20);
windowheight=(PLAYTOP+block*nrows+14);
window=XCreateSimpleWindow(display,RootWindow(display,screen),
0,0,windowwidth,windowheight,4,black,lightcolor);
/* make the icon */
icon_pixmap=XCreateBitmapFromData(display,window,
icon_bitmap_bits,icon_bitmap_width,icon_bitmap_height);
size_hints.flags=PPosition | PSize | PMinSize;
size_hints.min_width=windowwidth;
size_hints.min_height=windowheight;
#ifdef X11R3
size_hints.x=x;
size_hints.y=y;
size_hints.width=windowwidth;
size_hints.height=windowheight;
XSetStandardProperties(display,window,window_name,icon_name,
icon_pixmap,argv,argc,&size_hints);
#else
{XWMHints wm_hints;
XClassHint class_hints;
XTextProperty windowName, iconName;
if (XStringListToTextProperty(&window_name,1,&windowName)==0)
{fprintf(stderr,"%s: structure allocation for windowName failed.\n"
,progname);
exit(-1);
}
if (XStringListToTextProperty(&icon_name,1,&iconName)==0)
{fprintf(stderr,"%s: structure allocation for iconName failed.\n"
,progname);
exit(-1);
}
wm_hints.initial_state=NormalState;
wm_hints.input=True;
wm_hints.icon_pixmap=icon_pixmap;
wm_hints.flags=StateHint|IconPixmapHint|InputHint;
class_hints.res_name=progname;
class_hints.res_class="Basicwin";
XSetWMProperties(display,window,&windowName,&iconName,
argv,argc,&size_hints,&wm_hints,&class_hints);
}
#endif
/* pick the events to look for */
event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
XSelectInput(display,window,event_mask);
/* pick font: 9x15 is supposed to almost always be there,
the fixed font option put in by Peter Chang */
if ((font=XLoadQueryFont(display,"9x15"))==NULL)
if ((font=XLoadQueryFont(display,"fixed"))==NULL)
{fprintf(stderr,"%s: Sorry, having font problems.\n",progname);
exit(-1);
}
font_height=font->ascent+font->descent;
font_width=font->max_bounds.width;
/* make graphics contexts:
gc for black on background
gcpen for varying purposes */
gc=XCreateGC(display,window,0,NULL);
XSetFont(display,gc,font->fid);
XSetForeground(display,gc,black);
XSetBackground(display,gc,lightcolor);
/* the following might speed things up on some displays? */
XSetPlaneMask(display,gc,black|white|translate[0]|translate[1]|translate[2]
|translate[3]|translate[16]|translate[4]);
gcpen=XCreateGC(display,window,0,NULL);
XSetFont(display,gcpen,font->fid);
XSetForeground(display,gcpen,translate[pencolor]);
XSetBackground(display,gcpen,lightcolor);
/* show the window */
XMapWindow(display,window);
return;
}
void makebuttons()
/* makes buttons and reallocates everything */
{int i,j;
long event_mask;
XEvent report;
Cursor cursor;
/* first destroy any old buttons */
XDestroySubwindows(display,window);
/* now make the new buttons */
quitbutton =makebutton(windowwidth-92,windowheight-42 ,68,18);
pausebutton=makebutton(windowwidth-92,windowheight-60 ,68,18);
fillbutton =makebutton(windowwidth-92,windowheight-78 ,68,18);
srbutton =makebutton(windowwidth-92,windowheight-114,68,36);
speedbutton=makebutton(windowwidth-BUTTONWIDTH-4,
RIGHTBOXHEIGHT+8,BUTTONWIDTH,BUTTONHEIGHT);
rulebar =makebutton(BARLEFT,BARTOP,BARWIDTH,BARHEIGHT);
nhood =makebutton(20,BARHEIGHT+BARTOP+4,BARHEIGHT,BARHEIGHT);
rightbox =makebutton(windowwidth-BUTTONWIDTH-4,2,
BUTTONWIDTH,RIGHTBOXHEIGHT);
playground=XCreateSimpleWindow(display,window,
PLAYLEFT,PLAYTOP,block*ncols,
block*nrows,0,black,translate[4]);
event_mask=ExposureMask|ButtonReleaseMask|ButtonPressMask|
PointerMotionHintMask|ButtonMotionMask;
XSelectInput(display,playground,event_mask);
XMapWindow(display,playground);
/* wait for playground to be displayed befor proceeding */
i=1; /* a flag */
while (i)
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if (report.xexpose.window!=playground) i=0;
default:
break;
}
}
/* make image structure, first destroying any old one */
if (NULL!=spinimage)
{XDestroyImage(spinimage);
spinimage=NULL;
}
spinimage=XGetImage((Display *) display, (Drawable) playground,
0,0,block*ncols,block*nrows,
AllPlanes,ZPixmap);
if (NULL==spinimage)
{fprintf(stderr,"trouble creating image structure\n");
exit(-1);
}
/* make special cursors to be cute */
cursor=XCreateFontCursor(display,XC_sb_up_arrow);
XDefineCursor(display,playground,cursor);
cursor=XCreateFontCursor(display,XC_hand2);
XDefineCursor(display,rightbox,cursor);
XDefineCursor(display,nhood,cursor);
XDefineCursor(display,rulebar,cursor);
XDefineCursor(display,srbutton,cursor);
XDefineCursor(display,quitbutton,cursor);
XDefineCursor(display,pausebutton,cursor);
XDefineCursor(display,fillbutton,cursor);
XDefineCursor(display,speedbutton,cursor);
/* reallocate various arrays */
if (NULL!=work) free((char *) work);
if (NULL==( work=(unsigned char *) malloc(volume)))
{fprintf(stderr,"allocation problems\n");
cleanup();
}
/* clear working array */
for (i=0;i<volume;i++)
work[i]=0;
/* allocate a few bytes of extra space at ends of field array */
for (i=0;i<2;i++)
{if (NULL!=field[i]) free((char *) field[i]-BPW);
if (NULL==( field[i]= (unsigned char *) malloc(volume+2*BPW)))
{fprintf(stderr,"allocation problems\n");
cleanup();
}
/* clear initial state */
for (j=0;j<(volume+2*BPW);j++)
field[i][j]=0;
field[i]+=BPW;
}
/* set blockshift to log_2 of block */
blockshift=0;
i=block;
while(i>>=1)
blockshift++;
/* draw everything */
repaint();
return;
}
Window makebutton(int xoffset, int yoffset, int xsize, int ysize)
/* Puts button of specified dimensions on window. Nothing is drawn. */
{Window buttonwindow;
long event_mask;
buttonwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
xsize,ysize,0,black,lightcolor);
event_mask=ButtonPressMask|ExposureMask; /* look for mouse-button presses */
XSelectInput(display,buttonwindow,event_mask);
XMapWindow(display,buttonwindow);
return buttonwindow;
}
void drawbutton(Window buttonwindow, int xoffset, int yoffset,
int xsize, int ysize, char *text, int state)
/* Draw a button in buttonwindow of specified dimensions with text
centered. Color of button determined by sign of "state."
size of border determined by magnitude of state. */
{int textlength,i,j;
int cdark,clight,cup,cdown;
int cleft,cright,cbutton,ctext;
cup=lightcolor;
cdown=darkcolor;
cdark=black;
clight=white;
if (state<0)
{cbutton=cdown;
ctext=clight;
cleft=cdark;
cright=clight;
}
else
{cbutton=cup;
ctext=cdark;
cleft=clight;
cright=cdark;
}
j=abs(state);
XSetForeground(display,gcpen,cbutton);
XFillRectangle(display,buttonwindow,gcpen,xoffset+j,yoffset+j,
xsize-2*j,ysize-2*j);
XSetForeground(display,gcpen,cleft);
XFillRectangle(display,buttonwindow,gcpen,xoffset,yoffset,xsize,j);
XFillRectangle(display,buttonwindow,gcpen,xoffset,yoffset,j,ysize);
XSetForeground(display,gcpen,cright);
for (i=0;i<j;i++)
{XDrawLine(display,buttonwindow,gcpen,
xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
XDrawLine(display,buttonwindow,gcpen,
xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
}
if (NULL!=text)
{textlength=strlen(text);
XSetForeground(display,gcpen,ctext);
XDrawString(display,buttonwindow,gcpen,
xoffset+(xsize-textlength*font_width)/2,
yoffset+(ysize+font_height)/2-2,
text,textlength);
}
return;
}
void cleanup()
{XUnloadFont(display,font->fid);
XFreeGC(display,gc);
XFreeGC(display,gcpen);
XCloseDisplay(display);
XDestroyImage(spinimage);
if (NULL!=work) free((char *) work);
if (NULL!=field[0]) free((char *) field[0]-BPW);
if (NULL!=field[1]) free((char *) field[1]-BPW);
exit(1);
}
/* from here on is the stuff for saving and restoring from a gif file */
/* for info on how gif works, see ftp://network.ucsd.edu/graphics/GIF.shar.Z */
char *picturename="xautomalab.gif";
void decompress(),compress(),inittable();
int findcode();
void loadpic(data,xsize,ysize) /* load in a GIF image */
unsigned char *data; /* where the output data starts */
int xsize,ysize; /* output picture dimensions */
/* load in a gif image */
{int i,j,filesize,gwidth,gheight,gvolume;
unsigned char *ptr, *ptr1, *rawgif;
int colorbits,codesize;
FILE *infile;
if (NULL==(infile=fopen(picturename,"r")))
{fprintf(stderr,"couldn't open input file\n");
return;
}
/* find the file size */
fseek(infile, 0L, 2);
filesize = ftell(infile);
fseek(infile, 0L, 0);
/* make a place in memory for the file */
if (NULL==(rawgif= (unsigned char *) malloc(filesize) ))
{fprintf(stderr, "not enough memory to read gif file\n");
return;
}
ptr=rawgif;
/* read in the file */
if (fread(ptr, filesize, 1, infile) != 1)
{fprintf(stderr, "read failed\n");
free((char *) rawgif);
return;
}
fclose(infile);
/* check for GIF signature */
if (strncmp((char *) ptr,"GIF87a", 6))
{fprintf(stderr, "not a GIF87a file\n");
free((char *) rawgif);
return;
}
ptr+=6;
ptr+=4; /* skip over screen size */
colorbits=1+((*ptr)&0xf); /* how many bits of color */
ptr+=2; /* skip over background */
if (*ptr++) /* should be zero */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=(3*(1<<colorbits)); /* skip over colormap */
if (','!=(*ptr++)) /* image separator */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=4; /* skip over image offset */
gwidth=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
gheight=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
if (0x80&(*ptr++)) /* skip over local color map */
ptr+=(3*(1<<colorbits));
/* should now be at start of data */
codesize=(*ptr++);
/* make a place for the decompressed file */
gvolume=gwidth*gheight;
ptr1=(unsigned char *) malloc(gvolume);
decompress(codesize,ptr,ptr1,gvolume);
free((char *) rawgif);
/* map picture into data, allowing for different dimensions */
for (j=0;j<ysize;j++)
{if (j>=gheight) break;
for (i=0;i<xsize;i++)
{if (i>=gwidth) break;
data[i+j*xsize]=ptr1[i+j*gwidth];
}
}
free((char *) ptr1);
fixboundary();
showpic();
return;
}
void savepic(unsigned char *data, int xsize, int ysize)
/* save the field as a GIF image */
/* data points to where the input data starts */
{int i;
int colorbits=5,codesize=5; /* assume a 32 color image */
FILE *outfile;
if (NULL==(outfile=fopen(picturename,"w")))
{fprintf(stderr,"couldn't open output file\n");
return;
}
/* GIF signature */
fwrite("GIF87a",6,1,outfile);
/* screen descriptor */
stringbuffer[0]=xsize&0xff; /* screen width */
stringbuffer[1]=(xsize>>8)&0xff;
stringbuffer[2]=ysize&0xff; /* screen height */
stringbuffer[3]=(ysize>>8)&0xff;
stringbuffer[4]=(0x80) /* M=1; global color map follows */
|((colorbits-1)<<4) /* -1+ bits of color reslution */
|(colorbits-1); /* -1+bits per pixel in image */
stringbuffer[5]=0; /* background color */
stringbuffer[6]=0; /* should be zero */
fwrite(stringbuffer,7,1,outfile);
/* global color map */
for (i=0;i<(1<<colorbits);i++)
{colorcell.pixel=translate[i];
XQueryColor(display,cmap,&colorcell);
fputc(colorcell.red>>8,outfile);
fputc(colorcell.green>>8,outfile);
fputc(colorcell.blue>>8,outfile);
}
/* image descriptor */
stringbuffer[0]=','; /* image descriptor separator */
stringbuffer[1]=0; /* image offset */
stringbuffer[2]=0;
stringbuffer[3]=0;
stringbuffer[4]=0;
stringbuffer[5]=xsize&0xff; /* image width */
stringbuffer[6]=(xsize>>8)&0xff;
stringbuffer[7]=ysize&0xff; /* image height */
stringbuffer[8]=(ysize>>8)&0xff;
stringbuffer[9]=0; /* use global color map, no interlace */
fwrite(stringbuffer,10,1,outfile);
/* start of image data */
fputc(codesize,outfile);
compress(codesize,data,outfile,volume);
/* gif terminator */
fputc(';',outfile);
fclose(outfile);
return;
}
/* LZW compression */
/* hash function assumes TABLELENGTH is a power of 2 */
# define TABLELENGTH (1<<13)
char **addresses=NULL; /* where to find the string */
int *codes=NULL, /* the code value */
*linktonext=NULL, /* the next index in the hash chain */
*lengths=NULL, /* the length of the coded string */
*codeindex=NULL; /* the index for a given code */
int nextcode; /* the next unused code */
/* hashit is supposed to give a unique fairly random number in the table for
each length a and string b */
# define hashit(a,b) (51*a+53*(57*b[0]+59*(61*b[a-1]+b[a>>1])))&(TABLELENGTH-1)
void compress(initcodesize,ptr,outfile,size)
int initcodesize; /* the initial compression bits */
char * ptr; /* where the data comes from */
FILE * outfile; /* where the output goes */
int size; /* how much data */
{int currentcode,prefixcode=0,codesize,maxbits=12,maxcode;
int clearcode,eoicode,currentplace=0,length,blocksize=0,bitoffset;
int findcode();
unsigned long outputword;
unsigned char blockbuffer[256]; /* to hold data blocks before writing */
/* allocate space for hash tables */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(addresses=(char **) malloc(sizeof(char *)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* set up initial code table */
inittable(initcodesize);
clearcode=(1<<initcodesize);
eoicode=clearcode+1;
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
/* start with a clear code */
outputword=clearcode;
bitoffset=codesize;
/* now do the compressing */
while (currentplace<size)
{ /* check if codesize needs increasing */
if (nextcode>maxcode)
{codesize++;
maxcode=1<<codesize;
/* if too big, then reset compressor */
if (codesize>maxbits)
{if (bitoffset) outputword|=(clearcode<<bitoffset);
else outputword=clearcode;
bitoffset+=maxbits;
inittable(initcodesize);
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
}
}
/* look for an unstored string */
length=1;
/* the LZW stuff */
while (nextcode>
(currentcode=findcode(length,(char *)(ptr+currentplace))))
{prefixcode=currentcode;
length++;
if ((currentplace+length)>=size) break;
}
nextcode++;
currentplace+=(length-1);
/* output the prefix code */
if (bitoffset) outputword|=(prefixcode<<bitoffset);
else outputword=prefixcode;
bitoffset+=codesize;
/* output finished bytes to blocks */
while (bitoffset>=8)
{blockbuffer[blocksize]=outputword&0xff;
outputword>>=8;
bitoffset-=8;
blocksize++;
/* output filled block */
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
}
/* output the end of information code */
if (bitoffset) outputword|=(eoicode<<bitoffset);
else outputword=eoicode;
bitoffset+=codesize;
/* finish outputting the data */
while (bitoffset>=0)
{blockbuffer[blocksize]=(char) (outputword&0xff);
outputword>>=8;
bitoffset-=8;
blocksize++;
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
/* output the last block */
if (blocksize)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
}
/* a final zero block count */
fputc(0, outfile);
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void decompress(initcodesize,ptr,ptr1,size)
int initcodesize;
unsigned char *ptr, *ptr1; /* compressed data from ptr go to ptr1 */
int size; /* an upper limit purely as a check */
{int i,currentcode,codesize=0,maxbits=12,blocksize;
int clearcode,eoicode,codemask=0;
int bitoffset=0,indx,oldindx=0;
int currentplace=0,oldplace=0;
int findcode();
unsigned long inputword=0;
unsigned char *p1, *p2;
/* first deblock the data */
p1=p2=ptr;
blocksize=(*p1++);
while (blocksize)
{while (blocksize--)
(*p2++)=(*p1++); /* a wonderful example of how abstruse C can be */
blocksize=(*p1++);
}
/* set up initial code table */
currentcode=clearcode=(1<<initcodesize);
eoicode=clearcode+1;
/* allocate space for hash table */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(addresses=(char **) malloc(sizeof(char*)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentplace=eoicode;
}
while (currentcode!=eoicode)
{if (currentcode==clearcode) /* reset the decompressor */
{inittable(initcodesize);
codesize=initcodesize+1;
nextcode=eoicode+1;
codemask=(1<<codesize)-1;
oldindx=(-1);
}
else /* code represents data */
{indx=codeindex[currentcode]; /* where in table is currentcode */
if (indx>=0) /* it is there */
{ /* put it into the output */
for (i=0;i<lengths[indx];i++)
ptr1[currentplace+i]=addresses[indx][i];
if (oldindx>=0) /* first character treated differently */
{findcode(lengths[oldindx]+1,(char *) (ptr1+oldplace));
nextcode++; /* add new code to table */
}
oldplace=currentplace;
currentplace+=lengths[indx];
oldindx=indx;
}
else /* not in table yet; must be old code plus last=first character */
{for (i=0;i<lengths[oldindx];i++)
ptr1[currentplace+i]=addresses[oldindx][i];
ptr1[currentplace+lengths[oldindx]]=addresses[oldindx][0];
/* store new code */
findcode(lengths[oldindx]+1,(char *)(ptr1+currentplace));
oldplace=currentplace;
currentplace+=lengths[oldindx]+1;
oldindx=codeindex[nextcode++];
}
/* crude error checking */
if ((oldindx<0)||(currentplace>size))
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* check if codesize needs increasing */
if (nextcode>codemask)
if (codesize<maxbits)
{codesize++;
codemask=(1<<codesize)-1;
}
while (bitoffset<codesize) /* read some more data */
{if (bitoffset) inputword|=(((int)(*ptr++))<<bitoffset);
else inputword=(*ptr++);
bitoffset+=8;
}
/* strip off current code */
currentcode=inputword&codemask;
inputword>>=codesize;
bitoffset-=codesize;
if (currentcode>nextcode)
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void inittable(size)
int size;
{int i;
for (i=0;i<TABLELENGTH;i++)
{linktonext[i]=(-1);
codes[i]=(-1);
lengths[i]=(-1);
}
for (i=0;i<4096;i++)
codeindex[i]=(-1);
/* store initial codes for raw characters */
nextcode=0;
for (i=0;i<(1<<size);i++)
{stringbuffer[i]=i;
findcode(1,(char *) (stringbuffer+i));
nextcode++;
}
return;
}
int findcode(length,string)
char *string;
int length;
/* return code for string of given length;
if not found, store it and return nextcode */
{int i,j,indx,previousindex;
indx=hashit(length,string);
/* look for string in table */
previousindex=indx;
while (indx>0)
{if (lengths[indx]==length) /* is the length right ? */
{for (j=0;j<length;j++)
if ((string[j]) != (addresses[indx][j])) break;
if (j==length) return codes[indx];
}
previousindex=indx;
indx=linktonext[indx];
}
/* not found, so store it */
indx=previousindex;
i=indx;
while (codes[i]>=0) /* find an unused slot in table */
{i++;
if (i>=TABLELENGTH) i-=TABLELENGTH;
}
if (i!=indx)
{linktonext[indx]=i; /* link to it */
indx=i; /* move to it */
}
codes[indx]=nextcode; /* save the new code */
lengths[indx]=length;
addresses[indx]=string;
codeindex[nextcode]=indx;
return nextcode;
}
Raw
xautomalab.txt
XAutomaLab
by
Michael Creutz
Physics Department
Brookhaven National Laboratory
Upton, NY 11973
[email protected]
Play God over your own universe. With AutomaLab you control the
microphysics of a discrete world of cellular automata. This
two-dimensional land is an array of cells (initially 198 by 198)
displayed as pixels on an X-window screen. Each cell can be alive or
empty, with evolution occurring in discrete time steps. Empty cells
are grey, newborn ones are red, older living ones are blue, and cells
that have just died are green. When the tracer is turned on, old life
leaves a legacy by shading the background color.
In one time step, the fate of each cell depends on the number of its
living neighbors. Using Boolean gadgets, you control when a new live
cell will be born on an empty site and when a living one will survive.
Another gadget lets you pick which neighbors, up to the eight nearest
neighbors, are counted.
To be more precise, if a particular birth gadget is set and if that
number of selected neighbors to a dead cell are alive, then that cell
will become alive on the next time step. If a survivor gadget is set,
then a living cell with that number of live neighbors will survive the
updating.
With mouse selections you determine whether the boundary of your
universe is dead, alive, or periodic. A final boundary choice, called
"flow," is alive on the top and dead on the sides and bottom.
The quit, pause and clear gadgets should be self explanatory. When
the system is paused, a click outside the playfield or any specific
gadget will update the system a single step. Finally, you can toggle
individual cells on and off by pressing a mouse button and sketching
over your world.
The save button saves the current configuration as "xautomalab.gif," a
standard gif file that you can print or manipulate with any graphics
program that likes gif files. The restore gadget will reload a
previously stored configuration. The rule used to create a stored
configuration is not itself stored or restored. Note that for the
save feature to work you need write permission on the directory from
which you are running.
Since sketching when the cells are small is rather imprecise, the
save/restore buttons can be useful for creating special initial
configurations. Using the big block size, sketch your desired
configuration while the system is paused. Save it, and then switch to
a smaller block size. Finally reload the picture at this new
resolution. When your fancy spaceship works, be sure to save a copy
of the image.
The restore function will take any valid gif picture. Try copying
some unrelated picture to xautomalab.gif and then hit the restore
button. You can then watch your favorite rule dissolve the picture.
The cell size buttons determine how many screen pixels are used per
cell. This defaults to one. The window can also be resized to give
different lattice sizes.
With the eight neighborhood there are 18 birth/survivor buttons. This
means there are 2^18=262,144 possible rules in this case. Other
neighborhoods give many more, all selectable with the mouse. Indeed,
it is unlikely that you will be able to try them all.
In addition, the number of possible universes is further doubled using
the "xor past" gadget, which probably needs some explanation. When
this is activated, the new state is finally XOR'ed with the history
one time step back. Thus if a cell was alive in the past, the new
state is the opposite of what the birth and survivor gadgets want.
The reason for doing this is to produce reversible rules. If the
history and current states are interchanged, the system will go
backwards through the sequence of configurations from which it came.
An analogy is reversing all the momenta of a bunch of atoms. This
interchange is accomplished by the "reverse" gadget, which
interchanges young and moribund cells. Without the Xpast gadget set,
the reverse gadget does not appear.
If this seems confusing, try this: With the Xpast gadget set, clear
the screen, draw some small picture, select a random rule from the
birth and survivor gadgets, and let the system evolve until the screen
becomes a mess. Then hit the reverse button and watch the initial
picture reappear. Try repeating this experiment, but alter a single
pixel with the mouse at the time of the reversal.
The starting rule is one I particularly liked at the time I last
updated the program. It may change with further updates. Watch it
for a while and then try other things.
A particularly well known rule is Conway's classic cellular automaton
model "life." This uses the eight cell neighborhood and a new cell is
born for exactly three live neighbors, while a living cell dies with
less than 2 (lonely) or more than 4 (overcrowding) neighbors.
Another well-known rule is Fredkin's modulo two model which uses the 4
cell neighborhood. Here a state flips if it has an odd number of
active neighbors and is unchanged otherwise. Start with some small
picture and observe how the initial state is replicated.
If you want to try a rule running from a random start, run for a while
with some chaotic rule (i.e. most rules with births on one neighbor)
and then switch to your rule of choice.
Xautomalab is based on my earlier Amiga program Automalab, which
appeared on the May 1991 issue of Jumpdisk. That version used direct
access the Amiga graphics chips for speed. Comparing this X version
to the previous illustrates the awesome power of the Amiga Blitter.
If you have an Amiga and want a copy of Automalab, check in the www
site mentioned below.
The program is written in C and should compile on anything supporting
xwindows. It is best on a color display. With monochrome, only the
living cells show.
To compile "cc -O -o xautomalab xautomalab.c -lX11" should do it,
unless your X11 libraries are in a peculiar place.
This program and some related things are available via
"http://thy.phy.bnl.gov/www/xtoys.html" on the World Wide Web. It
is still evolving, and the latest version can be found there.
The colors are defined in the subroutine openwindow. An easy way to
change them is to use the chart of X11 color names, which should be in
any Xwindow programming guide.
This program with source is freely distributable.
Raw
xfire.c
/* Xfires, a forest fire simulator for X windows.
Michael Creutz
[email protected]
compiling: cc -O -o xfires -L/usr/X11R6/lib xfires.c -lX11
version of August 2006
The latest version is kept at "http://thy.phy.bnl.gov/www/xtoys/xtoys.html"
A text description of this program is there as well.
*/
# include <X11/Xlib.h>
# include <X11/Xutil.h>
# include <X11/Xos.h>
# include <X11/Xatom.h>
# include <X11/cursorfont.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <malloc.h>
/* size and position parameters for buttons, etc. */
# define PLAYTOP 116
# define PLAYLEFT 14
# define BUTTONWIDTH 68
# define BBWIDTH 96
# define SRHEIGHT (2*18)
# define BPW (sizeof(long))
/* lattice dimensions; ncols will be truncated to
a multiple of bytes per long word (BPW),
and then two units are lost for boundaries */
int nrows,ncols,volume,mask;
int block, /* size of cells displayed on the screen */
blockshift; /* log_2(block) */
unsigned char *field[2]={NULL,NULL}; /* pointers to old and new system */
int old=0,new=1; /* which field is current? */
# define barren 0
# define tree 1
# define fire 2
int paused=0, /* is the system paused? */
iteration=0; /* counting sweeps */
long mrand48(),lrand48(); /* generates a random word, lrand is non-negative */
double drand48(); /* generates a random double */
/* for fast but crude random number generation */
# define RSIZE 127
long births[RSIZE];
int randomizer=57,birthindex=0;
char stringbuffer[256]; /* generally useful */
/* things for regulating the updating speed */
int speed=0; /* updating delay proportional to speed */
int delay=0; /* counter for implementing speed control */
struct timeval timeout; /* timer for speed control */
/* various window stuff */
Display *display;
int screen;
static char *progname;
Window window,quitbutton,pausebutton,playground,blockbutton,srbutton,
speedbutton,makebutton();
XColor xcolor,colorcell;
Colormap cmap;
GC gc,gcpen;
int windowwidth,windowheight;
XFontStruct *font=NULL;
int font_height,font_width;
XSizeHints size_hints;
int darkcolor,lightcolor,black,white;
XImage *spinimage=NULL;
long translate[256]; /* for converting colors */
char * srtext[2]={"save","restore"};
void drawbutton(),openwindow(),makebuttons(),update(),repaint(),
cleanup(),showpic(),fixboundary(),loadpic(),savepic(),compress(),
decompress(),inittable();
int main(argc,argv)
int argc;
char **argv;
{unsigned int width, height,x,y;
int i;
XEvent report;
progname=argv[0];
/* initial lattice size */
block=1;
nrows=200/block;
ncols=(~(BPW-1))&(200/block);
volume=nrows*ncols;
/* set up array for fast random number generation */
for (i=0;i<RSIZE;i++)
births[i]=(lrand48()>>2);
/* open the window and make the buttons */
openwindow(argc,argv);
makebuttons();
/* loop forever, looking for events */
while(1)
{if (0==paused)
{if (delay) /* don't update yet */
{delay--;
/* this use of select() seems a kludge to me;
why can't usleep() be more standard? */
timeout.tv_sec=0;
timeout.tv_usec=100000; /* .1 sec per delay unit */
select(0,NULL,NULL,NULL,&timeout);
}
else
{delay=speed;
update();
}
}
while (paused|XPending(display))
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if ((report.xexpose.window)!=window) break; /* cuts down flashing,
but you might remove this line if things aren't being redrawn */
if (report.xexpose.count!=0) break; /* more in queue, wait for them */
repaint();
break;
case ConfigureNotify:
width=report.xconfigure.width;
height=report.xconfigure.height;
if ((width<size_hints.min_width)||(height<size_hints.min_height))
{fprintf(stderr,"%s: window too small to proceed.\n",progname);
cleanup();
}
else if ((width!=windowwidth)||(height!=windowheight))
{windowwidth=width;
windowheight=height;
XClearWindow(display,window);
nrows=(height-PLAYTOP-10)/block;
ncols=(~(BPW-1))&((width-PLAYLEFT-16)/block);
volume=nrows*ncols;
makebuttons();
showpic();
}
break;
case ButtonPress:
if (report.xbutton.window==quitbutton)
cleanup();
else if (report.xbutton.window==pausebutton)
{paused=1-paused;
drawbutton(pausebutton,0,0,BUTTONWIDTH,18,"pause",1-2*paused);
}
else if (report.xbutton.window==blockbutton)
{i=(1<<(4*report.xbutton.x/BBWIDTH));
if (i!=block) /* reset for new block size */
{block=i;
nrows=(windowheight-PLAYTOP-10)/block;
ncols=(~(BPW-1))&((windowwidth-PLAYLEFT-16)/block);
volume=nrows*ncols;
XSetForeground(display,gcpen,lightcolor);
XFillRectangle(display,window,gcpen,0,0,
windowwidth,windowheight);
makebuttons();
showpic();
}
}
else if (report.xbutton.window==playground) /* start a fire */
{x=report.xbutton.x/block;
y=report.xbutton.y/block;
field[old][x+ncols*y]=field[new][x+ncols*y]=fire;
fixboundary();
showpic();
}
else if (report.xbutton.window==srbutton) /* save or restore */
{y=(report.xbutton.y*2)/SRHEIGHT;
drawbutton(srbutton,0,y*(SRHEIGHT/2),BUTTONWIDTH,SRHEIGHT/2,
srtext[y],-1);
if (0==y)
savepic(field[old],ncols,nrows);
else if (1==y) /* load in gif image */
{loadpic(field[old],ncols,nrows);
for (i=0;i<volume;i++)
field[new][i]=1&(field[old][i]&=3);
showpic();
}
drawbutton(srbutton,0,y*(SRHEIGHT/2),BUTTONWIDTH,SRHEIGHT/2,
srtext[y],1);
}
else if (report.xbutton.window==speedbutton)
{ /* reset speed */
speed=10-(11*report.xbutton.x)/BBWIDTH;
if (speed<0) speed=0;
if (speed>10) speed=10;
delay=speed;
drawbutton(speedbutton,0,0,BBWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(BBWIDTH-2))/11,1,
(BBWIDTH-2)/11,16,"",2);
}
else
update(); /* do a sweep when mouse clicked */
break;
default:
break;
} /* end of switch */
} /* end of if XPending */
} /* end of while(1) */
} /* end of main */
void update()
{int i,nfires,ntrees,newtrees=volume/32;
unsigned char *pold,*pnew,*pup,*pdown,*pleft,*pright,*ptop;
/* grow new trees */
while (newtrees)
{newtrees--;
i=volume;
while (i>=volume)
{i=mask&births[birthindex];
births[birthindex]^=births[randomizer];
if ((++birthindex)>=RSIZE) birthindex=0;
if ((++randomizer)>=RSIZE) randomizer=0;
}
if (field[old][i]!=fire) field[new][i]=field[old][i]=tree;
}
/* spread fires */
pnew=field[new]+ncols;
ptop=field[old]+volume-ncols;
for (pold=field[old]+ncols; pold<ptop; pold++,pnew++)
{if (fire==*pold)
{pright=pnew+1;
pleft=pnew-1;
pup=pnew-ncols;
pdown=pnew+ncols;
*pright=fire & (*pright | (*pright<<1));
*pleft =fire & (*pleft | (*pleft <<1));
*pup =fire & (*pup | (*pup <<1));
*pdown =fire & (*pdown | (*pdown <<1));
*pold=*pnew=barren; /* old fire dies */
}
}
/* spread fires from top and bottom boundaries */
for (i=0;i<ncols;i++)
if (fire==field[old][i])
{if (tree==field[new][i+ncols]) field[new][i+ncols]=fire;
field[old][i]=field[new][i]=barren;
}
for (i=volume-ncols;i<volume;i++)
if (fire==field[old][i])
{if (tree==field[new][i-ncols]) field[new][i-ncols]=fire;
field[old][i]=field[new][i]=barren;
}
old=new;
new=1-new;
fixboundary();
showpic();
iteration++;
if (25<=iteration)
{ntrees=nfires=0;
for (i=ncols;i<volume-ncols;i++)
{ntrees+=(field[old][i]==tree);
nfires+=(field[old][i]==fire);
}
iteration=0;
sprintf(stringbuffer,"%d fires, %d trees ",nfires,ntrees);
XDrawImageString(display,window,gc,8,82,stringbuffer,strlen(stringbuffer));
if (0==nfires) field[old][(int) (volume*drand48())]=fire;
}
return;
}
void showpic() /* display the field */
{int row,col,i1,i2,color,j,j1,j2,blocktop=block;
char *picture=(*spinimage).data;
if (block>4) blocktop=block-1;
if (8==(*spinimage).depth)
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[old][row*ncols+col]];
j=block*(col+block*ncols*row);
if (color!=picture[j])
for (i1=0;i1<blocktop;i1++)
{j1=i1*block*ncols+j;
for (i2=0;i2<blocktop;i2++)
picture[j1+i2]=color;
}
}
else
{for (j=0;j<volume;j++)
picture[j]=translate[field[old][j]];
}
}
else /* depth is not 8, use xputpixel (this is really ugly) */
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[old][row*ncols+col]];
if (color!=XGetPixel(spinimage,j1=block*col,j2=block*row))
for (i2=0;i2<blocktop;i2++)
for (i1=0;i1<blocktop;i1++)
XPutPixel(spinimage,j1+i1,j2+i2,color);
}
else
for (row=0;row<nrows;row++)
{j1=row*ncols;
for (col=0;col<ncols;col++)
XPutPixel(spinimage,col,row,translate[field[old][j1+col]]);
}
}
XPutImage(display,playground,gc,spinimage,0,0,0,0,block*ncols,block*nrows);
XSync(display,0);
return;
}
void fixboundary()
/* copies edges for periodicity */
{int i;
for (i=0;i<ncols;i++)
{field[old][i]=field[old][volume-2*ncols+i];
field[old][volume-ncols+i]=field[old][ncols+i];
}
for (i=0;i<volume;i+=ncols)
{field[old][i]=field[old][i+ncols-2];
field[old][i+ncols-1]=field[old][i+1];
}
return;
}
void repaint()
/* this fixes the window up whenever it is uncovered */
{int i;
drawbutton(pausebutton,0,0,BUTTONWIDTH,18,"pause",1-2*paused);
drawbutton(quitbutton,0,0,BUTTONWIDTH,18,"quit",1);
drawbutton(window,PLAYLEFT-4,PLAYTOP-4,block*ncols+8,block*nrows+8,NULL,4);
drawbutton(speedbutton,0,0,BBWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(BBWIDTH-2))/11,1,
(BBWIDTH-2)/11,16,"",2);
for (i=0;i<2;i++)
drawbutton(srbutton,0,i*(SRHEIGHT/2),BUTTONWIDTH,SRHEIGHT/2,srtext[i],1);
/* fix the block buttons */
for (i=1;i<=4;i++)
{sprintf(stringbuffer,"%d",(1<<i)/2);
drawbutton(blockbutton,(i-1)*BBWIDTH/4,0,BBWIDTH/4,18,
stringbuffer,1-2*(i==(blockshift+1)));
}
/* write various strings */
XDrawString(display,window,gc,24,62,"cell size",9);
sprintf(stringbuffer,"%d by %d lattice ",ncols-2,nrows-2);
XDrawImageString(display,window,gc,14,100,
stringbuffer,strlen(stringbuffer));
XDrawString(display,window,gc,windowwidth-40,108
,"MJC",3);
showpic();
return;
}
void openwindow(argc,argv)
/* a lot of this is taken from basicwin in the Xlib Programming Manual */
int argc;
char **argv;
{char *window_name="Forest fires";
char *icon_name="fires";
long event_mask;
Pixmap icon_pixmap;
char *display_name=NULL;
int i;
# define icon_bitmap_width 16
# define icon_bitmap_height 16
static unsigned char icon_bitmap_bits[] = {
0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
/* open up the display */
if ((display=XOpenDisplay(display_name))==NULL)
{fprintf(stderr,"%s: cannot connect to X server %s\n",
progname,XDisplayName(display_name));
exit(-1);
}
screen=DefaultScreen(display);
cmap=DefaultColormap(display,screen);
darkcolor=black=BlackPixel(display,screen);
lightcolor=white=WhitePixel(display,screen);
translate[0]=white;
translate[1]=black;
translate[2]=lightcolor;
translate[3]=darkcolor;
if (XAllocNamedColor(display,cmap,"firebrick",&colorcell,&xcolor))
darkcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"wheat",&colorcell,&xcolor))
lightcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"black",&colorcell,&xcolor))
translate[barren]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"forest green",&colorcell,&xcolor))
translate[tree]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"yellow",&colorcell,&xcolor))
translate[fire]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"grey",&colorcell,&xcolor))
translate[3]=colorcell.pixel;
/* fill out the color table for future uses */
for(i=4;i<256;i++)
translate[i]=translate[i%4];
/* make the main window */
windowwidth=(block*ncols+PLAYLEFT+16);
windowheight=(PLAYTOP+block*nrows+10);
window=XCreateSimpleWindow(display,RootWindow(display,screen),
0,0,windowwidth,windowheight,4,black,lightcolor);
/* make the icon */
icon_pixmap=XCreateBitmapFromData(display,window,
icon_bitmap_bits,icon_bitmap_width,icon_bitmap_height);
/* hints for window manager */
size_hints.flags=PPosition | PSize | PMinSize;
size_hints.min_width=windowwidth;
size_hints.min_height=windowheight-160;
#ifdef X11R3
size_hints.x=x;
size_hints.y=y;
size_hints.width=windowwidth;
size_hints.height=windowheight;
XSetStandardProperties(display,window,window_name,icon_name,
icon_pixmap,argv,argc,&size_hints);
#else
{XWMHints wm_hints;
XClassHint class_hints;
XTextProperty windowName, iconName;
if (XStringListToTextProperty(&window_name,1,&windowName)==0)
{fprintf(stderr,"%s: structure allocation for windowName failed.\n"
,progname);
exit(-1);
}
if (XStringListToTextProperty(&icon_name,1,&iconName)==0)
{fprintf(stderr,"%s: structure allocation for iconName failed.\n"
,progname);
exit(-1);
}
wm_hints.initial_state=NormalState;
wm_hints.input=True;
wm_hints.icon_pixmap=icon_pixmap;
wm_hints.flags=StateHint|IconPixmapHint|InputHint;
class_hints.res_name=progname;
class_hints.res_class="Basicwin";
XSetWMProperties(display,window,&windowName,&iconName,
argv,argc,&size_hints,&wm_hints,&class_hints);
}
#endif
/* pick the events to look for */
event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
XSelectInput(display,window,event_mask);
/* pick font: 9x15 is supposed to almost always be there */
if ((font=XLoadQueryFont(display,"9x15"))==NULL)
if ((font=XLoadQueryFont(display,"fixed"))==NULL)
{fprintf(stderr,"%s: Sorry, having font problems.\n",progname);
exit(-1);
}
font_height=font->ascent+font->descent;
font_width=font->max_bounds.width;
/* make graphics contexts:
gc for black on white
gcpen for varying purposes */
gc=XCreateGC(display,window,0,NULL);
XSetFont(display,gc,font->fid);
XSetForeground(display,gc,black);
XSetBackground(display,gc,lightcolor);
gcpen=XCreateGC(display,window,0,NULL);
XSetFont(display,gcpen,font->fid);
XSetForeground(display,gcpen,darkcolor);
XSetBackground(display,gcpen,lightcolor);
/* show the window */
XMapWindow(display,window);
return;
}
void makebuttons()
{int i;
long event_mask;
XEvent report;
Cursor cursor;
/* first destroy any old buttons */
XDestroySubwindows(display,window);
/* now make the new buttons */
quitbutton =makebutton(4,4,BUTTONWIDTH,18);
pausebutton =makebutton(BUTTONWIDTH+8,4,BUTTONWIDTH,18);
blockbutton =makebutton(16,32,BBWIDTH,18);
srbutton =makebutton(windowwidth-BUTTONWIDTH-4,4,BUTTONWIDTH,18*2);
speedbutton=makebutton(windowwidth-BBWIDTH-8,
46,
BBWIDTH,18);
playground=XCreateSimpleWindow(display,window,
PLAYLEFT,PLAYTOP,block*ncols,block*nrows,0,black,translate[3]);
event_mask=ExposureMask|ButtonReleaseMask|ButtonPressMask|
PointerMotionHintMask|ButtonMotionMask;
XSelectInput(display,playground,event_mask);
XMapWindow(display,playground);
/* wait for playground to be displayed befor proceeding */
i=1; /* a flag */
while (i)
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if (report.xexpose.window!=playground) i=0;
default:
break;
}
}
/* make image structure */
if (NULL!=spinimage)
{XDestroyImage(spinimage);
spinimage=NULL;
}
spinimage=XGetImage((Display *) display, (Drawable) playground,
0,0,block*ncols,block*nrows,AllPlanes,ZPixmap);
if (NULL==spinimage)
{fprintf(stderr,"trouble creating image structure\n");
exit(-1);
}
/* make special cursors to be cute */
cursor=XCreateFontCursor(display,XC_sb_up_arrow);
XDefineCursor(display,playground,cursor);
cursor=XCreateFontCursor(display,XC_hand2);
XDefineCursor(display,blockbutton,cursor);
XDefineCursor(display,quitbutton,cursor);
XDefineCursor(display,pausebutton,cursor);
XDefineCursor(display,srbutton,cursor);
XDefineCursor(display,speedbutton,cursor);
/* reallocate various arrays */
for (i=0;i<2;i++)
{if (NULL!=field[i]) free((char *) field[i]);
if (NULL==( field[i]= (unsigned char *) malloc(volume)))
{fprintf(stderr,"allocation problems\n");
cleanup();
}
}
/* set blockshift to log_2 of block */
blockshift=0;
i=block;
while(i>>=1)
blockshift++;
/* mask used for random site selection */
mask=1;
while (mask<volume) mask=1|(mask<<1);
/* set initial state random */
for (i=0;i<volume;i++)
field[old][i]=field[new][i]=lrand48()&tree;
/* draw everything */
repaint();
return;
}
Window makebutton(xoffset,yoffset,xsize,ysize)
int xoffset,yoffset,xsize,ysize;
/* Puts button of specified dimensions on window. Nothing is drawn. */
{Window buttonwindow;
long event_mask;
buttonwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
xsize,ysize,0,black,lightcolor);
event_mask=ButtonPressMask|ExposureMask; /* look for mouse-button presses */
XSelectInput(display,buttonwindow,event_mask);
XMapWindow(display,buttonwindow);
return buttonwindow;
}
void drawbutton(buttonwindow,xoffset,yoffset,xsize,ysize,text,state)
Window buttonwindow;
int xoffset,yoffset,xsize,ysize,state;
char * text;
/* Draw a button in buttonwindow of specified dimensions with text
centered. Color of button determined by sign of "state."
size of border determined by magnitude. */
{int textlength,i,j;
int cdark,clight,cup,cdown;
int cleft,cright,cbutton,ctext;
cup=lightcolor;
cdown=darkcolor;
cdark=black;
clight=white;
if (state<0)
{cbutton=cdown;
ctext=clight;
cleft=cdark;
cright=clight;
}
else
{cbutton=cup;
ctext=cdark;
cleft=clight;
cright=cdark;
}
j=abs(state);
XSetForeground(display,gcpen,cbutton);
XFillRectangle(display,buttonwindow,gcpen,xoffset+j,yoffset+j,
xsize-2*j,ysize-2*j);
XSetForeground(display,gcpen,cleft);
XFillRectangle(display,buttonwindow,gcpen,xoffset,yoffset,xsize,j);
XFillRectangle(display,buttonwindow,gcpen,xoffset,yoffset,j,ysize);
XSetForeground(display,gcpen,cright);
for (i=0;i<j;i++)
{XDrawLine(display,buttonwindow,gcpen,
xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
XDrawLine(display,buttonwindow,gcpen,
xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
}
if (NULL!=text)
{textlength=strlen(text);
XSetForeground(display,gcpen,ctext);
XDrawString(display,buttonwindow,gcpen,
xoffset+(xsize-textlength*font_width)/2,
yoffset+(ysize+font_height)/2-2,
text,textlength);
}
return;
}
void cleanup()
{XUnloadFont(display,font->fid);
XFreeGC(display,gc);
XFreeGC(display,gcpen);
XCloseDisplay(display);
XDestroyImage(spinimage);
if (NULL!=field[0]) free((char *) field[0]);
if (NULL!=field[1]) free((char *) field[1]);
exit(1);
}
/* from here on is the stuff for saving and restoring from a gif file */
/* for info on how gif works, see ftp://network.ucsd.edu/graphics/GIF.shar.Z */
char *picturename="xfires.gif";
void loadpic(data,xsize,ysize) /* load GIF image to field */
unsigned char *data; /* where the output data starts */
int xsize,ysize; /* output picture dimensions */
{int i,j,filesize,gwidth,gheight,gvolume;
unsigned char *ptr, *ptr1, *rawgif;
int colorbits,codesize;
FILE *infile;
if (NULL==(infile=fopen(picturename,"r")))
{fprintf(stderr,"couldn't open input file\n");
return;
}
/* find the file size */
fseek(infile, 0L, 2);
filesize = ftell(infile);
fseek(infile, 0L, 0);
/* make a place in memory for the file */
if (NULL==(rawgif= (unsigned char *) malloc(filesize) ))
{fprintf(stderr, "not enough memory to read gif file\n");
return;
}
ptr=rawgif;
/* read in the file */
if (fread(ptr, filesize, 1, infile) != 1)
{fprintf(stderr, "read failed\n");
free((char *) rawgif);
return;
}
fclose(infile);
/* check for GIF signature */
if (strncmp((char *) ptr,"GIF87a", 6))
{fprintf(stderr, "not a GIF87a file\n");
free((char *) rawgif);
return;
}
ptr+=6;
ptr+=4; /* skip over screen size */
colorbits=1+((*ptr)&0xf); /* how many bits of color */
ptr+=2; /* skip over background */
if (*ptr++) /* should be zero */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=(3*(1<<colorbits)); /* skip over colormap */
if (','!=(*ptr++)) /* image separator */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=4; /* skip over image offset */
gwidth=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
gheight=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
if (0x80&(*ptr++)) /* skip over local color map */
ptr+=(3*(1<<colorbits));
/* should now be at start of data */
codesize=(*ptr++);
/* make a place for the decompressed file */
gvolume=gwidth*gheight;
ptr1=(unsigned char *) malloc(gvolume);
decompress(codesize,ptr,ptr1,gvolume);
free((char *) rawgif);
/* map picture into data, allowing for different dimensions */
for (j=0;j<ysize;j++)
{if (j>=gheight) break;
for (i=0;i<xsize;i++)
{if (i>=gwidth) break;
data[i+j*xsize]=ptr1[i+j*gwidth];
}
}
free((char *) ptr1);
fixboundary();
return;
}
void savepic(data,xsize,ysize) /* save the field as a GIF image */
unsigned char *data; /* where the input data starts */
int xsize,ysize; /* picture dimensions */
{int i;
int colorbits=5,codesize=5; /* assume a 32 color image */
FILE *outfile;
if (NULL==(outfile=fopen(picturename,"w")))
{fprintf(stderr,"couldn't open output file\n");
return;
}
/* GIF signature */
fwrite("GIF87a",6,1,outfile);
/* screen descriptor */
stringbuffer[0]=xsize&0xff; /* screen width */
stringbuffer[1]=(xsize>>8)&0xff;
stringbuffer[2]=ysize&0xff; /* screen height */
stringbuffer[3]=(ysize>>8)&0xff;
stringbuffer[4]=(0x80) /* M=1; global color map follows */
|((colorbits-1)<<4) /* -1+ bits of color reslution */
|(colorbits-1); /* -1+bits per pixel in image */
stringbuffer[5]=0; /* background color */
stringbuffer[6]=0; /* should be zero */
fwrite(stringbuffer,7,1,outfile);
/* global color map */
for (i=0;i<(1<<colorbits);i++)
{colorcell.pixel=translate[i];
XQueryColor(display,cmap,&colorcell);
fputc(colorcell.red>>8,outfile);
fputc(colorcell.green>>8,outfile);
fputc(colorcell.blue>>8,outfile);
}
/* image descriptor */
stringbuffer[0]=','; /* image descriptor separator */
stringbuffer[1]=0; /* image offset */
stringbuffer[2]=0;
stringbuffer[3]=0;
stringbuffer[4]=0;
stringbuffer[5]=xsize&0xff; /* image width */
stringbuffer[6]=(xsize>>8)&0xff;
stringbuffer[7]=ysize&0xff; /* image height */
stringbuffer[8]=(ysize>>8)&0xff;
stringbuffer[9]=0; /* use global color map, no interlace */
fwrite(stringbuffer,10,1,outfile);
/* start of image data */
fputc(codesize,outfile);
compress(codesize,data,outfile,volume);
/* gif terminator */
fputc(';',outfile);
fclose(outfile);
return;
}
/* LZW compression */
/* hash function assumes TABLELENGTH is a power of 2 */
# define TABLELENGTH (1<<13)
char **addresses=NULL; /* where to find the string */
int *codes=NULL, /* the code value */
*linktonext=NULL, /* the next index in the hash chain */
*lengths=NULL, /* the length of the coded string */
*codeindex=NULL; /* the index for a given code */
int nextcode; /* the next unused code */
/* hashit is supposed to give a unique fairly random number in the table for
each length a and string b */
# define hashit(a,b) (51*a+53*(57*b[0]+59*(61*b[a-1]+b[a>>1])))&(TABLELENGTH-1)
void compress(initcodesize,ptr,outfile,size)
int initcodesize; /* the initial compression bits */
char * ptr; /* where the data comes from */
FILE * outfile; /* where the output goes */
int size; /* how much data */
{int currentcode,prefixcode=0,codesize,maxbits=12,maxcode;
int clearcode,eoicode,currentplace=0,length,blocksize=0,bitoffset;
int findcode();
unsigned long outputword;
unsigned char blockbuffer[256]; /* to hold data blocks before writing */
/* allocate space for hash tables */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(addresses=(char **) malloc(sizeof(char *)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* set up initial code table */
inittable(initcodesize);
clearcode=(1<<initcodesize);
eoicode=clearcode+1;
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
/* start with a clear code */
outputword=clearcode;
bitoffset=codesize;
/* now do the compressing */
while (currentplace<size)
{ /* check if codesize needs increasing */
if (nextcode>maxcode)
{codesize++;
maxcode=1<<codesize;
/* if too big, then reset compressor */
if (codesize>maxbits)
{if (bitoffset) outputword|=(clearcode<<bitoffset);
else outputword=clearcode;
bitoffset+=maxbits;
inittable(initcodesize);
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
}
}
/* look for an unstored string */
length=1;
while (nextcode>
(currentcode=findcode(length,(char *)(ptr+currentplace))))
{prefixcode=currentcode;
length++;
if ((currentplace+length)>=size) break;
}
nextcode++;
currentplace+=(length-1);
/* output the prefix code */
if (bitoffset) outputword|=(prefixcode<<bitoffset);
else outputword=prefixcode;
bitoffset+=codesize;
/* output finished bytes to blocks */
while (bitoffset>=8)
{blockbuffer[blocksize]=outputword&0xff;
outputword>>=8;
bitoffset-=8;
blocksize++;
/* output filled block */
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
}
/* output the end of information code */
if (bitoffset) outputword|=(eoicode<<bitoffset);
else outputword=eoicode;
bitoffset+=codesize;
/* finish outputting the data */
while (bitoffset>=0)
{blockbuffer[blocksize]=(char) (outputword&0xff);
outputword>>=8;
bitoffset-=8;
blocksize++;
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
/* output the last block */
if (blocksize)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
}
/* a final zero block count */
fputc(0, outfile);
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void decompress(initcodesize,ptr,ptr1,size)
int initcodesize;
unsigned char *ptr, *ptr1; /* compressed data from ptr go to ptr1 */
int size; /* an upper limit purely as a check */
{int i,currentcode,codesize=0,maxbits=12,blocksize;
int clearcode,eoicode,codemask=0;
int bitoffset=0,indx,oldindx=0;
int currentplace=0,oldplace=0;
int findcode();
unsigned long inputword=0;
unsigned char *p1, *p2;
/* first deblock the data */
p1=p2=ptr;
blocksize=(*p1++);
while (blocksize)
{while (blocksize--)
(*p2++)=(*p1++); /* a wonderful example of how abstruse C can be */
blocksize=(*p1++);
}
/* set up initial code table */
currentcode=clearcode=(1<<initcodesize);
eoicode=clearcode+1;
/* allocate space for hash table */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(addresses=(char **) malloc(sizeof(char*)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentplace=eoicode;
}
while (currentcode!=eoicode)
{if (currentcode==clearcode) /* reset the decompressor */
{inittable(initcodesize);
codesize=initcodesize+1;
nextcode=eoicode+1;
codemask=(1<<codesize)-1;
oldindx=(-1);
}
else /* code represents data */
{indx=codeindex[currentcode]; /* where in table is currentcode */
if (indx>=0) /* it is there */
{ /* put it into the output */
for (i=0;i<lengths[indx];i++)
ptr1[currentplace+i]=addresses[indx][i];
if (oldindx>=0) /* first character treated differently */
{findcode(lengths[oldindx]+1,(char *) (ptr1+oldplace));
nextcode++; /* add new code to table */
}
oldplace=currentplace;
currentplace+=lengths[indx];
oldindx=indx;
}
else /* not in table yet; must be old code plus last=first character */
{for (i=0;i<lengths[oldindx];i++)
ptr1[currentplace+i]=addresses[oldindx][i];
ptr1[currentplace+lengths[oldindx]]=addresses[oldindx][0];
/* store new code */
findcode(lengths[oldindx]+1,(char *)(ptr1+currentplace));
oldplace=currentplace;
currentplace+=lengths[oldindx]+1;
oldindx=codeindex[nextcode++];
}
/* crude error checking */
if ((oldindx<0)||(currentplace>size))
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* check if codesize needs increasing */
if (nextcode>codemask)
if (codesize<maxbits)
{codesize++;
codemask=(1<<codesize)-1;
}
while (bitoffset<codesize) /* read some more data */
{if (bitoffset) inputword|=(((int)(*ptr++))<<bitoffset);
else inputword=(*ptr++);
bitoffset+=8;
}
/* strip off current code */
currentcode=inputword&codemask;
inputword>>=codesize;
bitoffset-=codesize;
if (currentcode>nextcode)
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void inittable(size)
int size;
{int i,findcode();
for (i=0;i<TABLELENGTH;i++)
{linktonext[i]=(-1);
codes[i]=(-1);
lengths[i]=(-1);
}
for (i=0;i<4096;i++)
codeindex[i]=(-1);
/* store initial codes for raw characters */
nextcode=0;
for (i=0;i<(1<<size);i++)
{stringbuffer[i]=i;
findcode(1,(char *) (stringbuffer+i));
nextcode++;
}
return;
}
int findcode(length,string)
char *string;
int length;
/* return code for string of given length;
if not found, store it and return nextcode */
{int i,j,indx,previousindex;
indx=hashit(length,string);
/* look for string in table */
previousindex=indx;
while (indx>0)
{if (lengths[indx]==length) /* is the length right ? */
{for (j=0;j<length;j++)
if ((string[j]) != (addresses[indx][j])) break;
if (j==length) return codes[indx];
}
previousindex=indx;
indx=linktonext[indx];
}
/* not found, so store it */
indx=previousindex;
i=indx;
while (codes[i]>=0) /* find an unused slot in table */
{i++;
if (i>=TABLELENGTH) i-=TABLELENGTH;
}
if (i!=indx)
{linktonext[indx]=i; /* link to it */
indx=i; /* move to it */
}
codes[indx]=nextcode; /* save the new code */
lengths[indx]=length;
addresses[indx]=string;
codeindex[nextcode]=indx;
return nextcode;
}
Raw
xfire.txt
Forest Fires
by
Michael Creutz
[email protected]
This program simulates forest fires.
On each site of a lattice is either
nothing, a tree, or a fire.
In one time step a fire spreads to
adjacent trees and leaves an
empty space. Trees are born
in a random matter with a probability
of approximately 1/32 per time step.
If no fires are active, one is
started at a random location. Fires
can also be started with the
mouse button.
An Amiga version of this program was
published in the Dec. 1993 issue of
Jumpdisk.
This program is freely distributable.
Raw
xising.c
/*
* xising
*
* Xwindow program to simulate the 2-d Ising model. The underlying
* approach uses many "demons" circulating around the lattice trying
* to flip spins.
*
* The latest version of this program can be found at
* http://thy.phy.bnl.gov/www/xtoys/xtoys.html
*
* Michael Creutz [email protected]
*
November 1999, version 2.11: corrects problems on alphas pointed out
by Kari Rummukainen
*
* compiling: cc -O -o xising -L/usr/X11R6/lib xising.c -lm -lX11
*
* On suns:
*
* cc -O -I/usr/openwin/include -L/usr/openwin/lib newxising.c -lm -lX11
*
* If you find some machine supporting X that this does not work on,
* please let me know.
*
*/
/* the usual includes */
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <X11/Xatom.h>
#include <X11/cursorfont.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
/*
* initial lattice dimensions: ncols will be truncated to a multiple
* of WORDSIZE nrows will be changed if a multiple of 29 since I let
* the demons hop 29 sites to get a bit of randomness
*/
int nrows = 256, ncols = 256;
long **field = NULL; /* to store the system; it will be
* dynamically allocated to a [ncols]
* by [nrows/WORDSIZE] array */
/* some convenient definitions */
#define WORDSIZE (CHAR_BIT*sizeof(long))
#define LEFTBIT (~LONG_MAX)
#define nextrow (row+1-nrows*(row>=nrows-1))
/* nshift is here to allow for different endianness */
#define nextcol (col+nshift*(1-hwords*(col==topcol)))
#define previousrow (row-1+nrows*(row<1))
#define previouscol (col-nshift*(1-hwords*(col==bottomcol)))
char stringbuffer[100];
static char *progname;
char bits[256]; /* for counting bits */
/* hwords will store the number of long words in a single row */
int hwords, volume;
/* bitcount counts the set bits in a word; call setup() before using */
#define bitcount8(i) (bits[(i)&255])
#define bitcount16(i) (bitcount8(i)+bitcount8((i)>>8))
/* check if long int is 32 bits */
#if ((0x7FFFFFFF)==(LONG_MAX)) /* 32 bits */
#define bitcount(i) (bitcount16(i)+bitcount16((i)>>16))
#define RND (random()^(random()<<1))
#else /* assume 64 bits */
#define bitcount32(i) (bitcount16(i)+bitcount16((i)>>16))
#define bitcount(i) (bitcount32(i)+bitcount32((i)>>32))
#define RND (random()^(random()<<16)^(random()<<40))
#endif
/* things to monitor temperature, etc. */
double beta; /* the inverse temperature */
long energycount = 0, volumecount = 0;
int heater = 0, cooler = 0;
int canonical = 1;
int boundary = 0, algorithm = 0, paused = 0, iteration = 0;
/* variables to allow for peculiar bit orderings on some machines */
int nshift, topcol, bottomcol;
/* for local algorithm; use WORDSIZE two-bit demons */
long demon0, demon1, checkermask;
long *work0 = NULL, *work1 = NULL;
/*
* bitprob represents the probability for the corresponding demon bit
* to be set multiplied by 1<<16 and rounded to an integer; used to
* speed up demon refreshing
*/
long bitprob0, bitprob1;
/* stuff for the cluster algorithm; use three demon bits */
#define DEMONBITS 3
int ndemons;
long **cluster = NULL; /* to store cluster shape */
long **sadx = NULL, **sady = NULL; /* labels sad demons */
int **unchecked = NULL; /* to reduce redundant calculations */
int boxtop, boxbottom; /* will enclose cluster */
long **demon = NULL; /* for demons in cluster algorithm */
long **newdemon = NULL;
int *activerow = NULL;
long demonindex = 0; /* for shuffling demons */
long clustergrowth = 0;
int direction = 1; /* for vertical sweeps through
* cluster */
/* things for regulating the updating speed */
int speed = 0; /* updating delay proportional to
* speed */
int delay = 0; /* counter for implementing speed
* control */
struct timeval timeout; /* timer for speed control */
/* various window stuff */
Display *display;
int screen;
Window window, quitbutton, pausebutton, algobutton, heatwindow,
boundwindow, speedbutton, canonbutton, makebutton();
GC gc, gcpen, gcclear;
XImage *spinimage = NULL, *clusterimage = NULL;
XFontStruct *font = NULL;
int font_height, font_width;
unsigned int windowwidth, windowheight;
XSizeHints size_hints;
int depth, darkcolor, lightcolor, black, white, spinup, spindown;
long event_mask;
/* all my functions */
void refreshdemons(int step); /* refreshes the cluster demons, step
* steps over some of them for only a
* partial refresh */
void drawbutton(), openwindow(), setup(), makebuttons(), growcluster(),
startcluster(), localupdate(), repaint(), setalg(),
setheater(), setboundary(), fixthermometer(), check(),
leftshift(), rightshift(), allocarray(), heatbath();
/* text for the buttons */
char *boundtext[4] = {"periodic", "antiperiodic",
"up boundary", "down"};
char *heatertext[3] = {"run free", "heat", "cool"};
char *algotext[2] = {"local", "cluster"};
/* general layout dimensions */
#define LEFTOFFSET 64
#define TOPOFFSET 124
#define HEATERHEIGHT 54
#define HEATERWIDTH 86
#define BOUNDHEIGHT 72
#define BOUNDWIDTH 116
#define BUTTONWIDTH 48
#define BUTTONHEIGHT 18
#define TTOP 128
#define THGT (nrows*2)
#define TBOTTOM (TTOP+THGT)
#define TLEFT 24
#define TWDTH 5
int main(argc, argv)
int argc;
char **argv;
{
unsigned int width, height, i;
XEvent report;
progname = argv[0];
openwindow(argc, argv);
setup();
makebuttons();
/* loop forever, looking for events */
while (1) {
if (clustergrowth) /* finish cluster regardless of
* algorithm */
while (clustergrowth) {
growcluster();
}
else if ((0 == XPending(display)) && (0 == paused)){
/* no events waiting */
if (delay) { /* don't update yet */
delay--;
/*
* this use of select() seems a kludge to me; why can't
* usleep() be more standard?
*/
timeout.tv_sec = 0;
timeout.tv_usec = 50000;/* .05 sec per delay unit */
select(0, NULL, NULL, NULL, &timeout);
} else {
delay = speed; /* reset delay counter */
if (1 == algorithm) {
startcluster();
} else
localupdate();
}
}
while (paused || XPending(display)) {
XNextEvent(display, &report); /* find out what X wants done */
switch (report.type) {
case Expose: /* the window has been exposed and
* needs redrawing */
if ((report.xexpose.window) != window)
break; /* cuts down flashing, but you might
* remove this line if things aren't
* being redrawn */
if (report.xexpose.count != 0)
break; /* more in queue, wait for them */
repaint();
break;
case ConfigureNotify: /* the window has been resized */
width = report.xconfigure.width;
height = report.xconfigure.height;
if ((width < size_hints.min_width) ||
(height < size_hints.min_height)) {
fprintf(stderr, "%s: window too small to proceed.\n", progname);
XUnloadFont(display, font->fid); /* I'm still too neat */
XFreeGC(display, gc); /* from my Amiga days */
XFreeGC(display, gcpen);
XFreeGC(display, gcclear);
XCloseDisplay(display);
exit(1);
}
if ((width != windowwidth) || (height != windowheight)) {
windowwidth = width;
windowheight = height;
/* new lattice dimensions */
ncols = width - 75;
nrows = (height - 164 - 24) / 2;
makebuttons();
}
break;
case ButtonPress:
if (report.xbutton.window == quitbutton) {
XUnloadFont(display, font->fid);
XFreeGC(display, gc);
XFreeGC(display, gcpen);
XFreeGC(display, gcclear);
XCloseDisplay(display);
exit(1);
} else if (report.xbutton.window == pausebutton) {
paused = 1 - paused;
drawbutton(pausebutton, 0, 0, BUTTONWIDTH, BUTTONHEIGHT,
"pause", 1 - 2 * paused);
} else if (report.xbutton.window == algobutton) {
/* toggle algorithm */
algorithm = 1 - algorithm;
/* 1 for cluster, 0 for local updating */
setalg();
for (i = 0; i < 2; i++)
drawbutton(algobutton, 0, i * 18, 68, 18,
algotext[i], 1 - 2 * (algorithm == i));
} else if (report.xbutton.window == heatwindow) {
/* adjust heater */
setheater((report.xbutton.y * 3) / HEATERHEIGHT);
} else if (report.xbutton.window == boundwindow)
/* set boundary conditions */
{
setboundary((report.xbutton.y * 4) / BOUNDHEIGHT);
} else if (report.xbutton.window == canonbutton) {
canonical = 1 - canonical;
if (canonical)
drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT,
"canonical", -1);
else
drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT,
"microcanonical", 1);
} else if (report.xbutton.window == speedbutton) {
/*
* reset speed; speed=0 is the fastest, 10 the slowest
*/
speed = 10 - (11 * (report.xbutton.x - 1)) / (BOUNDWIDTH - 2);
if (speed < 0)
speed = 0;
if (speed > 10)
speed = 10;
delay = speed;
drawbutton(speedbutton, 0, 0, BOUNDWIDTH, 18, "speed", -1);
drawbutton(speedbutton, 1 + ((10 - speed) * (BOUNDWIDTH - 2)) / 11,
1, (BOUNDWIDTH - 2) / 11, 16, "", 2);
} else { /* do an update on a random mouse
* press */
if (canonical &&
(report.xbutton.y > TOPOFFSET) &&
(report.xbutton.x < LEFTOFFSET)) {
beta = 0.5 / (0.6 + 1.0 * (TBOTTOM - report.xbutton.y) / THGT);
bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta));
bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta));
fixthermometer();
} else {
if (1 == algorithm)
startcluster();
else
localupdate();
}
}
break;
default:
break;
}
}
}
return 1;
}
void startcluster()
{ /* flip the old cluster and start a
* new one */
long xedge, yedge;
int row, col, bit;
long currentword;
long index0, index1;
/* flip previous cluster and adjust demons on boundary */
for (row = 0; row < nrows; row++) {
for (col = 0; col < hwords; col++) {
index0 = (demonindex + 2 * (col + row * hwords));
index1 = (index0 + 1);
while (index0 >= ndemons)
index0 -= ndemons; /* faster than using % */
while (index1 >= ndemons)
index1 -= ndemons;
if (0 == unchecked[row][col]) { /* don't bother if unchecked */
unchecked[row][col] = 1;/* reset unchecked */
currentword = cluster[row][col];
/* find cluster edges */
xedge = currentword ^ ((currentword << 1)
| (1 & (cluster[row][nextcol] >> (WORDSIZE - 1))));
yedge = currentword ^ cluster[nextrow][col];
/* update demons and field */
for (bit = 0; bit < DEMONBITS; bit++) {
demon[bit][index0]
= (newdemon[bit][index0] & xedge) |
((~xedge) & demon[bit][index0]);
demon[bit][index1]
= (newdemon[bit][index1] & yedge) |
((~yedge) & demon[bit][index1]);
}
field[row][col] ^= currentword;
} /* end of if checked */
/* do any cooling or heating required */
/*
* This doesn't work very well in the cluster algorithm; the
* demons hold the heat for a while before dumping it, thus the
* thermomenter doesn't read reliably. For major heating and
* cooling it seems better to go to the local algorithm.
*/
if ((!canonical) && (drand48() < 0.05)) {
if (heater) {
demon[0][index0] ^= (RND & RND);
demon[0][index1] ^= (RND & RND);
} else if (cooler) {
demon[0][index0] &= RND;
demon[0][index1] &= RND;
} else { /* rotate the demon's first bit to
* help scramble */
demon[0][index0] = (demon[0][index0] << 1)
| (1 & (demon[0][index0] >> (WORDSIZE - 1)));
demon[0][index1] = (demon[0][index1] << 1)
| (1 & (demon[0][index1] >> (WORDSIZE - 1)));
}
}
/* accumulate counts */
currentword = demon[0][index0];
energycount += bitcount(currentword);
currentword = demon[0][index1];
energycount += bitcount(currentword);
volumecount += WORDSIZE;
} /* end of col loop */
} /* end of row loop */
if (canonical) {
if (heater || cooler) {
beta += (cooler - heater) * .0005;
bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta));
bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta));
fixthermometer();
}
}
/* clear previous cluster */
for (row = 0; row < nrows; row++) {
activerow[row] = 0;
for (col = 0; col < hwords; col++)
cluster[row][col] = 0;
}
/* randomize starting demon location */
demonindex = lrand48() % ndemons;
/* pick random site to seed cluster */
row = lrand48() % nrows;
col = lrand48() % hwords;
clustergrowth = 1;
cluster[row][col] = 1 << (lrand48() % WORDSIZE);
activerow[row] = activerow[previousrow] = activerow[nextrow] = 1;
/* set up bounding box for cluster */
boxbottom = boxtop = row;
/* monitor temperature */
if (canonical) {
refreshdemons(5); /* refresh 1/10th of the demons */
} else if (volumecount >= volume) { /* update beta and
* thermometer */
beta = .5 * log((double) (-1. + volumecount / (0.5 * energycount)));
bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta));
bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta));
fixthermometer();
energycount = 0;
volumecount = 0;
} /* end of temperature monitoring */
/* copy lattice to window */
XPutImage(display, window, gc, spinimage, 0, 0,
LEFTOFFSET, TOPOFFSET, ncols, nrows);
/* display the new cluster */
XPutImage(display, window, gc, clusterimage, 0, 0,
LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows);
return;
} /* end of startcluster */
void growcluster()
/* this new version (July `97) is supposed to be more readable */
{
int rowloop, row, col, nrow, ncol, newstuff;
long edge, currentword, newsites;
clustergrowth = 0;
direction = -direction; /* change sweeping direction */
/* start where last activity was */
if (direction > 0)
row = boxtop;
else
row = boxbottom;
for (rowloop = 0; rowloop < 4 * nrows; rowloop++) {
/* multiple passes seem to help */
if (activerow[row]) {
activerow[row] = 0;
nrow = nextrow;
for (col = 0; col < hwords; col++) {
ncol = nextcol;
currentword = cluster[row][col];
/* look for new stuff if part of cluster is around */
newstuff = (0 !=
(currentword | cluster[nrow][col] | cluster[row][ncol]));
while (newstuff) {
newstuff = 0;
/* make sure necessary demons calculated */
if (unchecked[row][col])
check(row, col);
/* now put things from/to next column */
if (1 & (sadx[row][col])) {
if (1 & (currentword ^
((cluster[row][ncol]) >> (WORDSIZE - 1)))) {
newstuff = 1;
cluster[row][ncol] |= LEFTBIT;
currentword |= 1;
}
}
/* put things from/to next row */
edge = currentword ^ (cluster[nrow][col]);
if ((newsites = edge & sady[row][col])) {
/* single "=" intended */
currentword |= newsites;
cluster[nrow][col] |= newsites;
newstuff = 1;
activerow[nrow] = 1;
}
/* find cluster edge within current word */
edge = (currentword ^ (currentword << 1)) & (~1);
while ((newsites = edge & (sadx[row][col]))) {
/* need to grow */
currentword |= newsites | (newsites >> 1);
edge = (currentword ^ (currentword << 1)) & (~1);
newstuff = 1;
}
if (newstuff) {
cluster[row][col] = currentword;
activerow[row] = activerow[previousrow] = 1;
clustergrowth = 1;
}
} /* end of while newstuff */
} /* end of column loop */
} /* end of if activerow */
if (direction > 0) { /* alternate directions */
if (activerow[row]) /* back up */
boxtop = row = nextrow;
else
row = previousrow;
} else {
if (activerow[row]) /* back up */
boxbottom = row = previousrow;
else
row = nextrow;
}
if (speed) /* speed is not an issue; draw
* progress */
XPutImage(display, window, gc, clusterimage,
0, row, LEFTOFFSET, TOPOFFSET + 28 + nrows + row, ncols, 1);
} /* end of row loop */
/* display new cluster */
XPutImage(display, window, gc, clusterimage, 0, 0,
LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows);
return;
} /* end of growcluster */
void check(row, col)
int row, col;
{ /* calculate sadx and sady when
* unchecked */
/* demons are sad if they cannot flip a bond */
long right, bottom, carry0, carry1;
int bit;
long index0, index1;
right = field[row][col] ^ ((field[row][col] << 1)
| (1 & (field[row][nextcol] >> (WORDSIZE - 1))));
bottom = field[row][col] ^ field[nextrow][col];
if (boundary) { /* correct for antiperiodic
* boundaries */
if (row == nrows - 1)
bottom ^= (-1);
if (col == topcol)
right ^= 1;
}
index0 = (demonindex + 2 * (col + row * hwords)); /* start of xdemons */
index1 = (index0 + 1); /* start of ydemons */
while (index0 >= ndemons)
index0 -= ndemons;
while (index1 >= ndemons)
index1 -= ndemons;
/* copy demons to newdemon */
for (bit = 0; bit < DEMONBITS; bit++) {
newdemon[bit][index0] = demon[bit][index0];
newdemon[bit][index1] = demon[bit][index1];
}
/* subtract 1 from newdemon */
carry0 = carry1 = (-1);
for (bit = 0; bit < DEMONBITS; bit++) {
newdemon[bit][index0] ^= carry0;
newdemon[bit][index1] ^= carry1;
carry0 &= newdemon[bit][index0];
carry1 &= newdemon[bit][index1];
}
/* add two if neighbor antiparallel */
for (bit = 1; bit < DEMONBITS; bit++) {
newdemon[bit][index0] ^= right;
newdemon[bit][index1] ^= bottom;
right &= (~newdemon[bit][index0]);
bottom &= (~newdemon[bit][index1]);
}
carry0 ^= right;
carry1 ^= bottom;
/* make demons sad if overflow */
sadx[row][col] = carry0;
sady[row][col] = carry1;
unchecked[row][col] = 0;
return;
}
void setup()
{
int i, j, count;
/* initialize "bits" used for bitcounts */
for (i = 0; i < 256; i++) {
j = i;
count = j & 1;
while (j)
count += ((j >>= 1) & 1);
bits[i] = count; /* the number of set bits in i */
}
/* make checkermask have alternate bits set */
checkermask = 2;
/*
* Do this 32 times in case we are on a 64 bit machine. The lost
* time otherwise is highly insignificant.
*/
for (i = 0; i < 32; i++)
checkermask |= (checkermask << 2);
return;
}
void repaint()
/* this fixes the window up whenever it is uncovered */
{
int i;
XSetPlaneMask(display, gc, AllPlanes);
drawbutton(quitbutton, 0, 0, BUTTONWIDTH, BUTTONHEIGHT, "quit", 1);
drawbutton(pausebutton, 0, 0, BUTTONWIDTH, BUTTONHEIGHT,
"pause", 1 - 2 * paused);
drawbutton(speedbutton, 0, 0, BOUNDWIDTH, 18, "speed", -1);
drawbutton(speedbutton, 1 + ((10 - speed) * (BOUNDWIDTH - 2)) / 11, 1,
(BOUNDWIDTH - 2) / 11, 16, "", 2);
if (canonical)
drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT,
"canonical", -1);
else
drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT,
"microcanonical", 1);
/* draw the algorithm buttons */
for (i = 0; i < 2; i++)
drawbutton(algobutton, 0, i * 18, 68, 18,
algotext[i], 1 - 2 * (algorithm == i));
/* draw the heater buttons */
for (i = 0; i < 3; i++)
drawbutton(heatwindow, 0, i * HEATERHEIGHT / 3,
HEATERWIDTH, HEATERHEIGHT / 3,
heatertext[i], 1 - 2 * (i == (heater + 2 * cooler)));
/* draw the boundary condition buttons */
for (i = 0; i < 4; i++)
drawbutton(boundwindow, 0, i * BOUNDHEIGHT / 4,
BOUNDWIDTH, BOUNDHEIGHT / 4,
boundtext[i], 1 - 2 * (i == boundary));
/* draw thermometer */
drawbutton(window, TLEFT - 14, TTOP - 14,
TWDTH + 28, TBOTTOM - TTOP + 34, "", 2);
drawbutton(window, TLEFT - 3, TTOP - 2,
TWDTH + 6, TBOTTOM - TTOP + 4, "", 3);
XSetForeground(display, gcpen, black);
XDrawLine(display, window, gcpen, TLEFT - 3, TTOP - 2,
TLEFT + TWDTH + 2, TTOP - 2);
XDrawLine(display, window, gcpen, TLEFT - 3, TTOP - 2,
TLEFT - 3, TBOTTOM + 2);
XDrawLine(display, window, gcpen, TLEFT + TWDTH + 2, TBOTTOM + 2,
TLEFT + TWDTH + 2, TTOP - 2);
XFillArc(display, window, gcpen,
TLEFT - (3 * TWDTH) / 2 - 3, TBOTTOM - TWDTH - 3,
4 * TWDTH + 4, 4 * TWDTH + 4, 0, 64 * 360);
XSetForeground(display, gcpen, white);
XFillArc(display, window, gcpen,
TLEFT - (3 * TWDTH) / 2 - 2, TBOTTOM - TWDTH - 2,
4 * TWDTH, 4 * TWDTH, 0, 64 * 360);
XSetForeground(display, gcpen, darkcolor);
XFillArc(display, window, gcpen,
TLEFT - (3 * TWDTH) / 2, TBOTTOM - TWDTH,
4 * TWDTH - 2, 4 * TWDTH - 2, 0, 64 * 360);
XSetForeground(display, gcpen, spinup);
XFillArc(display, window, gcpen,
TLEFT - TWDTH / 2, TBOTTOM, TWDTH + 2, TWDTH + 2, 10, 64 * 250);
XSetForeground(display, gcpen, black);
for (i = TBOTTOM - 10; i > TTOP; i -= (THGT / 20)) {
XDrawLine(display, window, gcpen, TLEFT - 7, i, TLEFT - 4, i);
XDrawLine(display, window, gcpen, TLEFT + TWDTH + 3, i,
TLEFT + TWDTH + 6, i);
}
i = TBOTTOM - THGT * (-0.6 + 1.0 / log(1. + sqrt(2.)));
XDrawLine(display, window, gcpen, TLEFT - 15, i, TLEFT + TWDTH + 15, i);
XDrawString(display, window, gcpen, TLEFT + TWDTH + 15, i + 4, "T", 1);
XDrawString(display, window, gcpen, TLEFT + TWDTH + 23, i + 8, "c", 1);
fixthermometer();
/* write various strings */
sprintf(stringbuffer, "%d by %d lattice", ncols, nrows);
XDrawString(display, window, gcpen, 150, 92,
stringbuffer, strlen(stringbuffer));
XDrawString(display, window, gcpen, 20, 92, "beta=", 5);
XDrawString(display, window, gcpen, LEFTOFFSET, 114, "spins:", 6);
XDrawString(display, window, gcpen, LEFTOFFSET, 142 + nrows, "changes:", 8);
XDrawString(display, window, gcpen, TLEFT - 8, TBOTTOM + 40, "MJC", 3);
/* draw border and redraw images */
drawbutton(window, LEFTOFFSET - 4, TOPOFFSET - 4,
ncols + 8, nrows + 8, NULL, 4);
drawbutton(window, LEFTOFFSET - 4, TOPOFFSET + 28 + nrows - 4,
ncols + 8, nrows + 8, NULL, 4);
XPutImage(display, window, gc, spinimage, 0, 0,
LEFTOFFSET, TOPOFFSET, ncols, nrows);
XPutImage(display, window, gc, clusterimage, 0, 0,
LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows);
/* this might speed things up a little bit */
XSetPlaneMask(display, gc, spinup ^ spindown);
return;
}
void refreshdemons(int step)
/*
* refresh cluster demons at beta note: Calling this after each
* cluster update gives the canonical algorithm.
*/
/*
* To generate a bit set with probability p: write p as a binary
* fraction. Compare a random bit string with the digits of this
* fraction and find the first place at which they are the same.
* Output the bit that occupies that place. This is done here in
* parallel on WORDSIZE demons at a time.
*/
{
long morebits, randombits, bitprob, currentbit;
double bfactor, prob;
long i;
int bit, j;
bfactor = exp(-beta);
for (bit = 0; bit < DEMONBITS; bit++) {
bfactor *= bfactor; /* gives exp(-2 beta) for the first
* bit, exp(-4 beta) for the next,
* ... */
prob = bfactor / (1. + bfactor); /* prob that current demon
* bit is 1 */
bitprob = ((1 << 16) * prob); /* convert to 16 bit integer */
for (i = 0; i < ndemons; i += step) {
demon[bit][i] = 0;
morebits = (-1);
for (j = 0; j < 16; j++) {
randombits = RND; /* new set of random bits */
currentbit = -((bitprob >> (15 - j)) & 1);
/*
* set demon to current bit if it equals randombits and has
* not been fixed yet
*/
demon[bit][i] |= (currentbit & randombits & morebits);
/* shut off morebits where they are equal */
morebits &= (currentbit ^ randombits);
if (0 == morebits)
break; /* all bits decided */
}
}
}
return;
}
void heatbath()
/*
* refresh local demons at beta
*/
/*
* To generate a bit set with probability p: write p as a binary
* fraction. Compare a random bit string with the digits of this
* fraction and find the first place at which they are the same.
* Output the bit that occupies that place. This is done here in
* parallel on WORDSIZE demons at a time.
*/
{
long morebits, randombits, currentbit;
int j;
demon0 = demon1 = 0;
morebits = (-1);
for (j = 15; (j >= 0) && (morebits != 0); j--) {
randombits = RND; /* new set of random bits */
currentbit = -((bitprob0 >> j) & 1);
/* current bit is a word filled with the j'th bit of bitprob0 */
/*
* set demon to current bit if it equals randombits and has not
* been fixed yet
*/
demon0 |= (currentbit & randombits & morebits);
/* shut off morebits where they are equal */
morebits &= (currentbit ^ randombits);
if (0 == morebits)
break; /* all bits decided */
}
/* now do demon1 */
morebits = (-1);
for (j = 15; (j >= 0) && (morebits != 0); j--) {
randombits = RND;
currentbit = -((bitprob1 >> j) & 1);
demon1 |= (currentbit & randombits & morebits);
morebits &= (currentbit ^ randombits);
if (0 == morebits)
break;
}
return;
}
void localupdate()
/* the local updating routine */
{
long *up, *down; /* pointers to above and below rows */
long temp0, temp1, reject, spins, left, right, top, bottom;
int row = 0, col, rowloop;
checkermask = (~checkermask);
for (rowloop = 0; rowloop < nrows; rowloop++) {
row = (row + 29); /* jump ahead 29 rows; allows for
* energy to move long distances */
while (row >= nrows)
row -= nrows;
/* get neighbors */
up = field[previousrow];
down = field[nextrow];
leftshift(field[row], work0);
rightshift(field[row], work1);
/* to keep speed up, only refresh occasionally */
if ((canonical) && ((row & 3) == 3))
heatbath();
for (col = 0; col < hwords; col++) {
spins = field[row][col];
/* determine antiparallel neighbors */
right = work0[col] ^ spins;
left = work1[col] ^ spins;
top = up[col] ^ spins;
bottom = down[col] ^ spins;
if (0 == row)
switch (boundary) {
case 1:
top ^= (~0);
break;
case 2:
top = (~0) ^ spins;
break;
case 3:
top = spins;
break;
}
else if (nrows - 1 == row)
switch (boundary) {
case 1:
bottom ^= (~0);
break;
case 2:
bottom = (~0) ^ spins;
break;
case 3:
bottom = spins;
break;
}
/* add up top, bottom, and left antiparallel spins */
temp0 = top ^ bottom ^ left;
temp1 = (top & (bottom | left)) | (bottom & left);
/* add in right */
reject = temp0 & temp1 & right;
temp1 ^= (temp0 & right);
temp0 ^= right;
/* add demon0 */
reject ^= (temp0 & temp1 & demon0);
temp1 ^= (temp0 & demon0);
temp0 ^= demon0;
/* add in demon1 */
reject ^= (temp1 & demon1);
temp1 ^= demon1;
/* subtract two */
temp1 = ~temp1;
reject ^= temp1;
/* only change one color sites */
reject |= checkermask;
/* make changes */
demon0 = (demon0 & reject) | (temp0 & (~reject));
demon1 = (demon1 & reject) | (temp1 & (~reject));
field[row][col] = (spins ^ (~reject));
/* put changes into cluster image */
cluster[row][col] &= checkermask;
cluster[row][col] |= (~reject);
energycount += bitcount(demon0 & (~checkermask));
}
/* heat or cool */
if ( !canonical ) {
if (heater)
demon0 |= (RND & RND & RND & RND & RND & RND);
else if (cooler)
demon1 &= (RND | RND | RND | RND);
}
}
/* heat or cool for canonical case */
if (canonical && (heater || cooler)) {
beta += (cooler - heater) * .000005 * nrows;
bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta));
bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta));
fixthermometer();
}
/* monitor temperature every 50 "half" iterations */
iteration++;
if (iteration == 50) { /* update beta and thermometer */
if (!canonical) {
beta = .25 * log((double) (-1. + 50 * volume / (2. * energycount)));
bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta));
bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta));
fixthermometer();
}
iteration = 0;
energycount = 0;
}
/* copy lattice to window */
if (!XPending(display)){
XPutImage(display, window, gc, spinimage, 0, 0,
LEFTOFFSET, TOPOFFSET, ncols, nrows);
XPutImage(display, window, gc, clusterimage, 0, 0,
LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows);
XSync(display,0);
}
return;
}
void fixthermometer()
/* fill the thermometer with fluid */
{
int temp;
sprintf(stringbuffer, "%.3f", beta);
XDrawImageString(display, window, gcpen, LEFTOFFSET, 92,
stringbuffer, strlen(stringbuffer));
temp = THGT * (-0.6 + 0.5 / beta);
if (temp < 0)
temp = 0;
if (temp > THGT)
temp = THGT;
XSetForeground(display, gcpen, darkcolor);
XFillRectangle(display, window, gcpen, TLEFT, TBOTTOM - temp, TWDTH, temp);
XFillRectangle(display, window, gcclear, TLEFT, TTOP, TWDTH, THGT - temp);
XSetForeground(display, gcpen, black);
return;
}
void rightshift(source, destination)
/* shifts a row of cells to the right by one bit */
long *source, *destination;
{
int col;
long carry = 0;
for (col = bottomcol; col != (topcol + nshift); col += nshift) {
destination[col] = (~LEFTBIT) & (source[col] >> 1);
if (carry)
destination[col] |= LEFTBIT;
carry = source[col] & 1;
}
if (carry)
destination[bottomcol] |= LEFTBIT;
switch (boundary) {
case 1:
destination[bottomcol] ^= LEFTBIT;
break;
case 2:
destination[bottomcol] |= LEFTBIT;
break;
case 3:
destination[bottomcol] &= (~LEFTBIT);
break;
}
return;
}
void leftshift(source, destination)
/* shifts a row of cells one bit to the left */
long *source, *destination;
{
int col;
long carry = 0;
for (col = topcol; col != (bottomcol - nshift); col -= nshift) {
destination[col] = (source[col] << 1);
if (carry)
destination[col] |= 1;
carry = source[col] >> (WORDSIZE - 1);
}
if (carry)
destination[topcol] |= 1;
switch (boundary) {
case 1:
destination[topcol] ^= 1;
break;
case 2:
destination[topcol] |= 1;
break;
case 3:
destination[topcol] &= (~1);
break;
}
return;
}
void setalg()
/* take preparatory action if algorithm button hit */
{
int row, col;
energycount = volumecount = iteration = 0;
if (1 == algorithm) { /* reset cluster stuff */
if (boundary > 1) {
setboundary(0);
}
refreshdemons(1);
clustergrowth = 0;
for (col = 0; col < hwords; col++)
for (row = 0; row < nrows; row++) {
cluster[row][col] = 0;
unchecked[row][col] = 1;
}
}
return;
}
void setheater(value)
/* acts on heater buttons */
int value;
{
switch (value) {
case 0:
heater = cooler = 0;
break;
case 1:
heater = 1 - heater;
break;
case 2:
cooler = 1 - cooler;
}
drawbutton(heatwindow, 0, 0,
HEATERWIDTH, HEATERHEIGHT / 3, heatertext[0],
2 * (heater | cooler) - 1);
drawbutton(heatwindow, 0, (HEATERHEIGHT) / 3,
HEATERWIDTH, HEATERHEIGHT / 3, heatertext[1], 1 - 2 * heater);
drawbutton(heatwindow, 0, (HEATERHEIGHT * 2) / 3,
HEATERWIDTH, HEATERHEIGHT / 3, heatertext[2], 1 - 2 * cooler);
return;
}
void setboundary(value)
/* take action if boundary buttons hit */
int value;
{
drawbutton(boundwindow, 0, boundary * BOUNDHEIGHT / 4,
BOUNDWIDTH, BOUNDHEIGHT / 4, boundtext[boundary], 1);
/* don't allow fixed boundaries with cluster algorithm */
if ((algorithm == 1) && (value > 1))
fprintf(stderr,
"fixed boundaries not implemented for cluster algorithm\n");
else
boundary = value;
drawbutton(boundwindow, 0, boundary * BOUNDHEIGHT / 4,
BOUNDWIDTH, BOUNDHEIGHT / 4, boundtext[boundary], -1);
return;
}
void openwindow(argc, argv)
/*
* a lot of this is taken from the basicwin program in the Xlib
* Programming Manual
*/
int argc;
char **argv;
{
char *window_name = "2-d Ising Model";
char *icon_name = "Ising";
Pixmap icon_pixmap;
char *display_name = NULL;
XColor xcolor, colorcell;
Colormap cmap;
/* icon */
#define icon_bitmap_width 16
#define icon_bitmap_height 16
static unsigned char icon_bitmap_bits[] = {
0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
/* open up the display */
if ((display = XOpenDisplay(display_name)) == NULL) {
fprintf(stderr, "%s: cannot connect to X server %s\n",
progname, XDisplayName(display_name));
exit(-1);
}
screen = DefaultScreen(display);
depth = DefaultDepth(display, screen);
cmap = DefaultColormap(display, screen);
/* depth=1; *//* uncomment to test monochrome display */
spindown = darkcolor = black = BlackPixel(display, screen);
spinup = lightcolor = white = WhitePixel(display, screen);
if (depth > 1) { /* color? This is not the right way
* to do it, but ... */
if (XAllocNamedColor(display, cmap, "coral", &colorcell, &xcolor))
darkcolor = colorcell.pixel;
if (XAllocNamedColor(display, cmap, "turquoise", &colorcell, &xcolor))
lightcolor = colorcell.pixel;
if (XAllocNamedColor(display, cmap, "grey", &colorcell, &xcolor))
spinup = colorcell.pixel;
if (XAllocNamedColor(display, cmap, "black", &colorcell, &xcolor))
spindown = colorcell.pixel;
}
/* make the main window */
windowwidth = ncols + 75;
windowheight = 2 * nrows + 164 + 24;
window = XCreateSimpleWindow(display, RootWindow(display, screen),
0, 0, windowwidth, windowheight, 4, BlackPixel(display, screen),
lightcolor);
/* make the icon */
icon_pixmap = XCreateBitmapFromData(display, window,
icon_bitmap_bits, icon_bitmap_width,
icon_bitmap_height);
size_hints.flags = PPosition | PSize | PMinSize;
size_hints.min_width = windowwidth;
size_hints.min_height = 210;
#ifdef X11R3
size_hints.x = 0;
size_hints.y = 0;
size_hints.width = windowwidth;
size_hints.height = windowheight;
XSetStandardProperties(display, window, window_name, icon_name,
icon_pixmap, argv, argc, &size_hints);
#else
{
XWMHints wm_hints;
XClassHint class_hints;
XTextProperty windowName, iconName;
if (XStringListToTextProperty(&window_name, 1, &windowName) == 0) {
fprintf(stderr, "%s: structure allocation for windowName failed.\n"
,progname);
exit(-1);
}
if (XStringListToTextProperty(&icon_name, 1, &iconName) == 0) {
fprintf(stderr, "%s: structure allocation for iconName failed.\n"
,progname);
exit(-1);
}
wm_hints.initial_state = NormalState;
wm_hints.input = True;
wm_hints.icon_pixmap = icon_pixmap;
wm_hints.flags = StateHint | IconPixmapHint | InputHint;
class_hints.res_name = progname;
class_hints.res_class = "Basicwin";
XSetWMProperties(display, window, &windowName, &iconName,
argv, argc, &size_hints, &wm_hints, &class_hints);
}
#endif
/* make a default cursor */
XDefineCursor(display,window, XCreateFontCursor(display,XC_hand2));
/* pick the events to look for */
event_mask = ExposureMask | ButtonPressMask | StructureNotifyMask;
XSelectInput(display, window, event_mask);
/* pick font: 9x15 is supposed to almost always be there */
if ((font = XLoadQueryFont(display, "9x15")) == NULL) {
fprintf(stderr, "%s: Cannot open 9x15 font\n", progname);
exit(-1);
}
font_height = font->ascent + font->descent;
font_width = font->max_bounds.width;
/*
* make graphics contexts: gc for black on white, gcpen for
* various, defaults to black on lightcolor, gclear does the
* reverse, clears thermometer top
*/
gc = XCreateGC(display, window, 0, NULL);
XSetFont(display, gc, font->fid);
XSetForeground(display, gc, spinup);
XSetBackground(display, gc, spindown);
gcpen = XCreateGC(display, window, 0, NULL);
XSetFont(display, gcpen, font->fid);
XSetForeground(display, gcpen, black);
XSetBackground(display, gcpen, lightcolor);
gcclear = XCreateGC(display, window, 0, NULL);
XSetFont(display, gcclear, font->fid);
XSetForeground(display, gcclear, spinup);
XSetBackground(display, gcclear, black);
/* show the window */
XMapWindow(display, window);
return;
}
void allocarray(array, xsize, ysize, datasize)
char ***array;
int xsize, ysize, datasize;
/*
* dynamically allocates a two dimensional (*array)[xize][ysize] with
* datasize bytes per element
*/
{
int i;
if ((*array) != ((char **) NULL)) /* free previous allocation */
free(*array);
/* allocate memory for array of xsize pointers plus data */
*array = (char **) malloc(sizeof(char *) * xsize + datasize * xsize * ysize);
if (*array == (char **) NULL) {
fprintf(stderr, "%s: memory allocation problems\n", progname);
exit(-1);
}
/* initialize the pointers */
for (i = 0; i < xsize; i++)
*((*array) + i) = ((char *) (*array + xsize) + (i * ysize * datasize));
return;
}
void makebuttons()
/*
* creates the buttons and allocates arrays; call after any resizing
*/
{
int row, col;
XDestroySubwindows(display, window);
XSetForeground(display, gcpen, lightcolor);
XFillRectangle(display, window, gcpen, 0, 0, windowwidth, windowheight);
quitbutton = makebutton(4, 4, BUTTONWIDTH, BUTTONHEIGHT);
pausebutton = makebutton(58, 4, BUTTONWIDTH, BUTTONHEIGHT);
boundwindow = makebutton(windowwidth - BOUNDWIDTH - 4, 4,
BOUNDWIDTH, BOUNDHEIGHT);
heatwindow = makebutton(windowwidth - HEATERWIDTH - BOUNDWIDTH - 12, 4,
HEATERWIDTH, HEATERHEIGHT);
algobutton = makebutton(18, 36, 68, 36);
speedbutton = makebutton(windowwidth - BOUNDWIDTH - 4,
BOUNDHEIGHT + 26, BOUNDWIDTH, 18);
canonbutton = makebutton(LEFTOFFSET + 4 + (ncols - font_width * 15) / 2,
windowheight - 26, font_width * 15, BUTTONHEIGHT);
/* reset various lattice dimensions */
hwords = ncols / WORDSIZE;
ncols = hwords * WORDSIZE; /* truncate ncols to a multiple of
* wordsize */
if (0 == (nrows % 29))
nrows--; /* so demon jump of 29 rows works ok */
volume = nrows * ncols;
ndemons = 2 * nrows * hwords;
/* allocate various arrays */
if (NULL != spinimage) {
(*spinimage).data = NULL;
XDestroyImage(spinimage);
}
if (NULL != clusterimage) {
(*clusterimage).data = NULL;
XDestroyImage(clusterimage);
}
allocarray(&field, nrows, hwords, sizeof(long));
allocarray(&cluster, nrows, hwords, sizeof(long));
allocarray(&sadx, nrows, hwords, sizeof(long));
allocarray(&sady, nrows, hwords, sizeof(long));
allocarray(&unchecked, nrows, hwords, sizeof(long));
allocarray(&demon, DEMONBITS, ndemons, sizeof(long));
allocarray(&newdemon, DEMONBITS, ndemons, sizeof(long));
if (NULL != work0)
free(work0);
if (NULL != work1)
free(work1);
if (NULL != activerow)
free(activerow);
work0 = (long *) malloc(sizeof(long) * hwords);
work1 = (long *) malloc(sizeof(long) * hwords);
activerow = (int *) malloc(sizeof(int) * nrows);
if ((NULL == work0) || (NULL == work1)) {
fprintf(stderr, "%s: memory allocation problems\n", progname);
exit(-1);
}
energycount = volumecount = iteration = 0;
/* make image structures */
spinimage = XCreateImage(display, (Visual *) & window, 1, XYBitmap, 0,
(char *) *field, ncols, nrows, 32, 0);
clusterimage = XCreateImage(display, (Visual *) & window, 1, XYBitmap, 0,
(char *) *cluster, ncols, nrows, 32, 0);
if ((NULL == spinimage) || (NULL == clusterimage)) {
fprintf(stderr, "%s: memory allocation problems\n",
progname);
exit(-1);
}
(*spinimage).bitmap_unit = WORDSIZE;
(*clusterimage).bitmap_unit = WORDSIZE;
/* test for byte order on client */
**field = 1;
if (*(char *) (*field)) {
(*spinimage).byte_order = LSBFirst;
(*clusterimage).byte_order = LSBFirst;
/* printf("byte_order=LSBFirst\n"); */
} else {
(*spinimage).byte_order = MSBFirst;
(*clusterimage).byte_order = MSBFirst;
/* printf("byte_order=MSBFirst\n"); */
}
/* find out in what order bits are put on screen */
**field = 0xff;
if (XGetPixel(spinimage, 0, 0)) { /* bits are backward */
nshift = (-1);
topcol = 0;
bottomcol = hwords - 1;
/* printf("bitmap_bit_order=LSBFirst\n"); */
} else {
nshift = 1;
topcol = hwords - 1;
bottomcol = 0;
}
/* set initial state partially random and clear cluster array */
for (col = 0; col < hwords; col++)
for (row = 0; row < nrows; row++) {
field[row][col] = ~(RND * (drand48() > .77));
cluster[row][col] = 0;
unchecked[row][col] = 1;
}
/* equilibrate initial demons at critical temperature */
beta = 0.5 * log(1.0 + sqrt(2.0)); /* critical coupling */
bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta));
bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta));
refreshdemons(1);
heatbath();
return;
}
Window makebutton(xoffset, yoffset, xsize, ysize)
int xoffset, yoffset, xsize, ysize;
/*
* Puts button of specified dimensions on window. Nothing is drawn.
*/
{
Window buttonwindow;
long event_mask;
Cursor cursor;
buttonwindow = XCreateSimpleWindow(display, window, xoffset, yoffset,
xsize, ysize, 0, black, lightcolor);
event_mask = ButtonPressMask | ExposureMask; /* look for mouse-button
* presses */
XSelectInput(display, buttonwindow, event_mask);
/* use a hand cursor for buttons */
cursor=XCreateFontCursor(display,XC_hand2);
XDefineCursor(display,buttonwindow,cursor);
XMapWindow(display, buttonwindow);
return buttonwindow;
}
void drawbutton(buttonwindow, xoffset, yoffset, xsize, ysize, text, state)
Window buttonwindow;
int xoffset, yoffset, xsize, ysize, state;
char *text;
/*
* Draw a button in buttonwindow of specified dimensions with text
* centered. Color of button determined by sign of "state." size of
* border determined by magnitude.
*/
{
int textlength, i, j;
int cdark, clight, cup, cdown;
int cleft, cright, cbutton, ctext;
cup = lightcolor;
cdown = darkcolor;
cdark = black;
clight = white;
if (state < 0) {
cbutton = cdown;
ctext = clight;
cleft = cdark;
cright = clight;
} else {
cbutton = cup;
ctext = cdark;
cleft = clight;
cright = cdark;
}
if (1 == depth)
clight = cleft = cright = black;
j = abs(state);
XSetForeground(display, gcpen, cbutton);
XFillRectangle(display, buttonwindow, gcpen, xoffset + j, yoffset + j,
xsize - 2 * j, ysize - 2 * j);
XSetForeground(display, gcpen, cleft);
XFillRectangle(display, buttonwindow, gcpen, xoffset, yoffset, xsize, j);
XFillRectangle(display, buttonwindow, gcpen, xoffset, yoffset, j, ysize);
XSetForeground(display, gcpen, cright);
for (i = 0; i < j; i++) {
XDrawLine(display, buttonwindow, gcpen,
xoffset + i, yoffset + ysize - i - 1,
xoffset + xsize - i - 1, yoffset + ysize - i - 1);
XDrawLine(display, buttonwindow, gcpen,
xoffset + xsize - i - 1, yoffset + i,
xoffset + xsize - i - 1, yoffset + ysize - i - 1);
}
if (NULL != text) {
textlength = strlen(text);
XSetForeground(display, gcpen, ctext);
XDrawString(display, buttonwindow, gcpen,
xoffset + (xsize - textlength * font_width) / 2,
yoffset + (ysize + font_height) / 2 - 2,
text, textlength);
}
XSetForeground(display, gcpen, darkcolor);
return;
}
Raw
xising.txt
Xising
Michael Creutz
[email protected]
WHAT IT IS:
This is an Xwindow program to illustrate a Monte Carlo simulation of a
simple statistical mechanical system. The Ising model demonstrates
ferromagnetism, displaying a second order phase transition from a
disordered state at high temperature to an ordered state when cool.
The program allows dynamical observation of this system under several
simple Monte Carlo algorithms.
The system is a two dimensional array of spins, each of which is
represented by a bit in a bit map displayed on the screen. Black and
grey represent the two possible states of a given spin. The energy of
the system is entirely determined by nearest neighbor pairs. If two
neighbors are the same, then the energy has one value, and if they
differ, it is higher by two units.
The program source is freely distributable. The latest version is
kept on our WWW server at
"http://thy.phy.bnl.gov/www/xtoys/xtoys.html"
along with some related programs.
COMPILING:
The program should run on anything supporting Xwindows. If you find
an exception, please let me know. To compile, try the xtoys Makefile.
Otherwise, "cc -O -o xising -L/usr/X11R6/lib xising.c -lm -lX11"
usually works. On Sun workstations the X stuff is in an
unconventional place; try "cc -O -I/usr/openwin/include
-L/usr/openwin/lib newxising.c -lm -lX11". Replacing "cc" with "gcc"
may also help. If either of these doesn't work try using the makefile
found at the above web site. If you still have problems, possibly the
X11 includes are not being found. Then you need to compile with a -I
option to where them and possibly change the -lX11 to help the linker.
Remember, as with any Xwindow program, you need to have "xhost" and
"DISPLAY" set up properly to run it.
WHAT THE SCREEN SHOWS:
The display shows the spins in the image labeled "spins." This is
above another bit map labeled "changes," representing the spins being
changed under the current algorithm. The display has a thermometer to
indicate the temperature and various buttons for controlling the
updating. The inverse temperature, beta, is displayed above the
thermometer. The critical value for beta is exactly known to be
0.5*log(1+sqrt(2))=.44068.... The critical temperature is marked on
the thermometer.
One set of three buttons, labeled "run free," "heat," and "cool," are
for adjusting the temperature. Normal running is with the
energy/temperature conserved, but to heat or cool the system toggle
the appropriate button. In a canonical mode (see below) the
temperature can also be adjusted by clicking on the thermometer.
Another set of buttons controls the boundary conditions. The choices
include "periodic," wherein the spin neighbors at an edge are those
on the opposite edge, and "antiperiodic," which is similar but with the
the wrapped neighbors inverted. The other two boundary choices
are constant white or black. The latter two are not permitted when
the cluster algorithm, discussed below, is in operation.
The algorithm button toggles between two alternatives. The program
starts with a "local" updating scheme. This can be toggled into a
variation of the "cluster" algorithm of R.H. Swendsen and J.-S. Wang,
Phys. Rev. Letters 58, 86 (1987), where a large block of spins is
grown and then flipped in unison.
At the bottom of the display is a button that switches between a
canonical and a microcanonical updating. The latter case is described
in some detail in my paper Phys. Rev. Letters 50, 1411 (1983). A set
of "demons" circulates around the lattice trying to flip spins. Each
carries a two bit sack of energy ranging from 0 to 16 units in steps
of 4. Any energy change associated with a spin flip is compensated by
a change in this sack. If the demon's sack cannot accommodate the
change, the flip is rejected. In this mode the temperature is not
fixed, but calculated from the average demon energy. Note the thermal
fluctuations on the thermometer when in a microcanonical mode. The
microcanonical version of the cluster algorithm is discussed in
Phys. Rev. Letters 69, 1002 (1992).
The program attains its speed by updating spins one word at a time
using multispin coding and bit manipulation. The canonical modes are
obtained by allowing the demons to be occasionally refreshed by
visiting a heat bath. When in one of these, a mouse click near the
thermometer will adjust the temperature to the value pointed at by the
cursor. Clicking elsewhere in the main window when the system is
paused with any algorithm does a single step. The remaining buttons,
"quit" and "pause" should be self explanatory, as should the "speed"
slider. The lattice dimensions can be adjusted by resizing the
window.
SOME EXPERIMENTS:
After starting the program, press the heat button and observe how the
domains get small and the acceptance, as shown in the "changes"
display, gets large. Then press the cool button until the
temperature, as displayed in the thermometer, is below the critical
value. Watch the domains grow as the system magnetizes. Note how the
acceptance is largest at the domain boundaries.
In many cases a single domain will grow to dominate the system. If,
however, bands of different phases wrap around the lattice in either a
horizontal or vertical direction, then the system can have a hard time
deciding which phase will dominate and it can remain metastable for a
long time.
Switching the boundary conditions to antiperiodic forces the system to
have at least one domain wall. Switching between black and white
edges allows one to create large included domains which gradually
shrink away.
Returning the system to near the critical temperature, switch to
the cluster algorithm. In this case a few iterations quickly give
independent configurations. Heat the system and observe how
the typical cluster size becomes quite small. Cooling the system
below the critical temperature gives single clusters covering most
of the system, which then flashes between dominantly spin up and
down.
To illustrate the power of the cluster approach, use the local
algorithm to heat the system to a high temperature and then rapidly
quench it to somewhat below the critical value. Before the local
approach has had time to have the smaller domains dissolve in the
dominant one, change to the cluster approach. Note how quickly the
cluster sweeps anneal out the included domains.
Raw
xpotts.c
/* Xpotts, a Potts model simulator for X windows.
Michael Creutz
[email protected]
compiling: cc -O -L/usr/X11R6/lib xpotts.c -lm -lX11
latest version kept at:
http://thy.phy.bnl.gov/www/xtoys/xtoys.html
version 1.5
August 2006
minor modifications to eliminate gcc -Wall warnings, March 1999
*/
# include <X11/Xlib.h>
# include <X11/Xutil.h>
# include <X11/Xos.h>
# include <X11/Xatom.h>
# include <X11/cursorfont.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <malloc.h>
# include <math.h>
/* dimensions for control boxes */
# define HEATERHEIGHT 80
# define HEATERWIDTH 96
# define MHEIGHT 80
# define MWIDTH 96
# define QWIDTH 162
# define PLAYLEFT 20
# define PLAYTOP 80
# define XHEIGHT 40
# define BUTTONWIDTH 66
# define SRHEIGHT (40)
# define BBWIDTH 96
# define INFOLEFT (block*ncols+56)
typedef unsigned char spin;
spin *field=NULL; /* pointer to system */
int nrows,ncols,volume; /* lattice dimensions (plus two for boundaries) */
int q=8; /* the `q' state potts model, default to 8 */
int block, /* size of cells displayed on the screen */
blockshift; /* log_2(block) */
int heatervalue=2, magvalue=3; /* start with cooling and free magnitization */
int demon=0,mdemon=0; /* the energy and magnetism demons, respectively */
# define DEMONMASK 0xf /* how big the demon can get */
int iteration=0,energycount=0,mcount=0,volumecount=0;
long mrand48(),lrand48(); /* generates a random word, lrand is non-negative */
int paused=0; /* is the system paused? */
char stringbuffer[256]; /* generally useful */
void drawbutton(),openwindow(),makebuttons(),update(),repaint(),cleanup(),
showpic(),setheater(),setmag(),setq(),savepic(),loadpic();
int energy();
/* things for regulating the updating speed */
int speed=0; /* updating delay proportional to speed */
int delay=0; /* counter for implementing speed control */
struct timeval timeout; /* timer for speed control */
/* various window stuff */
Display *display;
int screen;
static char *progname;
Window window,quitbutton,pausebutton,srbutton,heatwindow,magwindow,
speedbutton,qbutton,blockbutton,makebutton();
XColor xcolor,colorcell;
Colormap cmap;
GC gc;
int windowwidth,windowheight;
XFontStruct *font=NULL;
int font_height,font_width;
XSizeHints size_hints;
int darkcolor,lightcolor,black,white;
XImage *spinimage=NULL;
long translate[256]; /* for converting colors */
/* text for some of the buttons */
char *srtext[2]={"save","restore"};
char *heatertext[4]={"conserve E","heat","cool","free E"};
char *magtext[4]={"conserve M","increase M","decrease M","free M"};
char *qtext[3]={"<","q",">"};
int main(argc,argv)
int argc;
char **argv;
{unsigned int width, height;
int i,y;
XEvent report;
progname=argv[0];
/* initial lattice size */
block=1;
nrows=192/block;
ncols=192/block;
ncols&=(~(sizeof(long)-1));
if (!((nrows-2)%29)) nrows-=1; /* because demons hop 29 rows */
if (!((ncols-2)%29)) ncols-=sizeof(long);
volume=nrows*ncols;
openwindow(argc,argv);
makebuttons();
/* loop forever, looking for events */
while(1)
{if ((0==XPending(display))&&(0==paused))
{if (delay) /* don't update yet */
{delay--;
/* this use of select() seems a kludge to me;
why can't usleep() be more standard? */
timeout.tv_sec=0;
timeout.tv_usec=100000; /* .1 sec per delay unit */
select(0,NULL,NULL,NULL,&timeout);
}
else
{delay=speed;
update();
}
}
else
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if ((report.xexpose.window)!=window) break; /* cuts down flashing,
but you might remove this line if things aren't being redrawn */
if (report.xexpose.count!=0) break; /* more in queue, wait for them */
repaint();
break;
case ConfigureNotify:
width=report.xconfigure.width;
height=report.xconfigure.height;
if ((width<size_hints.min_width)||(height<size_hints.min_height))
{fprintf(stderr,"%s: window too small to proceed.\n",progname);
cleanup();
}
else if ((width!=windowwidth)||(height!=windowheight))
{windowwidth=width;
windowheight=height;
nrows=(windowheight-PLAYTOP-20)/block;
ncols=(windowwidth-82-190)/block;
ncols&=(~(sizeof(long)-1));
if (!((nrows-2)%29)) nrows-=2;
if (!((ncols-2)%29)) ncols-=sizeof(long);
volume=nrows*ncols;
makebuttons();
showpic();
}
break;
case ButtonPress:
if (report.xbutton.window==quitbutton)
cleanup();
else if (report.xbutton.window==pausebutton)
{paused=1-paused;
drawbutton(pausebutton,0,0,BUTTONWIDTH,18,"pause",1-2*paused);
}
else if (report.xbutton.window==heatwindow)
{setheater((report.xbutton.y*4)/HEATERHEIGHT);
}
else if (report.xbutton.window==magwindow)
{setmag((report.xbutton.y*4)/MHEIGHT);
}
else if (report.xbutton.window==qbutton)
{setq((report.xbutton.x*3)/QWIDTH);
}
else if (report.xbutton.window==srbutton)
{y=(report.xbutton.y*2)/SRHEIGHT;
drawbutton(srbutton,0,y*(SRHEIGHT/2),BUTTONWIDTH,SRHEIGHT/2,
srtext[y],-1);
if (0==y)
savepic(field,ncols,nrows);
else if (1==y) /* load in gif image */
{loadpic(field,ncols,nrows);
for (i=0;i<volume;i++)
while (((int) field[i])>=q) field[i]-=q; /* mod q the image */
showpic();
}
drawbutton(srbutton,0,y*(SRHEIGHT/2),BUTTONWIDTH,SRHEIGHT/2,
srtext[y],1);
}
else if (report.xbutton.window==blockbutton)
{i=(1<<(4*report.xbutton.x/HEATERWIDTH));
if (i!=block) /* reset for new block size */
{block=i;
nrows=(windowheight-PLAYTOP-20)/block;
ncols=(windowwidth-82-190)/block;
ncols&=(~(sizeof(long)-1));
if (!((nrows-2)%29)) nrows-=2;
if (!((ncols-2)%29)) ncols-=sizeof(long);
volume=nrows*ncols;
makebuttons();
showpic();
}
}
else if (report.xbutton.window==speedbutton)
{ /* reset speed */
speed=10-(11*report.xbutton.x)/HEATERWIDTH;
if (speed<0) speed=0;
if (speed>10) speed=10;
delay=speed;
drawbutton(speedbutton,0,0,HEATERWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(HEATERWIDTH-2))/11,1,
(HEATERWIDTH-2)/11,16,"",2);
}
else
update(); /* do a sweep when mouse clicked */
break;
default:
break;
} /* end of switch */
} /* end of if XPending */
} /* end of while(1) */
} /* end of main */
void update()
{int newdemon,oldenergy,newenergy,direction;
int rowloop,colloop,newmdemon;
double beta,hfield,E,M;
int i,row,col;
spin i1,i2,i3,i4,trial,oldspin,*ptr,*ptr1;
row=col=1;
for (rowloop=2;rowloop<nrows;rowloop++)
{row+=29;
while (row>=nrows-1) row-=(nrows-2);
direction=1+lrand48()%(q-1);
if (heatervalue) /* if heater on, adjust demon */
if (8&lrand48()) /* half the time */
{if (2==heatervalue) demon=0;
if ((1==heatervalue)&&(demon<DEMONMASK)) demon++;
if (3==heatervalue) demon=DEMONMASK>>1;
}
if (magvalue) /* adjust magnetic demon */
if (0==(60&lrand48())) /* one in 16 chance since 60 has 4 bits set */
{if (2==magvalue) mdemon=0;
if (1==magvalue) mdemon=DEMONMASK;
if (3==magvalue) mdemon=DEMONMASK>>1;
}
ptr1=field+ncols*row;
for (colloop=2;colloop<ncols;colloop++)
{col+=29;
while (col>=ncols-1) col-=(ncols-2);
ptr=ptr1+col;
oldspin=(*ptr);
i1=(*(ptr+1)); /* right neighbor */
i2=(*(ptr-1)); /* left neighbor */
i3=(*(ptr+ncols)); /* top neighbor */
i4=(*(ptr-ncols)); /* bottom neighbor */
trial=oldspin+direction;
if (((int) trial)>=q) trial-=q;
oldenergy=(oldspin==i1)
+(oldspin==i2)
+(oldspin==i3)
+(oldspin==i4);
newenergy=(trial==i1)
+(trial==i2)
+(trial==i3)
+(trial==i4);
newdemon=demon-oldenergy+newenergy;
newmdemon=mdemon-(0==trial)+(0==oldspin);
if (0==((newdemon|newmdemon)&(~DEMONMASK))) /* accept changes */
{demon=newdemon;
mdemon=newmdemon;
*ptr=trial;
direction=q-direction;
/* fix vertical boundary */
if (1==col)
*(ptr+ncols-2)=(*ptr);
else if ((ncols-2)==col)
*(ptr-ncols+2)=(*ptr);
}
volumecount++;
energycount+=(2&demon); /* use second bit so works for q=2 */
mcount+=(1&mdemon);
}
/* fix horizontal boundary */
if (1==row)
for (i=0;i<ncols;i++)
field[volume-ncols+i]=field[i+ncols];
else if ((nrows-2)==row)
for (i=0;i<ncols;i++)
field[i]=field[volume-2*ncols+i];
}
iteration++;
if (20==iteration) /* display information */
{if ((2*volumecount>energycount)&&(0!=energycount))
{beta=.5*log((double)(2*volumecount/(1.0*energycount)-1.0));
sprintf(stringbuffer," beta=%6.3f ",beta);
}
else
sprintf(stringbuffer," beta divergent");
XDrawImageString(display,window,gc,INFOLEFT+6,150
,stringbuffer,strlen(stringbuffer));
if ((volumecount>mcount)&&(0!=mcount))
{hfield=(-log((double) (volumecount/(1.0*mcount)-1.0)));
sprintf(stringbuffer," h =%6.3f ",hfield);
}
else
sprintf(stringbuffer," h divergent ");
XDrawImageString(display,window,gc,INFOLEFT+6,168
,stringbuffer,strlen(stringbuffer));
iteration=volumecount=energycount=mcount=0;
energy(&E,&M);
sprintf(stringbuffer," E =%6.3f ",E);
XDrawImageString(display,window,gc,INFOLEFT+6,186
,stringbuffer,strlen(stringbuffer));
sprintf(stringbuffer," M =%6.3f ",M);
XDrawImageString(display,window,gc,INFOLEFT+6,202
,stringbuffer,strlen(stringbuffer));
}
showpic();
return;
}
int energy(e,m)
double *e,*m;
{int i,row,col;
int etot=0,mtot=0;
for (row=1;row<nrows-1;row++)
for (col=1;col<ncols-1;col++)
{i=col+ncols*row;
etot+= (field[i]==field[i+1])
+(field[i]==field[i+ncols]);
mtot+=(0==field[i]);
}
*e=(1.0-etot/(2.0*(nrows-2)*(ncols-2)))/(1-1.0/q);
*m=(mtot/(1.0*(nrows-2)*(ncols-2))-1.0/q)/(1-1.0/q);
return etot;
}
void showpic() /* display the field */
{int row,col,i1,i2,color,j,j1,j2,blocktop=block;
char *picture=(*spinimage).data;
if (block>4) blocktop=block-1;
if (8==(*spinimage).depth)
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[row*ncols+col]];
j=block*(col+block*ncols*row);
if (color!=picture[j])
for (i1=0;i1<blocktop;i1++)
{j1=i1*block*ncols+j;
for (i2=0;i2<blocktop;i2++)
picture[j1+i2]=color;
}
}
else
{for (j=0;j<volume;j++)
picture[j]=translate[field[j]];
}
}
else /* depth is not 8, use xputpixel (this is really ugly) */
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[row*ncols+col]];
if (color!=XGetPixel(spinimage,j1=block*col,j2=block*row))
for (i2=0;i2<blocktop;i2++)
for (i1=0;i1<blocktop;i1++)
XPutPixel(spinimage,j1+i1,j2+i2,color);
}
else
for (row=0;row<nrows;row++)
{j1=row*ncols;
for (col=0;col<ncols;col++)
XPutPixel(spinimage,col,row,translate[field[j1+col]]);
}
}
XPutImage(display,window,gc,spinimage,0,0,PLAYLEFT,PLAYTOP,
block*ncols,block*nrows);
return;
}
void repaint()
/* this fixes the window up whenever it is uncovered */
{int i;
XSetForeground(display,gc,translate[14]);
XFillRectangle(display,window,gc,0,0,windowwidth,windowheight);
XSetForeground(display,gc,darkcolor);
drawbutton(pausebutton,0,0,BUTTONWIDTH,18,"pause",1-2*paused);
drawbutton(quitbutton,0,0,BUTTONWIDTH,18,"quit",1);
drawbutton(window,PLAYLEFT-4,PLAYTOP-4,block*ncols+8,block*nrows+8,NULL,4);
drawbutton(window,INFOLEFT,132,208,120,NULL,2);
drawbutton(speedbutton,0,0,HEATERWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(HEATERWIDTH-2))/11,1,
(HEATERWIDTH-2)/11,16,"",2);
for (i=0;i<2;i++)
drawbutton(srbutton,0,i*(SRHEIGHT/2),BUTTONWIDTH,SRHEIGHT/2,srtext[i],1);
/* fix the block buttons */
for (i=1;i<=4;i++)
{sprintf(stringbuffer,"%d",(1<<i)/2);
drawbutton(blockbutton,(i-1)*BBWIDTH/4,0,BBWIDTH/4,18,
stringbuffer,1-2*(i==(blockshift+1)));
}
/* fix the q button */
for (i=0;i<3;i+=2)
drawbutton(qbutton,i*(QWIDTH/3),0,QWIDTH/3,18,qtext[i],1);
setq(1);
/* draw the heater buttons */
for (i=0;i<4;i++)
drawbutton(heatwindow,0,i*HEATERHEIGHT/4,HEATERWIDTH,HEATERHEIGHT/4,
heatertext[i],1-2*(i==heatervalue));
/* draw the magnitization buttons */
for (i=0;i<4;i++)
drawbutton(magwindow,0,i*MHEIGHT/4,MWIDTH,MHEIGHT/4,
magtext[i],1-2*(i==magvalue));
/* write various strings */
sprintf(stringbuffer,"%d by %d lattice ",ncols-2,nrows-2);
XDrawImageString(display,window,gc,INFOLEFT+10,242
,stringbuffer,strlen(stringbuffer));
XDrawString(display,window,gc,
windowwidth-3*HEATERWIDTH/2-4,HEATERHEIGHT+42,
"cell size",9);
XDrawString(display,window,gc,windowwidth-50,windowheight-5,"MJC",3);
showpic();
return;
}
void setheater(value)
/* acts when heater buttons hit */
int value;
{value&=3;
drawbutton(heatwindow,0,(HEATERHEIGHT*heatervalue)/4,
HEATERWIDTH,HEATERHEIGHT/4,heatertext[heatervalue],1);
heatervalue=value;
drawbutton(heatwindow,0,(HEATERHEIGHT*heatervalue)/4,
HEATERWIDTH,HEATERHEIGHT/4,heatertext[heatervalue],-1);
return;
}
void setmag(value)
/* acts when magnitization buttons hit */
int value;
{value&=3;
drawbutton(magwindow,0,(MHEIGHT*magvalue)/4,MWIDTH,MHEIGHT/4,
magtext[magvalue],1);
magvalue=value;
drawbutton(magwindow,0,(MHEIGHT*magvalue)/4,MWIDTH,MHEIGHT/4,
magtext[magvalue],-1);
return;
}
void setq(value)
int value;
{int i;
if ((0==value)&&(q>2))
{q--;
for (i=0;i<volume;i++)
if (((int) field[i])>=q) field[i]--;
showpic();
}
if ((2==value)&&(q<256)) q++;
sprintf(stringbuffer,"q=%d",q);
drawbutton(qbutton,QWIDTH/3,0,QWIDTH/3,18,stringbuffer,1);
sprintf(stringbuffer,"(critical beta=%6.3f)",log((double)(1+sqrt(1.*q))));
XDrawImageString(display,window,gc,INFOLEFT+6,224
,stringbuffer,strlen(stringbuffer));
return;
}
void openwindow(argc,argv)
/* a lot of this is taken from basicwin in the Xlib Programming Manual */
int argc;
char **argv;
{char *window_name="Potts";
char *icon_name="Potts";
long event_mask;
Pixmap icon_pixmap;
char *display_name=NULL;
int i;
# define icon_bitmap_width 16
# define icon_bitmap_height 16
static unsigned char icon_bitmap_bits[] = {
0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
/* open up the display */
if ((display=XOpenDisplay(display_name))==NULL)
{fprintf(stderr,"%s: cannot connect to X server %s\n",
progname,XDisplayName(display_name));
exit(-1);
}
screen=DefaultScreen(display);
cmap=DefaultColormap(display,screen);
darkcolor=black=BlackPixel(display,screen);
lightcolor=white=WhitePixel(display,screen);
translate[0]=white;
translate[1]=black;
translate[2]=lightcolor;
translate[3]=darkcolor;
for(i=4;i<256;i++)
translate[i]=translate[i%4];
if (XAllocNamedColor(display,cmap,"salmon",&colorcell,&xcolor))
darkcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"wheat",&colorcell,&xcolor))
lightcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"red",&colorcell,&xcolor))
translate[0]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"blue",&colorcell,&xcolor))
translate[1]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"green",&colorcell,&xcolor))
translate[2]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"cyan",&colorcell,&xcolor))
translate[3]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"orange",&colorcell,&xcolor))
translate[4]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"purple",&colorcell,&xcolor))
translate[5]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"yellow",&colorcell,&xcolor))
translate[6]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"pink",&colorcell,&xcolor))
translate[7]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"brown",&colorcell,&xcolor))
translate[8]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"grey",&colorcell,&xcolor))
translate[9]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"turquoise",&colorcell,&xcolor))
translate[10]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"gold",&colorcell,&xcolor))
translate[11]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"magenta",&colorcell,&xcolor))
translate[12]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"navy",&colorcell,&xcolor))
translate[13]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"tan",&colorcell,&xcolor))
translate[14]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"violet",&colorcell,&xcolor))
translate[15]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"white",&colorcell,&xcolor))
translate[16]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"black",&colorcell,&xcolor))
translate[17]=colorcell.pixel;
/* feel free to type in more colors, I got bored */
for(i=18;i<256;i++)
translate[i]=translate[i%18];
/* make the main window */
windowwidth=block*ncols+82+190;
windowheight=PLAYTOP+block*nrows+20;
window=XCreateSimpleWindow(display,RootWindow(display,screen),
0,0,windowwidth,windowheight,4,black,lightcolor);
/* make the icon */
icon_pixmap=XCreateBitmapFromData(display,window,
icon_bitmap_bits,icon_bitmap_width,icon_bitmap_height);
size_hints.flags=PPosition | PSize | PMinSize;
size_hints.min_width=windowwidth;
size_hints.min_height=windowheight;
#ifdef X11R3
size_hints.x=x;
size_hints.y=y;
size_hints.width=windowwidth;
size_hints.height=windowheight;
XSetStandardProperties(display,window,window_name,icon_name,
icon_pixmap,argv,argc,&size_hints);
#else
{XWMHints wm_hints;
XClassHint class_hints;
XTextProperty windowName, iconName;
if (XStringListToTextProperty(&window_name,1,&windowName)==0)
{fprintf(stderr,"%s: structure allocation for windowName failed.\n"
,progname);
exit(-1);
}
if (XStringListToTextProperty(&icon_name,1,&iconName)==0)
{fprintf(stderr,"%s: structure allocation for iconName failed.\n"
,progname);
exit(-1);
}
wm_hints.initial_state=NormalState;
wm_hints.input=True;
wm_hints.icon_pixmap=icon_pixmap;
wm_hints.flags=StateHint|IconPixmapHint|InputHint;
class_hints.res_name=progname;
class_hints.res_class="Basicwin";
XSetWMProperties(display,window,&windowName,&iconName,
argv,argc,&size_hints,&wm_hints,&class_hints);
}
#endif
/* pick the events to look for */
event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
XSelectInput(display,window,event_mask);
/* pick font: 9x15 is supposed to almost always be there */
if ((font=XLoadQueryFont(display,"9x15"))==NULL)
if ((font=XLoadQueryFont(display,"fixed"))==NULL)
{fprintf(stderr,"%s: Sorry, having font problems.\n",progname);
exit(-1);
}
font_height=font->ascent+font->descent;
font_width=font->max_bounds.width;
/* make graphics context: */
gc=XCreateGC(display,window,0,NULL);
XSetFont(display,gc,font->fid);
XSetForeground(display,gc,darkcolor);
XSetBackground(display,gc,lightcolor);
/* show the window */
XMapWindow(display,window);
return;
}
void makebuttons()
{int i;
XEvent report;
Cursor cursor;
/* first destroy any old buttons */
XDestroySubwindows(display,window);
/* now make the new buttons */
quitbutton =makebutton(12,4,BUTTONWIDTH,18);
pausebutton =makebutton(QWIDTH-BUTTONWIDTH,4,BUTTONWIDTH,18);
srbutton =makebutton(QWIDTH+14,4,BUTTONWIDTH,SRHEIGHT);
heatwindow =makebutton(windowwidth-HEATERWIDTH-6,4,HEATERWIDTH,HEATERHEIGHT);
magwindow =makebutton(windowwidth-MWIDTH-HEATERWIDTH-14,4,MWIDTH,MHEIGHT);
qbutton =makebutton(6,26,QWIDTH,18);
blockbutton =makebutton(windowwidth-3*HEATERWIDTH/2-10,HEATERHEIGHT+10,
BBWIDTH,18);
speedbutton=makebutton(40,48,HEATERWIDTH,18);
/* wait for window to be displayed befor proceeding */
i=1; /* a flag */
while (i)
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if (report.xexpose.window!=window) i=0;
default:
break;
}
}
/* destroy old image structure */
if (NULL!=spinimage)
{XDestroyImage(spinimage);
spinimage=NULL;
}
/* clear out the play area and create the new image */
XFillRectangle(display,window,gc,0,0,windowwidth,windowheight);
spinimage=XGetImage((Display *) display, (Drawable) window,
PLAYLEFT,PLAYTOP,block*ncols,block*nrows,
AllPlanes,ZPixmap);
if (NULL==spinimage)
{fprintf(stderr,"trouble creating image structure\n");
exit(-1);
}
/* make special cursors to be cute */
cursor=XCreateFontCursor(display,XC_hand2);
XDefineCursor(display,pausebutton,cursor);
XDefineCursor(display,quitbutton,cursor);
XDefineCursor(display,srbutton,cursor);
XDefineCursor(display,blockbutton,cursor);
XDefineCursor(display,heatwindow,cursor);
XDefineCursor(display,magwindow,cursor);
XDefineCursor(display,qbutton,cursor);
XDefineCursor(display,blockbutton,cursor);
XDefineCursor(display,speedbutton,cursor);
/* reallocate field array */
if (NULL!=field) free((char *) field);
if (NULL==( field= (spin *) malloc(sizeof(spin)*volume)))
{fprintf(stderr,"allocation problems\n");
cleanup();
}
/* set initial state random */
for (i=0;i<volume;i++)
field[i]=lrand48()%q;
/* set blockshift to log_2 of block */
blockshift=0;
i=block;
while(i>>=1)
blockshift++;
/* draw everything */
repaint();
return;
}
void cleanup()
{XUnloadFont(display,font->fid);
XFreeGC(display,gc);
XCloseDisplay(display);
XDestroyImage(spinimage);
if (NULL!=field) free((char *) field);
exit(1);
}
Window makebutton(xoffset,yoffset,xsize,ysize)
int xoffset,yoffset,xsize,ysize;
/* Puts button of specified dimensions on window. Nothing is drawn. */
{Window buttonwindow;
long event_mask;
buttonwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
xsize,ysize,0,black,lightcolor);
event_mask=ButtonPressMask|ExposureMask; /* look for mouse-button presses */
XSelectInput(display,buttonwindow,event_mask);
XMapWindow(display,buttonwindow);
return buttonwindow;
}
void drawbutton(buttonwindow,xoffset,yoffset,xsize,ysize,text,state)
Window buttonwindow;
int xoffset,yoffset,xsize,ysize,state;
char * text;
/* Draw a button in buttonwindow of specified dimensions with text
centered. Color of button determined by sign of "state."
size of border determined by magnitude. */
{int textlength,i,j;
int cdark,clight,cup,cdown;
int cleft,cright,cbutton,ctext;
cup=lightcolor;
cdown=darkcolor;
cdark=black;
clight=white;
if (state<0)
{cbutton=cdown;
ctext=clight;
cleft=cdark;
cright=clight;
}
else
{cbutton=cup;
ctext=cdark;
cleft=clight;
cright=cdark;
}
j=abs(state);
XSetForeground(display,gc,cbutton);
XFillRectangle(display,buttonwindow,gc,xoffset+j,yoffset+j,
xsize-2*j,ysize-2*j);
XSetForeground(display,gc,cleft);
XFillRectangle(display,buttonwindow,gc,xoffset,yoffset,xsize,j);
XFillRectangle(display,buttonwindow,gc,xoffset,yoffset,j,ysize);
XSetForeground(display,gc,cright);
for (i=0;i<j;i++)
{XDrawLine(display,buttonwindow,gc,
xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
XDrawLine(display,buttonwindow,gc,
xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
}
if (NULL!=text)
{textlength=strlen(text);
XSetForeground(display,gc,ctext);
XDrawString(display,buttonwindow,gc,
xoffset+(xsize-textlength*font_width)/2,
yoffset+(ysize+font_height)/2-2,
text,textlength);
}
XSetForeground(display,gc,cdark);
return;
}
/* from here on is the stuff for saving and restoring from a gif file */
/* for info on how gif works, see ftp://network.ucsd.edu/graphics/GIF.shar.Z */
char *picturename="xpotts.gif";
void compress(),decompress(),inittable();
int findcode();
void loadpic(data,xsize,ysize) /* load in a GIF image */
unsigned char *data; /* where the output data starts */
int xsize,ysize; /* output picture dimensions */
/* load in a gif image */
{int i,j,filesize,gwidth,gheight,gvolume;
unsigned char *ptr, *ptr1, *rawgif;
int colorbits,codesize;
FILE *infile;
if (NULL==(infile=fopen(picturename,"r")))
{fprintf(stderr,"couldn't open input file\n");
return;
}
/* find the file size */
fseek(infile, 0L, 2);
filesize = ftell(infile);
fseek(infile, 0L, 0);
/* make a place in memory for the file */
if (NULL==(rawgif= (unsigned char *) malloc(filesize) ))
{fprintf(stderr, "not enough memory to read gif file\n");
return;
}
ptr=rawgif;
/* read in the file */
if (fread(ptr, filesize, 1, infile) != 1)
{fprintf(stderr, "read failed\n");
free((char *) rawgif);
return;
}
fclose(infile);
/* check for GIF signature */
if (strncmp((char *) ptr,"GIF87a", 6))
{fprintf(stderr, "not a GIF87a file\n");
free((char *) rawgif);
return;
}
ptr+=6;
ptr+=4; /* skip over screen size */
colorbits=1+((*ptr)&0xf); /* how many bits of color */
ptr+=2; /* skip over background */
if (*ptr++) /* should be zero */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=(3*(1<<colorbits)); /* skip over colormap */
if (','!=(*ptr++)) /* image separator */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=4; /* skip over image offset */
gwidth=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
gheight=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
if (0x80&(*ptr++)) /* skip over local color map */
ptr+=(3*(1<<colorbits));
/* should now be at start of data */
codesize=(*ptr++);
/* make a place for the decompressed file */
gvolume=gwidth*gheight;
ptr1=(unsigned char *) malloc(gvolume);
decompress(codesize,ptr,ptr1,gvolume);
free((char *) rawgif);
/* map picture into data, allowing for different dimensions */
for (j=0;j<ysize;j++)
{if (j>=gheight) break;
for (i=0;i<xsize;i++)
{if (i>=gwidth) break;
data[i+j*xsize]=ptr1[i+j*gwidth];
}
}
free((char *) ptr1);
return;
}
void savepic(data,xsize,ysize) /* save the field as a GIF image */
unsigned char *data; /* where the input data starts */
int xsize,ysize; /* picture dimensions */
{int i;
int colorbits=5,codesize=5; /* assume a 32 color image */
FILE *outfile;
if (NULL==(outfile=fopen(picturename,"w")))
{fprintf(stderr,"couldn't open output file\n");
return;
}
/* GIF signature */
fwrite("GIF87a",6,1,outfile);
/* screen descriptor */
stringbuffer[0]=xsize&0xff; /* screen width */
stringbuffer[1]=(xsize>>8)&0xff;
stringbuffer[2]=ysize&0xff; /* screen height */
stringbuffer[3]=(ysize>>8)&0xff;
stringbuffer[4]=(0x80) /* M=1; global color map follows */
|((colorbits-1)<<4) /* -1+ bits of color reslution */
|(colorbits-1); /* -1+bits per pixel in image */
stringbuffer[5]=0; /* background color */
stringbuffer[6]=0; /* should be zero */
fwrite(stringbuffer,7,1,outfile);
/* global color map */
for (i=0;i<(1<<colorbits);i++)
{colorcell.pixel=translate[i];
XQueryColor(display,cmap,&colorcell);
fputc(colorcell.red,outfile);
fputc(colorcell.green,outfile);
fputc(colorcell.blue,outfile);
}
/* image descriptor */
stringbuffer[0]=','; /* image descriptor separator */
stringbuffer[1]=0; /* image offset */
stringbuffer[2]=0;
stringbuffer[3]=0;
stringbuffer[4]=0;
stringbuffer[5]=xsize&0xff; /* image width */
stringbuffer[6]=(xsize>>8)&0xff;
stringbuffer[7]=ysize&0xff; /* image height */
stringbuffer[8]=(ysize>>8)&0xff;
stringbuffer[9]=0; /* use global color map, no interlace */
fwrite(stringbuffer,10,1,outfile);
/* start of image data */
fputc(codesize,outfile);
compress(codesize,data,outfile,volume);
/* gif terminator */
fputc(';',outfile);
fclose(outfile);
return;
}
/* LZW compression */
/* hash function assumes TABLELENGTH is a power of 2 */
# define TABLELENGTH (1<<13)
char **addresses=NULL; /* where to find the string */
int *codes=NULL, /* the code value */
*linktonext=NULL, /* the next index in the hash chain */
*lengths=NULL, /* the length of the coded string */
*codeindex=NULL; /* the index for a given code */
int nextcode; /* the next unused code */
/* hashit is supposed to give a unique fairly random number in the table for
each length a and string b */
# define hashit(a,b) (51*a+53*(57*b[0]+59*(61*b[a-1]+b[a>>1])))&(TABLELENGTH-1)
void compress(initcodesize,ptr,outfile,size)
int initcodesize; /* the initial compression bits */
char * ptr; /* where the data comes from */
FILE * outfile; /* where the output goes */
int size; /* how much data */
{int currentcode=0,prefixcode=0,codesize,maxbits=12,maxcode;
int clearcode,eoicode,currentplace=0,length,blocksize=0,bitoffset;
int findcode();
unsigned long outputword;
unsigned char blockbuffer[256]; /* to hold data blocks before writing */
/* allocate space for hash tables */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(addresses=(char **) malloc(sizeof(char *)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* set up initial code table */
inittable(initcodesize);
clearcode=(1<<initcodesize);
eoicode=clearcode+1;
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
/* start with a clear code */
outputword=clearcode;
bitoffset=codesize;
/* now do the compressing */
while (currentplace<size)
{ /* check if codesize needs increasing */
if (nextcode>maxcode)
{codesize++;
maxcode=1<<codesize;
/* if too big, then reset compressor */
if (codesize>maxbits)
{if (bitoffset) outputword|=(clearcode<<bitoffset);
else outputword=clearcode;
bitoffset+=maxbits;
inittable(initcodesize);
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
}
}
/* look for an unstored string */
length=1;
while (nextcode>
(currentcode=findcode(length,(char *)(ptr+currentplace))))
{prefixcode=currentcode;
length++;
if ((currentplace+length)>=size) break;
}
nextcode++;
currentplace+=(length-1);
/* output the prefix code */
if (bitoffset) outputword|=(prefixcode<<bitoffset);
else outputword=prefixcode;
bitoffset+=codesize;
/* output finished bytes to blocks */
while (bitoffset>=8)
{blockbuffer[blocksize]=outputword&0xff;
outputword>>=8;
bitoffset-=8;
blocksize++;
/* output filled block */
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
}
/* output the end of information code */
if (bitoffset) outputword|=(eoicode<<bitoffset);
else outputword=eoicode;
bitoffset+=codesize;
/* finish outputting the data */
while (bitoffset>=0)
{blockbuffer[blocksize]=(char) (outputword&0xff);
outputword>>=8;
bitoffset-=8;
blocksize++;
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
/* output the last block */
if (blocksize)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
}
/* a final zero block count */
fputc(0, outfile);
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void decompress(initcodesize,ptr,ptr1,size)
int initcodesize;
unsigned char *ptr, *ptr1; /* compressed data from ptr go to ptr1 */
int size; /* an upper limit purely as a check */
{int i,currentcode,codesize=0,maxbits=12,blocksize;
int clearcode,eoicode,codemask=0;
int bitoffset=0,indx,oldindx=0;
int currentplace=0,oldplace=0;
int findcode();
unsigned long inputword=0;
unsigned char *p1, *p2;
/* first deblock the data */
p1=p2=ptr;
blocksize=(*p1++);
while (blocksize)
{while (blocksize--)
(*p2++)=(*p1++); /* a wonderful example of how abstruse C can be */
blocksize=(*p1++);
}
/* set up initial code table */
currentcode=clearcode=(1<<initcodesize);
eoicode=clearcode+1;
/* allocate space for hash table */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(addresses=(char **) malloc(sizeof(char*)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentplace=eoicode;
}
while (currentcode!=eoicode)
{if (currentcode==clearcode) /* reset the decompressor */
{inittable(initcodesize);
codesize=initcodesize+1;
nextcode=eoicode+1;
codemask=(1<<codesize)-1;
oldindx=(-1);
}
else /* code represents data */
{indx=codeindex[currentcode]; /* where in table is currentcode */
if (indx>=0) /* it is there */
{ /* put it into the output */
for (i=0;i<lengths[indx];i++)
ptr1[currentplace+i]=addresses[indx][i];
if (oldindx>=0) /* first character treated differently */
{findcode(lengths[oldindx]+1,(char *) (ptr1+oldplace));
nextcode++; /* add new code to table */
}
oldplace=currentplace;
currentplace+=lengths[indx];
oldindx=indx;
}
else /* not in table yet; must be old code plus last=first character */
{for (i=0;i<lengths[oldindx];i++)
ptr1[currentplace+i]=addresses[oldindx][i];
ptr1[currentplace+lengths[oldindx]]=addresses[oldindx][0];
/* store new code */
findcode(lengths[oldindx]+1,(char *)(ptr1+currentplace));
oldplace=currentplace;
currentplace+=lengths[oldindx]+1;
oldindx=codeindex[nextcode++];
}
/* crude error checking */
if ((oldindx<0)||(currentplace>size))
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* check if codesize needs increasing */
if (nextcode>codemask)
if (codesize<maxbits)
{codesize++;
codemask=(1<<codesize)-1;
}
while (bitoffset<codesize) /* read some more data */
{if (bitoffset) inputword|=(((int)(*ptr++))<<bitoffset);
else inputword=(*ptr++);
bitoffset+=8;
}
/* strip off current code */
currentcode=inputword&codemask;
inputword>>=codesize;
bitoffset-=codesize;
if (currentcode>nextcode)
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void inittable(size)
int size;
{int i;
for (i=0;i<TABLELENGTH;i++)
{linktonext[i]=(-1);
codes[i]=(-1);
lengths[i]=(-1);
}
for (i=0;i<4096;i++)
codeindex[i]=(-1);
/* store initial codes for raw characters */
nextcode=0;
for (i=0;i<(1<<size);i++)
{stringbuffer[i]=i;
findcode(1,(char *) (stringbuffer+i));
nextcode++;
}
return;
}
int findcode(length,string)
char *string;
int length;
/* return code for string of given length;
if not found, store it and return nextcode */
{int i,j,indx,previousindex;
indx=hashit(length,string);
/* look for string in table */
previousindex=indx;
while (indx>0)
{if (lengths[indx]==length) /* is the length right ? */
{for (j=0;j<length;j++)
if ((string[j]) != (addresses[indx][j])) break;
if (j==length) return codes[indx];
}
previousindex=indx;
indx=linktonext[indx];
}
/* not found, so store it */
indx=previousindex;
i=indx;
while (codes[i]>=0) /* find an unused slot in table */
{i++;
if (i>=TABLELENGTH) i-=TABLELENGTH;
}
if (i!=indx)
{linktonext[indx]=i; /* link to it */
indx=i; /* move to it */
}
codes[indx]=nextcode; /* save the new code */
lengths[indx]=length;
addresses[indx]=string;
codeindex[nextcode]=indx;
return nextcode;
}
Raw
xpotts.txt
XPotts
Michael Creutz
[email protected]
WHAT IT IS:
This is an Xwindow program to illustrate a Monte Carlo simulation of
the Potts model. This grew from my earlier program XIsing. If you
have trouble understanding what is going on here, look at XIsing first.
This program allows dynamical observation of the Potts system under a
simple microcanonical Monte Carlo algorithms.
The system is a two dimensional array of spins, each of which is
represented by one of q possible states. The states are represented
by various colors displayed on the screen. The energy is determined
by nearest neighbor pairs as well as the value of an applied field.
If two neighbors are the same, then the energy has one value, and if
they differ, it is higher by one unit. The applied field acts only
on color 0, which is red, and gives it a slightly different energy.
The source is freely distributable.
COMPILING:
The program should run on anything supporting Xwindows with color.
If you find an exception, please let me know. It will not
run on a monochrome system.
To compile, try the xtoys Makefile. Otherwise, "cc -O xpotts.c -lm
-lX11" should work. If not, possibly the X11 includes are not being
found. Then you need to compile with a -I option to where they are
and possibly change the -lX11 to help the linker. Remember, as with
any Xwindow program, you need to have "xhost" and "DISPLAY" set up
properly.
WHAT THE SCREEN SHOWS:
The spins are displayed in a box surrounded by various gadgets
to control the evolution. Every 20 sweeps, various quantities
are calculated and displayed to the right of the lattice. These
include the inverse temperature, beta, the applied field h, the
system energy E, and the magnitization as defined from state 0.
The normalization is such that at h=0 the energy E runs from 1
to 0 as beta runs from 0 to infinity. At beta=0, the magnetization
M runs from 0 to 1 as h runs from zero to infinity.
The critical value for beta at h=0 is exactly known to be
log(1+sqrt(q)). This value is displayed on the screen below
the measured quantities.
One set of four buttons, labeled "conserve E," "heat," "cool,"
and "free E" are for adjusting the temperature. Normal running
is with the energy conserved, but to heat or cool the system
press the appropriate button. The "free E" button is similar to the
heater; the difference is that the heat button will not stop at infinite
temperature but will continue to negative beta.
Another set of buttons controls the magnetization, defined as the
number of sites with spins having value 0. A conventional simulation
would have this unconstrained, and is represented by the "free M"
button. Constraining M is the way this program implements a magnetic
field.
The basic algorithm is a microcanonical approach described in some
detail in my paper Phys. Rev. Letters 50, 1411 (1983). A set of "demons"
circulates around the lattice trying to flip spins. Each carries a
sack of energy. Any energy change associated with a spin flip is
compensated by a change in this sack. If the demon's sack cannot
accommodate the change, the flip is rejected. For the control of the
magnitization, a second sack contains a limited supply of state 0 spins
from which any changes to of from such state must be taken.
The behavior under this algorithm is quite close to that of a conventional
Metropolis et al. simulation. This program is substantially slower
than XIsing, which did all computations in parallel using multi-spin
coding.
The "quit" button and the "pause" "run" toggle are self explanatory.
The number of states allowed per site, q, can be adjusted via
the box just below them. The allowed range is from 2 to 255, although
only 18 different colors are used. As q is decreased, the highest spins
are replaced by the next value down.
The program can be run with one or two arguments. The first argument
is the initial value of q, which defaults to 8. The second is used
to make the display bigger. Its value is the size of the blocks
of pixels used. The best way to see what this argument does is to try
it, i.e. enter xpotts 2 4 to get the Ising model with big blocks.
It defaults to 1, but can be as large as 4. The larger
sizes slow things down a lot on a SUN, but may work better on other
machines.
SUGGESTED EXPERIMENTS:
The system starts at infinite temperature with the cooler on.
At large q, including the default of 8, there is a first order
transition. Note how as the system cools, blocks of fixed color
appear. The high temperature phase appears in more random regions
between them. As it cools, the system tends to spend a lot of
time near the critical temperature as the transition to the
nucleating low temperature phase releases latent heat.
If you switch to the conserve E mode before the low temperature phase
has completely taken over, the system will be phase separated
and the temperature should fluctuate around the critical value.
With the magnetism buttons you can use the spin zero magnetized
state dominated to compress the other phases. Note how this raises
the temperature at constant energy. Releasing the magnetic
pressure, you can watch the other states bounce back.
At q<=4 the transition is second order; observe how the fluctuations
in the domains are much larger near the critical temperature. The case
q=5 is supposed to be first order with a small latent heat. Can
you tell?
The heat button allows one to go to negative temperature.
Going down to q=2 gives the Ising model. It is hard to
see, but in this case leaving the heat button on will lead to
antiferromagnetic domains in a checkerboard pattern. Note how
resiliant this state is to an applied field. At higher q the
negative temperature system has lots of degeneracy and looks
like a mess.
Turning on the increase M button will favor color zero. At low
temperature alternate use of this and the decrease M button
allows one to pass through a first order transition from one
magnitization state to another. Leaving on the decrease M button
will freeze out color zero, giving the model effectively at
q-1. Try this for the q=3 model and watch the Ising model appear.
Raw
xsand.c
/* Xsand, a sandpile simulator for X windows.
Michael Creutz
[email protected]
compiling: cc -O -o xsand -L/usr/X11R6/lib xsand.c -lX11
version 1.6
August 2006
The latest version is kept at "http://thy.phy.bnl.gov/www/xtoys/xtoys.html"
A text description of this program is there as well.
*/
/* because of the unisys patent on LZW, it may be illegal to
distribute this program with the following "if" modified to "while"*/
# include <X11/Xlib.h>
# include <X11/Xutil.h>
# include <X11/Xos.h>
# include <X11/Xatom.h>
# include <X11/cursorfont.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <malloc.h>
# define BPW (sizeof(long))
/* size and position parameters for buttons, etc. */
# define PLAYTOP 116
# define PLAYLEFT 10
# define BOUNDHEIGHT 72
# define TRACEHEIGHT 36
# define BOUNDWIDTH 112
# define BARHEIGHT 15
# define BARWIDTH 136
# define BARTOP 8
# define BARLEFT 40
# define SRHEIGHT (3*18)
# define BUTTONWIDTH 68
int nrows,ncols,volume; /* lattice dimensions (plus two for boundaries) */
/* ncols will be truncated to
a multiple of bytes per long word (BPW),
and then two units are lost for boundaries */
int block, /* size of cells displayed on the screen */
blockshift; /* log_2(block) */
unsigned char *field[2]={NULL,NULL}; /* pointers to old and new system */
int old=0,new=1; /* which field is current? */
int paused=0, /* is the system paused? */
boundary, /* 0,1,2,3 for open, live, periodic, flow */
pencolor=4, /* for sketching */
trace, /* is the tracer on? */
colshift, /* how to shift to get the next column */
active=1,
autodouble=0;
long mask0x13,mask0x01,mask0x03,mask0x10; /* some useful masks */
char stringbuffer[256]; /* generally useful */
/* things for regulating the updating speed */
int speed=0; /* updating delay proportional to speed */
int delay=0; /* counter for implementing speed control */
struct timeval timeout; /* timer for speed control */
/* various window stuff */
Display *display;
int screen;
static char *progname;
Window window,quitbutton,pausebutton,playground,srbutton,blockbutton,
boundwindow,colorbar,fillbutton,tracewindow,doublebutton,adbutton,
speedbutton,makebutton();
XColor xcolor,colorcell;
Colormap cmap;
GC gc,gcpen;
int windowwidth,windowheight;
XFontStruct *font=NULL;
int font_height,font_width;
XSizeHints size_hints;
int darkcolor,lightcolor,black,white;
XImage *spinimage=NULL;
long translate[256]; /* for converting colors */
char * boundtext[4]={"open","sandy","periodic","flow"};
char * srtext[3]={"save","restore","+ saved"};
void drawbutton(),settrace(),openwindow(),makebuttons(),update(),repaint(),
cleanup(),fixboundary(),showpic(),setboundary(),sketch(),savepic(),
addpic(),setpen(),parallelupdate(),decompress(),compress(),inittable();
int main(argc,argv)
int argc;
char **argv;
{unsigned int width, height;
int i,x,y;
long bitcheck;
XEvent report;
progname=argv[0];
/* initial lattice size */
block=1;
nrows=200/block;
ncols=(~(BPW-1))&(200/block);
volume=nrows*ncols;
trace=0;
boundary=0;
/* open the window and make the buttons */
openwindow(argc,argv);
makebuttons();
/* set the masks for parallel updating */
for (i=0;i<sizeof(long);i++)
{mask0x13=0x13|(mask0x13<<8);
mask0x01=0x01|(mask0x01<<8);
mask0x03=0x03|(mask0x03<<8);
mask0x10=0x10|(mask0x10<<8);
}
/* check for byte order in a word */
bitcheck=1;
if (* (char*) (&bitcheck) )
colshift=(-1); /* LSBfirst */
else
colshift=1; /* MSBfirst */
/* loop forever, looking for events */
while(1)
{if (0==paused)
{if (delay) /* don't update yet */
{delay--;
/* this use of select() seems a kludge to me;
why can't usleep() be more standard? */
timeout.tv_sec=0;
timeout.tv_usec=100000; /* .1 sec per delay unit */
select(0,NULL,NULL,NULL,&timeout);
}
else
{delay=speed;
update();
}
}
if (paused|XPending(display)|((active|autodouble)==0))
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if ((report.xexpose.window)!=window) break; /* cuts down flashing,
but you might remove this line if things aren't being redrawn */
if (report.xexpose.count!=0) break; /* more in queue, wait for them */
repaint();
break;
case ConfigureNotify:
width=report.xconfigure.width;
height=report.xconfigure.height;
if ((width<size_hints.min_width)||(height<size_hints.min_height))
{fprintf(stderr,"%s: window too small to proceed.\n",progname);
cleanup();
}
else if ((width!=windowwidth)||(height!=windowheight))
{windowwidth=width;
windowheight=height;
XSetForeground(display,gcpen,lightcolor);
XFillRectangle(display,window,gcpen,0,0,windowwidth,windowheight);
nrows=(height-PLAYTOP-10)/block;
ncols=(~(BPW-1))&((width-PLAYLEFT-BOUNDWIDTH-20)/block);
volume=nrows*ncols;
makebuttons();
fixboundary();
showpic();
}
break;
case ButtonPress:
if (report.xbutton.window==quitbutton)
cleanup();
else if (report.xbutton.window==pausebutton)
{paused=1-paused;
drawbutton(pausebutton,0,0,BUTTONWIDTH,18,"pause",1-2*paused);
}
else if (report.xbutton.window==fillbutton)
{drawbutton(fillbutton,0,0,BUTTONWIDTH,18,"fill",-1);
for (i=0;i<volume;i++)
field[old][i]=pencolor;
fixboundary();
showpic();
drawbutton(fillbutton,0,0,BUTTONWIDTH,18,"fill",1);
}
else if (report.xbutton.window==boundwindow)
{setboundary((report.xbutton.y*4)/BOUNDHEIGHT);
showpic();
}
else if (report.xbutton.window==tracewindow)
settrace((report.xbutton.y*2)/TRACEHEIGHT);
else if (report.xbutton.window==playground) /* sketch */
{x=report.xbutton.x/block;
y=report.xbutton.y/block;
sketch(x,y);
}
else if (report.xbutton.window==srbutton) /* save or restore */
{y=(report.xbutton.y*3)/SRHEIGHT;
drawbutton(srbutton,0,y*(SRHEIGHT/3),BUTTONWIDTH,SRHEIGHT/3,
srtext[y],-1);
if (0==y)
savepic(field[old],ncols,nrows);
else if (1==y) /* load in gif image */
{for (i=0;i<volume;i++)
field[old][i]=0;
addpic(field[old],ncols,nrows);
showpic();
}
else /* add gif image mod 8 to existing field */
{addpic(field[old],ncols,nrows);
for (i=0;i<volume;i++)
field[old][i]&=7;
showpic();
}
drawbutton(srbutton,0,y*(SRHEIGHT/3),BUTTONWIDTH,SRHEIGHT/3,
srtext[y],1);
}
else if (report.xbutton.window==blockbutton)
{i=(1<<(4*report.xbutton.x/BOUNDWIDTH));
if (i!=block) /* reset for new block size */
{block=i;
nrows=(windowheight-PLAYTOP-10)/block;
ncols=(~(BPW-1))&((windowwidth-PLAYLEFT-BOUNDWIDTH-20)/block);
volume=nrows*ncols;
XSetForeground(display,gcpen,lightcolor);
XFillRectangle(display,window,gcpen,0,0,
windowwidth,windowheight);
makebuttons();
fixboundary();
showpic();
}
}
else if (report.xbutton.window==colorbar)
{setpen((report.xbutton.x*8)/BARWIDTH);
}
else if (report.xbutton.window==adbutton)
{autodouble=1-autodouble;
drawbutton(adbutton,0,0,BUTTONWIDTH,18,"auto-d",1-2*autodouble);
}
else if (report.xbutton.window==doublebutton)
{drawbutton(doublebutton,0,0,BUTTONWIDTH,18,"double",-1);
for (i=0;i<volume;i++)
field[old][i]=7&(field[old][i]<<1);
active=1;
fixboundary();
showpic();
drawbutton(doublebutton,0,0,BUTTONWIDTH,18,"double",1);
}
else if (report.xbutton.window==speedbutton)
{ /* reset speed */
speed=10-(11*report.xbutton.x)/BOUNDWIDTH;
if (speed<0) speed=0;
if (speed>10) speed=10;
delay=speed;
drawbutton(speedbutton,0,0,BOUNDWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(BOUNDWIDTH-2))/11,1,
(BOUNDWIDTH-2)/11,16,"",2);
}
else
update(); /* do a sweep when mouse clicked */
break;
default:
break;
} /* end of switch */
} /* end of if XPending */
} /* end of while(1) */
} /* end of main */
void sketch(x,y) /* sketch in playground starting at point x,y */
int x,y;
{int sketching,i,oldx,oldy,deltax,deltay,newx,newy;
Window root,child;
XEvent sketchreport;
unsigned int keys_buttons;
int window_x,window_y;
newx=x;
newy=y;
field[old][x+ncols*y]=pencolor;
sketching=1;
while (sketching)
{showpic();
XNextEvent(display,&sketchreport);
switch (sketchreport.type)
{case MotionNotify:
oldx=newx;
oldy=newy;
if (!XQueryPointer(display,playground,
&root,&child,&window_x,&window_y,
&newx,&newy,&keys_buttons))
{sketching=0;
break;
}
if (blockshift)
{newx>>=blockshift;
newy>>=blockshift;
}
if ((newx>ncols)|(newy>nrows)|(newx<0)|(newy<0))
{sketching=0;
break;
}
deltax=abs(newx-oldx);
deltay=abs(newy-oldy);
if (deltax>deltay)
for (i=1;i<=deltax;i++)
field [old][oldx+i*(newx-oldx)/deltax
+ncols*(oldy+i*(newy-oldy)/deltax)]=pencolor;
else
for (i=1;i<=deltay;i++)
field [old][oldx+i*(newx-oldx)/deltay
+ncols*(oldy+i*(newy-oldy)/deltay)]=pencolor;
break;
case ButtonRelease:
default:
sketching=0;
break;
}
}
fixboundary();
showpic();
return;
}
void update()
{int i;
parallelupdate(field[old],field[new]);
old=new;
new=1-new;
fixboundary();
showpic();
if ((0==active)&&autodouble)
{for (i=0;i<volume;i++)
field[old][i]=7&(field[old][i]<<1);
fixboundary();
showpic();
}
return;
}
void parallelupdate(source,destination)
/* do several spins at a time */
long * source;
long * destination;
{int i,maxi=volume/BPW;
active=0;
for (i=0;i<maxi;i++)
{destination[i]=source[i]&mask0x13;
active|=(source[i]=(source[i]>>2)&mask0x01);
}
if (trace)
for (i=0;i<maxi;i++)
destination[i]|=source[i]<<4;
maxi=(volume-ncols)/BPW;
for (i=ncols/BPW;i<maxi;i++)
destination[i]+=(source[i ]>>8)
+(source[i-colshift]<<(8*(BPW-1)))
+(source[i ]<<8)
+(source[i+colshift]>>(8*(BPW-1)))
+ source[i+ncols/BPW]
+ source[i-ncols/BPW];
return;
}
void showpic() /* display the field */
{int row,col,i1,i2,color,j,j1,j2,blocktop=block;
char *picture=(*spinimage).data;
if (block>4) blocktop=block-1;
if (8==(*spinimage).depth)
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[old][row*ncols+col]];
j=block*(col+block*ncols*row);
if (color!=picture[j])
for (i1=0;i1<blocktop;i1++)
{j1=i1*block*ncols+j;
for (i2=0;i2<blocktop;i2++)
picture[j1+i2]=color;
}
}
else
{for (j=0;j<volume;j++)
picture[j]=translate[field[old][j]];
}
}
else /* depth is not 8, use xputpixel (this is really ugly) */
{if (block>1) /* I wish I knew how to do this faster */
for (row=0;row<nrows;row++)
for (col=0;col<ncols;col++)
{color=translate[field[old][row*ncols+col]];
if (color!=XGetPixel(spinimage,j1=block*col,j2=block*row))
for (i2=0;i2<blocktop;i2++)
for (i1=0;i1<blocktop;i1++)
XPutPixel(spinimage,j1+i1,j2+i2,color);
}
else
for (row=0;row<nrows;row++)
{j1=row*ncols;
for (col=0;col<ncols;col++)
XPutPixel(spinimage,col,row,translate[field[old][j1+col]]);
}
}
XPutImage(display,playground,gc,spinimage,0,0,0,0,block*ncols,block*nrows);
return;
}
void fixboundary()
{int i,bv;
if (boundary<2) /* dead or alive */
{bv=4*boundary;
for (i=0;i<ncols;i++)
{field[old][i]=bv;
field[old][volume-1-i]=bv;
}
for (i=0;i<volume;i+=ncols)
{field[old][i]=bv;
field[old][i+ncols-1]=bv;
}
}
else if (2==boundary) /* periodic */
{for (i=0;i<ncols;i++)
{field[old][i]=field[old][volume-2*ncols+i];
field[old][volume-ncols+i]=field[old][ncols+i];
}
for (i=0;i<volume;i+=ncols)
{field[old][i]=field[old][i+ncols-2];
field[old][i+ncols-1]=field[old][i+1];
}
}
else /* flow: alive on top, dead on sides and bottom */
{for (i=0;i<ncols;i++)
{field[old][i]=4;
field[old][volume-1-i]=0;
}
for (i=0;i<volume;i+=ncols)
{field[old][i]=0;
field[old][i+ncols-1]=0;
}
}
return;
}
void repaint()
/* this fixes the window up whenever it is uncovered */
{int i;
drawbutton(pausebutton,0,0,BUTTONWIDTH,18,"pause",1-2*paused);
drawbutton(quitbutton,0,0,BUTTONWIDTH,18,"quit",1);
drawbutton(fillbutton,0,0,BUTTONWIDTH,18,"fill",1);
drawbutton(doublebutton,0,0,BUTTONWIDTH,18,"double",1);
drawbutton(adbutton,0,0,BUTTONWIDTH,18,"auto-d",1-2*autodouble);
drawbutton(window,PLAYLEFT-4,PLAYTOP-4,block*ncols+8,block*nrows+8,NULL,4);
drawbutton(speedbutton,0,0,BOUNDWIDTH,18,"speed",-1);
drawbutton(speedbutton,1+((10-speed)*(BOUNDWIDTH-2))/11,1,
(BOUNDWIDTH-2)/11,16,"",2);
for (i=0;i<3;i++)
drawbutton(srbutton,0,i*(SRHEIGHT/3),BUTTONWIDTH,SRHEIGHT/3,srtext[i],1);
/* fix the block buttons */
for (i=1;i<=4;i++)
{sprintf(stringbuffer,"%d",(1<<i)/2);
drawbutton(blockbutton,(i-1)*BOUNDWIDTH/4,0,BOUNDWIDTH/4,18,
stringbuffer,1-2*(i==(blockshift+1)));
}
XDrawString(display,window,gc,BARLEFT+24,BARTOP+68,"cell size",9);
/* draw the boundary condition buttons */
for (i=0;i<4;i++)
drawbutton(boundwindow,0,i*BOUNDHEIGHT/4,BOUNDWIDTH,BOUNDHEIGHT/4,
boundtext[i],1-2*(i==boundary));
XDrawString(display,window,gc,
windowwidth-BOUNDWIDTH,BOUNDHEIGHT+font_height+2," Boundary ",12);
/* draw the trace buttons */
drawbutton(tracewindow,0,0,BOUNDWIDTH,TRACEHEIGHT/2,
"tracer",1-2*trace);
drawbutton(tracewindow,0,TRACEHEIGHT/2,BOUNDWIDTH,TRACEHEIGHT/2,
"clear trace",1);
/* fix the colorbar */
for (i=0;i<8;i++)
{drawbutton(colorbar,i*BARWIDTH/8,0,
BARWIDTH/8,BARHEIGHT,NULL,2-4*(i==pencolor));
XSetForeground(display,gcpen,translate[i]);
XFillRectangle(display,colorbar,gcpen,i*BARWIDTH/8+2,2,
BARWIDTH/8-4,BARHEIGHT-4);
sprintf(stringbuffer,"%d",i);
XDrawString(display,window,gc,
4+BARLEFT+BARWIDTH*i/8,BARTOP+28,stringbuffer,1);
}
/* write various strings */
sprintf(stringbuffer,"%d by %d lattice ",ncols-2,nrows-2);
XDrawImageString(display,window,gc,BARLEFT-14,100,
stringbuffer,strlen(stringbuffer));
XDrawString(display,window,gc,windowwidth-44,windowheight-5
,"MJC",3);
showpic();
return;
}
void setboundary(value)
/* take action if boundary buttons hit */
int value;
{drawbutton(boundwindow,0,
(BOUNDHEIGHT*boundary)/4,BOUNDWIDTH,BOUNDHEIGHT/4,
boundtext[boundary],1);
boundary=value;
fixboundary();
drawbutton(boundwindow,0,
(BOUNDHEIGHT*boundary)/4,BOUNDWIDTH,BOUNDHEIGHT/4,
boundtext[boundary],-1);
return;
}
void settrace(value)
/* take action if trace buttons hit */
int value;
{int i;
if (1==value) /* clear the trace bit */
{drawbutton(tracewindow,0,TRACEHEIGHT/2,BOUNDWIDTH,TRACEHEIGHT/2,
"clear trace",-1);
for (i=0;i<volume;i++)
field[old][i]&=0xf;
showpic();
drawbutton(tracewindow,0,TRACEHEIGHT/2,BOUNDWIDTH,TRACEHEIGHT/2,
"clear trace",1);
}
else
{trace=1-trace;
drawbutton(tracewindow,0,0,BOUNDWIDTH,TRACEHEIGHT/2,
"tracer",1-2*trace);
}
return;
}
void setpen(value)
/* take action if pen buttons hit */
int value;
{drawbutton(colorbar,pencolor*BARWIDTH/8,0,
BARWIDTH/8,BARHEIGHT,NULL,2);
XSetForeground(display,gcpen,translate[pencolor]);
XFillRectangle(display,colorbar,gcpen,pencolor*BARWIDTH/8+2,2,
BARWIDTH/8-4,BARHEIGHT-4);
pencolor=value;
drawbutton(colorbar,pencolor*BARWIDTH/8,0,
BARWIDTH/8,BARHEIGHT,NULL,-2);
XSetForeground(display,gcpen,translate[pencolor]);
XFillRectangle(display,colorbar,gcpen,pencolor*BARWIDTH/8+2,2,
BARWIDTH/8-4,BARHEIGHT-4);
return;
}
void openwindow(argc,argv)
/* a lot of this is taken from basicwin in the Xlib Programming Manual */
int argc;
char **argv;
{char *window_name="Sand";
char *icon_name="Sand";
long event_mask;
Pixmap icon_pixmap;
char *display_name=NULL;
int i;
# define icon_bitmap_width 16
# define icon_bitmap_height 16
static unsigned char icon_bitmap_bits[] = {
0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
/* open up the display */
if ((display=XOpenDisplay(display_name))==NULL)
{fprintf(stderr,"%s: cannot connect to X server %s\n",
progname,XDisplayName(display_name));
exit(-1);
}
screen=DefaultScreen(display);
cmap=DefaultColormap(display,screen);
darkcolor=black=BlackPixel(display,screen);
lightcolor=white=WhitePixel(display,screen);
translate[0]=white;
translate[1]=black;
translate[2]=lightcolor;
translate[3]=darkcolor;
if (XAllocNamedColor(display,cmap,"salmon",&colorcell,&xcolor))
darkcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"wheat",&colorcell,&xcolor))
lightcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"grey",&colorcell,&xcolor))
translate[0]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"red",&colorcell,&xcolor))
translate[1]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"blue",&colorcell,&xcolor))
translate[2]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"green",&colorcell,&xcolor))
translate[3]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"orange",&colorcell,&xcolor))
translate[4]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"gold",&colorcell,&xcolor))
translate[5]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"yellow",&colorcell,&xcolor))
translate[6]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"white",&colorcell,&xcolor))
translate[7]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"turquoise",&colorcell,&xcolor))
translate[8]=colorcell.pixel;
/* fill out the color table for future uses */
for(i=9;i<256;i++)
translate[i]=translate[i%8];
for (i=16;i<20;i++)
translate[i]=translate[8]; /* for tracing */
/* make the main window */
windowwidth=(block*ncols+PLAYLEFT+BOUNDWIDTH+20);
windowheight=(PLAYTOP+block*nrows+10);
window=XCreateSimpleWindow(display,RootWindow(display,screen),
0,0,windowwidth,windowheight,4,black,lightcolor);
/* make the icon */
icon_pixmap=XCreateBitmapFromData(display,window,
icon_bitmap_bits,icon_bitmap_width,icon_bitmap_height);
size_hints.flags=PPosition | PSize | PMinSize;
size_hints.min_width=windowwidth;
size_hints.min_height=windowheight;
#ifdef X11R3
size_hints.x=x;
size_hints.y=y;
size_hints.width=windowwidth;
size_hints.height=windowheight;
XSetStandardProperties(display,window,window_name,icon_name,
icon_pixmap,argv,argc,&size_hints);
#else
{XWMHints wm_hints;
XClassHint class_hints;
XTextProperty windowName, iconName;
if (XStringListToTextProperty(&window_name,1,&windowName)==0)
{fprintf(stderr,"%s: structure allocation for windowName failed.\n"
,progname);
exit(-1);
}
if (XStringListToTextProperty(&icon_name,1,&iconName)==0)
{fprintf(stderr,"%s: structure allocation for iconName failed.\n"
,progname);
exit(-1);
}
wm_hints.initial_state=NormalState;
wm_hints.input=True;
wm_hints.icon_pixmap=icon_pixmap;
wm_hints.flags=StateHint|IconPixmapHint|InputHint;
class_hints.res_name=progname;
class_hints.res_class="Basicwin";
XSetWMProperties(display,window,&windowName,&iconName,
argv,argc,&size_hints,&wm_hints,&class_hints);
}
#endif
/* pick the events to look for */
event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
XSelectInput(display,window,event_mask);
/* pick font: 9x15 is supposed to almost always be there,
the fixed font option put in by Peter Chang */
if ((font=XLoadQueryFont(display,"9x15"))==NULL)
if ((font=XLoadQueryFont(display,"fixed"))==NULL)
{fprintf(stderr,"%s: Sorry, having font problems.\n",progname);
exit(-1);
}
font_height=font->ascent+font->descent;
font_width=font->max_bounds.width;
/* make graphics contexts:
gc for black on white
gcpen for varying purposes */
gc=XCreateGC(display,window,0,NULL);
XSetFont(display,gc,font->fid);
XSetForeground(display,gc,black);
XSetBackground(display,gc,lightcolor);
gcpen=XCreateGC(display,window,0,NULL);
XSetFont(display,gcpen,font->fid);
XSetForeground(display,gcpen,translate[pencolor]);
XSetBackground(display,gcpen,lightcolor);
/* show the window */
XMapWindow(display,window);
return;
}
void makebuttons()
{int i,j;
long event_mask;
XEvent report;
Cursor cursor;
/* first destroy any old buttons */
XDestroySubwindows(display,window);
/* now make the new buttons */
quitbutton =makebutton(windowwidth-92,windowheight-38 ,BUTTONWIDTH,18);
pausebutton =makebutton(windowwidth-92,windowheight-38- 18,BUTTONWIDTH,18);
fillbutton =makebutton(windowwidth-92,windowheight-38-2*18,BUTTONWIDTH,18);
doublebutton=makebutton(windowwidth-92,windowheight-38-3*18,BUTTONWIDTH,18);
adbutton =makebutton(windowwidth-92,windowheight-38-4*18,BUTTONWIDTH,18);
srbutton =makebutton(windowwidth-92,windowheight-38-7*18,BUTTONWIDTH,18*3);
boundwindow =makebutton(windowwidth-BOUNDWIDTH-4,4,BOUNDWIDTH,BOUNDHEIGHT);
tracewindow =makebutton(windowwidth-BOUNDWIDTH-4,BOUNDHEIGHT+25,
BOUNDWIDTH,TRACEHEIGHT);
colorbar =makebutton(BARLEFT,BARTOP,BARWIDTH,BARHEIGHT);
blockbutton =makebutton(BARLEFT+12,BARTOP+36,BOUNDWIDTH,18);
speedbutton=makebutton(windowwidth-BOUNDWIDTH-4,
BOUNDHEIGHT+TRACEHEIGHT+30,
BOUNDWIDTH,18);
playground=XCreateSimpleWindow(display,window,
PLAYLEFT,PLAYTOP,block*ncols,block*nrows,0,black,translate[8]);
event_mask=ExposureMask|ButtonReleaseMask|ButtonPressMask|
PointerMotionHintMask|ButtonMotionMask;
XSelectInput(display,playground,event_mask);
XMapWindow(display,playground);
/* wait for playground to be displayed befor proceeding */
i=1; /* a flag */
while (i)
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if (report.xexpose.window!=playground) i=0;
default:
break;
}
}
/* make image structure */
if (NULL!=spinimage)
{XDestroyImage(spinimage);
spinimage=NULL;
}
spinimage=XGetImage((Display *) display, (Drawable) playground,
0,0,block*ncols,block*nrows,AllPlanes,ZPixmap);
if (NULL==spinimage)
{fprintf(stderr,"trouble creating image structure\n");
exit(-1);
}
/* make special cursors to be cute */
cursor=XCreateFontCursor(display,XC_sb_up_arrow);
XDefineCursor(display,playground,cursor);
cursor=XCreateFontCursor(display,XC_hand2);
XDefineCursor(display,doublebutton,cursor);
XDefineCursor(display,colorbar,cursor);
XDefineCursor(display,boundwindow,cursor);
XDefineCursor(display,tracewindow,cursor);
XDefineCursor(display,srbutton,cursor);
XDefineCursor(display,blockbutton,cursor);
XDefineCursor(display,fillbutton,cursor);
XDefineCursor(display,adbutton,cursor);
XDefineCursor(display,fillbutton,cursor);
XDefineCursor(display,quitbutton,cursor);
XDefineCursor(display,pausebutton,cursor);
XDefineCursor(display,speedbutton,cursor);
/* reallocate various arrays */
/* allocate a few bytes of extra space at ends of field array */
for (i=0;i<2;i++)
{if (NULL!=field[i]) free((char *) field[i]-BPW);
if (NULL==( field[i]= (unsigned char *) malloc(volume+2*BPW)))
{fprintf(stderr,"allocation problems\n");
cleanup();
}
/* clear initial state */
for (j=0;j<(volume+2*BPW);j++)
field[i][j]=pencolor;
field[i]+=BPW;
}
/* set blockshift to log_2 of block */
blockshift=0;
i=block;
while(i>>=1)
blockshift++;
/* draw everything */
repaint();
return;
}
Window makebutton(xoffset,yoffset,xsize,ysize)
int xoffset,yoffset,xsize,ysize;
/* Puts button of specified dimensions on window. Nothing is drawn. */
{Window buttonwindow;
long event_mask;
buttonwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
xsize,ysize,0,black,lightcolor);
event_mask=ButtonPressMask|ExposureMask; /* look for mouse-button presses */
XSelectInput(display,buttonwindow,event_mask);
XMapWindow(display,buttonwindow);
return buttonwindow;
}
void drawbutton(buttonwindow,xoffset,yoffset,xsize,ysize,text,state)
Window buttonwindow;
int xoffset,yoffset,xsize,ysize,state;
char * text;
/* Draw a button in buttonwindow of specified dimensions with text
centered. Color of button determined by sign of "state."
size of border determined by magnitude. */
{int textlength,i,j;
int cdark,clight,cup,cdown;
int cleft,cright,cbutton,ctext;
cup=lightcolor;
cdown=darkcolor;
cdark=black;
clight=white;
if (state<0)
{cbutton=cdown;
ctext=clight;
cleft=cdark;
cright=clight;
}
else
{cbutton=cup;
ctext=cdark;
cleft=clight;
cright=cdark;
}
j=abs(state);
XSetForeground(display,gcpen,cbutton);
XFillRectangle(display,buttonwindow,gcpen,xoffset+j,yoffset+j,
xsize-2*j,ysize-2*j);
XSetForeground(display,gcpen,cleft);
XFillRectangle(display,buttonwindow,gcpen,xoffset,yoffset,xsize,j);
XFillRectangle(display,buttonwindow,gcpen,xoffset,yoffset,j,ysize);
XSetForeground(display,gcpen,cright);
for (i=0;i<j;i++)
{XDrawLine(display,buttonwindow,gcpen,
xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
XDrawLine(display,buttonwindow,gcpen,
xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
}
if (NULL!=text)
{textlength=strlen(text);
XSetForeground(display,gcpen,ctext);
XDrawString(display,buttonwindow,gcpen,
xoffset+(xsize-textlength*font_width)/2,
yoffset+(ysize+font_height)/2-2,
text,textlength);
}
return;
}
void cleanup()
{XUnloadFont(display,font->fid);
XFreeGC(display,gc);
XFreeGC(display,gcpen);
XCloseDisplay(display);
XDestroyImage(spinimage);
if (NULL!=field[0]) free((char *) field[0]-BPW);
if (NULL!=field[1]) free((char *) field[1]-BPW);
exit(1);
}
/* from here on is the stuff for saving and restoring from a gif file */
/* for info on how gif works, see ftp://network.ucsd.edu/graphics/GIF.shar.Z */
char *picturename="xsand.gif";
void addpic(data,xsize,ysize) /* add in a GIF image to field */
unsigned char *data; /* where the output data starts */
int xsize,ysize; /* output picture dimensions */
/* load in a gif image */
{int i,j,filesize,gwidth,gheight,gvolume;
unsigned char *ptr, *ptr1, *rawgif;
int colorbits,codesize;
FILE *infile;
if (NULL==(infile=fopen(picturename,"r")))
{fprintf(stderr,"couldn't open input file\n");
return;
}
/* find the file size */
fseek(infile, 0L, 2);
filesize = ftell(infile);
fseek(infile, 0L, 0);
/* make a place in memory for the file */
if (NULL==(rawgif= (unsigned char *) malloc(filesize) ))
{fprintf(stderr, "not enough memory to read gif file\n");
return;
}
ptr=rawgif;
/* read in the file */
if (fread(ptr, filesize, 1, infile) != 1)
{fprintf(stderr, "read failed\n");
free((char *) rawgif);
return;
}
fclose(infile);
/* check for GIF signature */
if (strncmp((char *) ptr,"GIF87a", 6))
{fprintf(stderr, "not a GIF87a file\n");
free((char *) rawgif);
return;
}
ptr+=6;
ptr+=4; /* skip over screen size */
colorbits=1+((*ptr)&0xf); /* how many bits of color */
ptr+=2; /* skip over background */
if (*ptr++) /* should be zero */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=(3*(1<<colorbits)); /* skip over colormap */
if (','!=(*ptr++)) /* image separator */
{fprintf(stderr, "corrupt GIF file\n");
free((char *) rawgif);
return;
}
ptr+=4; /* skip over image offset */
gwidth=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
gheight=(*ptr)+0x100*(*(ptr+1));
ptr+=2;
if (0x80&(*ptr++)) /* skip over local color map */
ptr+=(3*(1<<colorbits));
/* should now be at start of data */
codesize=(*ptr++);
/* make a place for the decompressed file */
gvolume=gwidth*gheight;
ptr1=(unsigned char *) malloc(gvolume);
decompress(codesize,ptr,ptr1,gvolume);
free((char *) rawgif);
/* map picture into data, allowing for different dimensions */
for (j=0;j<ysize;j++)
{if (j>=gheight) break;
for (i=0;i<xsize;i++)
{if (i>=gwidth) break;
data[i+j*xsize]+=ptr1[i+j*gwidth];
}
}
free((char *) ptr1);
fixboundary();
return;
}
void savepic(data,xsize,ysize) /* save the field as a GIF image */
unsigned char *data; /* where the input data starts */
int xsize,ysize; /* picture dimensions */
{int i;
int colorbits=5,codesize=5; /* assume a 32 color image */
FILE *outfile;
if (NULL==(outfile=fopen(picturename,"w")))
{fprintf(stderr,"couldn't open output file\n");
return;
}
/* GIF signature */
fwrite("GIF87a",6,1,outfile);
/* screen descriptor */
stringbuffer[0]=xsize&0xff; /* screen width */
stringbuffer[1]=(xsize>>8)&0xff;
stringbuffer[2]=ysize&0xff; /* screen height */
stringbuffer[3]=(ysize>>8)&0xff;
stringbuffer[4]=(0x80) /* M=1; global color map follows */
|((colorbits-1)<<4) /* -1+ bits of color reslution */
|(colorbits-1); /* -1+bits per pixel in image */
stringbuffer[5]=0; /* background color */
stringbuffer[6]=0; /* should be zero */
fwrite(stringbuffer,7,1,outfile);
/* global color map */
for (i=0;i<(1<<colorbits);i++)
{colorcell.pixel=translate[i];
XQueryColor(display,cmap,&colorcell);
fputc(colorcell.red>>8,outfile);
fputc(colorcell.green>>8,outfile);
fputc(colorcell.blue>>8,outfile);
}
/* image descriptor */
stringbuffer[0]=','; /* image descriptor separator */
stringbuffer[1]=0; /* image offset */
stringbuffer[2]=0;
stringbuffer[3]=0;
stringbuffer[4]=0;
stringbuffer[5]=xsize&0xff; /* image width */
stringbuffer[6]=(xsize>>8)&0xff;
stringbuffer[7]=ysize&0xff; /* image height */
stringbuffer[8]=(ysize>>8)&0xff;
stringbuffer[9]=0; /* use global color map, no interlace */
fwrite(stringbuffer,10,1,outfile);
/* start of image data */
fputc(codesize,outfile);
compress(codesize,data,outfile,volume);
/* gif terminator */
fputc(';',outfile);
fclose(outfile);
return;
}
/* LZW compression */
/* hash function assumes TABLELENGTH is a power of 2 */
# define TABLELENGTH (1<<13)
char **addresses=NULL; /* where to find the string */
int *codes=NULL, /* the code value */
*linktonext=NULL, /* the next index in the hash chain */
*lengths=NULL, /* the length of the coded string */
*codeindex=NULL; /* the index for a given code */
int nextcode; /* the next unused code */
/* hashit is supposed to give a unique fairly random number in the table for
each length a and string b */
# define hashit(a,b) (51*a+53*(57*b[0]+59*(61*b[a-1]+b[a>>1])))&(TABLELENGTH-1)
void compress(initcodesize,ptr,outfile,size)
int initcodesize; /* the initial compression bits */
char * ptr; /* where the data comes from */
FILE * outfile; /* where the output goes */
int size; /* how much data */
{int currentcode,prefixcode=0,codesize,maxbits=12,maxcode;
int clearcode,eoicode,currentplace=0,length,blocksize=0,bitoffset;
int findcode();
unsigned long outputword;
unsigned char blockbuffer[256]; /* to hold data blocks before writing */
/* allocate space for hash tables */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
if (NULL==(addresses=(char **) malloc(sizeof(char *)*TABLELENGTH)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentplace=size;
}
/* set up initial code table */
inittable(initcodesize);
clearcode=(1<<initcodesize);
eoicode=clearcode+1;
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
/* start with a clear code */
outputword=clearcode;
bitoffset=codesize;
/* now do the compressing */
while (currentplace<size)
{ /* check if codesize needs increasing */
if (nextcode>maxcode)
{codesize++;
maxcode=1<<codesize;
/* if too big, then reset compressor */
if (codesize>maxbits)
{if (bitoffset) outputword|=(clearcode<<bitoffset);
else outputword=clearcode;
bitoffset+=maxbits;
inittable(initcodesize);
codesize=initcodesize+1;
maxcode=1<<codesize;
nextcode=eoicode+1;
}
}
/* look for an unstored string */
length=1;
while (nextcode>
(currentcode=findcode(length,(char *)(ptr+currentplace))))
{prefixcode=currentcode;
length++;
if ((currentplace+length)>=size) break;
}
nextcode++;
currentplace+=(length-1);
/* output the prefix code */
if (bitoffset) outputword|=(prefixcode<<bitoffset);
else outputword=prefixcode;
bitoffset+=codesize;
/* output finished bytes to blocks */
while (bitoffset>=8)
{blockbuffer[blocksize]=outputword&0xff;
outputword>>=8;
bitoffset-=8;
blocksize++;
/* output filled block */
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
}
/* output the end of information code */
if (bitoffset) outputword|=(eoicode<<bitoffset);
else outputword=eoicode;
bitoffset+=codesize;
/* finish outputting the data */
while (bitoffset>=0)
{blockbuffer[blocksize]=(char) (outputword&0xff);
outputword>>=8;
bitoffset-=8;
blocksize++;
if (blocksize>=254)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
blocksize=0;
}
}
/* output the last block */
if (blocksize)
{fputc((char) blocksize, outfile);
fwrite(blockbuffer,blocksize,1,outfile);
}
/* a final zero block count */
fputc(0, outfile);
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void decompress(initcodesize,ptr,ptr1,size)
int initcodesize;
unsigned char *ptr, *ptr1; /* compressed data from ptr go to ptr1 */
int size; /* an upper limit purely as a check */
{int i,currentcode,codesize=0,maxbits=12,blocksize;
int clearcode,eoicode,codemask=0;
int bitoffset=0,indx,oldindx=0;
int currentplace=0,oldplace=0;
int findcode();
unsigned long inputword=0;
unsigned char *p1, *p2;
/* first deblock the data */
p1=p2=ptr;
blocksize=(*p1++);
while (blocksize)
{while (blocksize--)
(*p2++)=(*p1++); /* a wonderful example of how abstruse C can be */
blocksize=(*p1++);
}
/* set up initial code table */
currentcode=clearcode=(1<<initcodesize);
eoicode=clearcode+1;
/* allocate space for hash table */
if (NULL==(codes=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(linktonext=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(lengths=(int *) malloc(sizeof(int)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentcode=eoicode;
}
/* need one extra place in codeindex for overflow before resetting: */
if (NULL==(codeindex=(int *) malloc(sizeof(int)*4097)))
{fprintf(stderr,"compress: trouble allocating tables\n");
currentcode=eoicode;
}
if (NULL==(addresses=(char **) malloc(sizeof(char*)*TABLELENGTH)))
{fprintf(stderr,"decompress: trouble allocating tables\n");
currentplace=eoicode;
}
while (currentcode!=eoicode)
{if (currentcode==clearcode) /* reset the decompressor */
{inittable(initcodesize);
codesize=initcodesize+1;
nextcode=eoicode+1;
codemask=(1<<codesize)-1;
oldindx=(-1);
}
else /* code represents data */
{indx=codeindex[currentcode]; /* where in table is currentcode */
if (indx>=0) /* it is there */
{ /* put it into the output */
for (i=0;i<lengths[indx];i++)
ptr1[currentplace+i]=addresses[indx][i];
if (oldindx>=0) /* first character treated differently */
{findcode(lengths[oldindx]+1,(char *) (ptr1+oldplace));
nextcode++; /* add new code to table */
}
oldplace=currentplace;
currentplace+=lengths[indx];
oldindx=indx;
}
else /* not in table yet; must be old code plus last=first character */
{for (i=0;i<lengths[oldindx];i++)
ptr1[currentplace+i]=addresses[oldindx][i];
ptr1[currentplace+lengths[oldindx]]=addresses[oldindx][0];
/* store new code */
findcode(lengths[oldindx]+1,(char *)(ptr1+currentplace));
oldplace=currentplace;
currentplace+=lengths[oldindx]+1;
oldindx=codeindex[nextcode++];
}
/* crude error checking */
if ((oldindx<0)||(currentplace>size))
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* check if codesize needs increasing */
if (nextcode>codemask)
if (codesize<maxbits)
{codesize++;
codemask=(1<<codesize)-1;
}
while (bitoffset<codesize) /* read some more data */
{if (bitoffset) inputword|=(((int)(*ptr++))<<bitoffset);
else inputword=(*ptr++);
bitoffset+=8;
}
/* strip off current code */
currentcode=inputword&codemask;
inputword>>=codesize;
bitoffset-=codesize;
if (currentcode>nextcode)
{fprintf(stderr,"gif file appears to be corrupt\n");
break;
}
}
/* deallocate tables */
if (NULL!=codes) free((char *) codes);
if (NULL!=linktonext) free((char *) linktonext);
if (NULL!=lengths) free((char *) lengths);
if (NULL!=codeindex) free((char *) codeindex);
if (NULL!=addresses) free((char *) addresses);
codes=linktonext=lengths=codeindex=NULL;
addresses=(char **) NULL;
return;
}
void inittable(size)
int size;
{int i,findcode();
for (i=0;i<TABLELENGTH;i++)
{linktonext[i]=(-1);
codes[i]=(-1);
lengths[i]=(-1);
}
for (i=0;i<4096;i++)
codeindex[i]=(-1);
/* store initial codes for raw characters */
nextcode=0;
for (i=0;i<(1<<size);i++)
{stringbuffer[i]=i;
findcode(1,(char *) (stringbuffer+i));
nextcode++;
}
return;
}
int findcode(length,string)
char *string;
int length;
/* return code for string of given length;
if not found, store it and return nextcode */
{int i,j,indx,previousindex;
indx=hashit(length,string);
/* look for string in table */
previousindex=indx;
while (indx>0)
{if (lengths[indx]==length) /* is the length right ? */
{for (j=0;j<length;j++)
if ((string[j]) != (addresses[indx][j])) break;
if (j==length) return codes[indx];
}
previousindex=indx;
indx=linktonext[indx];
}
/* not found, so store it */
indx=previousindex;
i=indx;
while (codes[i]>=0) /* find an unused slot in table */
{i++;
if (i>=TABLELENGTH) i-=TABLELENGTH;
}
if (i!=indx)
{linktonext[indx]=i; /* link to it */
indx=i; /* move to it */
}
codes[indx]=nextcode; /* save the new code */
lengths[indx]=length;
addresses[indx]=string;
codeindex[nextcode]=indx;
return nextcode;
}
Raw
xsand.txt
XSand
by
Michael Creutz
[email protected]
This program simulates the sandpile automaton of Bak, Tang, and
Wiesenfeld. This is a simple model originally presented to illustrate
the concept of self organized criticality. This program generated the
pictures on page 46 of Briggs' book "Fractals -- The Patterns of
Chaos." Additonal pictures from this program appear in my article
with Per Bak appearing in "Fractals and Disordered Systems"
(Springer-Verlag, 1994).
The updating rule is extremely simple. On each cell of an initially
198 by 198 lattice is an integer amount of sand between 0 to 7,
inclusive. If this value exceeds 3, that cell is regarded as
"unstable" and for the next time step it takes four of its grains of
sand and places one on each of its neighbors. The updating is
simultaneous for all cells. The total amount of sand is conserved
except at the boundary.
The basic idea of self organized criticality is that after lots of
random addition of sand followed by relaxation, the system will
automatically enter a critical state where the size of an avalanche
created by additional sand addition is unpredictable without actually
running the process. The sizes of the ensuing avalanches
statistically have a power law distribution without any characteristic
scale.
Feel free to distribute these files to anyone interested. The latest
version can be found at
"http://thy.phy.bnl.gov/www/xtoys/xtoys.html".
A version for the Amiga can also be found there.
GADGETS
The screen opened by the program has lots of gadgets on it. They are:
Pause: Does the obvious. When paused, a click in the main window
outside any other gadget does a single updating step.
Cell size: Makes bigger "pixels". The overall lattice size can be
adjusted by resizing the window. The size in the x direction is
rounded to two less than a multiple of the number of bytes in a long
word. The initial 198 by 198 is the minimum allowed so the buttons
will fit nicely. Resizing fills the lattice with the selected pen
color.
Fill: Fills the lattice with the current pen color (height).
Tracer: When this is active those sites which "tumble" by being stable
are flagged and colored "cornflower blue" when they become stable. In
this way you can trace where an avalanche has gone.
Clear Trace: This clears the flagged sites to their normal color.
Double: This doubles all heights modulo 8. This is a convenient way
to quickly add lots of sand.
Auto-D: When this is active, the system will automatically double in
height whenever all sites become stable.
A color bar. This shows the colors by which the various heights are
shown. These squares are also gadgets, and by clicking on one of them
the corresponding color becomes the current pen color, shown in the
small box near the color bar. If any mouse button is pressed
while over the lattice, you can sketch with this color. To start an
avalanche, just select a color larger than 3.
The save button saves the current configuration as "xsand.gif," a
standard gif file that you can print or manipulate with any graphics
program that likes gif files. The restore gadget will reload a
previously stored configuration. Note that for the save feature to
work you need write permission on the directory from which you are
running. The "+ saved" button adds the saved configuration to the
presently displayed one, modulo 8 on each cell.
The other gadgets let you change the boundary conditions. The program
starts with open boundaries and sand is lost for any tumbling on the
edge. Periodic boundaries connect the left and right sides as well as
the top and bottom. In this case sand is never lost and it may be
that avalanches never end. Sandy boundaries represent an infinite
source of sand just outside the lattice. Every step a grain is tossed
on each edge site and two on each corner site. Try the flow boundary
condition to see what it does.
SOME EXPERIMENTS
1. With the mouse, scribble some sand randomly on the system. Then
repeatedly hit the "d" key to fill the lattice with a random mess.
Wait a few minutes until the system stabilizes and activity ceases.
Now you should be in the critical ensemble. Make the active color 4,
and click the mouse over the lattice. This will start an avalanche,
which unpredictably might be large or small.
After a few avalanches, turn on the trace button and make some more.
Now you can follow where the avalanche has passed.
Note that the avalanche regions always wind up simply connected, with
no untumbled islands left over. This is a theorem, and is true for
any state in the critical ensemble, but not an arbitrary state.
Select height 0 or 1 for the pen, scribble over a small region, and
then go back to making some more avalanches. Now it should be easy to
make islands, because by removing some sand you have left the critical
ensemble.
2. Select height 0 and clear the system. Then make the boundaries
sandy until the system fills up and nothing changes any more. Switch
back to open boundaries to let the excess sand run off. Try doubling
the final state and letting it relax back. Note that it returns
exactly to itself. This state is the unique one in the critical
ensemble with this property. Indeed, group theoretically this state
represents the identity.
3. Clear the system. Run one or two steps with sandy boundaries, and
then go back to open boundaries. Toggle the autodouble button on and
watch the show. This also yields the identity state.
4. Fill the system with height 2. Now draw a picture with height 3.
Put the boundaries sandy for a few moments, and then open them up
again. Note how the picture is restored after the avalanche ends.
This is a property of any state in the critical ensemble. Indeed,
this is a way to test that a state is critical. Try drawing some more
with height 0 or 1. Depending on the picture, the above experiment
may or may not mess up your picture.
5. Save the identity from experiment 2. Draw a picture using only
heights 2 and 3, as in the previous experiment. Use the "+ saved"
button to add in the identity. After a while your picture should
magically reappear.
6. Fill the system with height 2 and go to periodic boundaries.
Sketch a bit with height 4. With enough sketching, the avalanches
will no longer stop. The resulting dynamics can be hypnotic.
REFERENCES
P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987).
D. Dhar, Phys. Rev. Lett. 64, 1613 (1990).
M. Creutz, Computers in Physics 5, 198 (1991).
Raw
xwaves.c
/* Michael Creutz */
/* [email protected] */
/* to compile: */
/* cc -O xwaves.c -lX11 -lm */
/* January 2000: corrected array bounds on xpoints */
/* August 2006: button errors fixed */
# include <X11/Xlib.h>
# include <X11/Xutil.h>
# include <X11/Xos.h>
# include <X11/Xatom.h>
# include <X11/cursorfont.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <malloc.h>
# include <unistd.h>
# include <math.h>
#include <sys/time.h>
#include <sys/types.h>
#include <unistd.h>
struct timeval timeout;
/* some colors to use */
#define WHITE 0
#define BLACK 1
#define RED 2
#define BLUE 3
#define GREEN 4
#define YELLOW 5
#define MAGENTA 6
#define CYAN 7
#define ORANGE 8
#define PURPLE 9
#define PINK 10
#define BROWN 11
#define GREY 12
#define TURQUOISE 13
#define GOLD 14
#define NAVY 15
#define TAN 16
#define VIOLET 17
/* NSITES = number of lattice sites */
# define NSITES 256
double phi[NSITES],cpi[NSITES]; /* position space field and momentum */
double fphi[NSITES], fcpi[NSITES]; /* fourier space momentum */
double cn[NSITES],sn[NSITES]; /* for cosine and sin of angles */
double omega[NSITES]; /* frequency for given component */
double comega[NSITES],somega[NSITES];/* cos and sin of dt*omega */
double pi; /* 3.14.... */
int xcenter,ycenter,vx,vy; /* center of drifting blob */
XPoint xpoints[NSITES+7]; /* for drawing polygons */
# define dt (.02*(1-.8*slow))
void drawbutton(),openwindow(),makebuttons(),update(),repaint(),
cleanup(),drawit(),drawwave(),drawblob(),fastfourier(),
xtof(),ftox(),makedispersion();
/* dimensions for control boxes */
# define BUTTONWIDTH 72
# define BUTTONHEIGHT 18
/* initial button settings */
int paused=0; /* is the system paused? */
int damped=0; /* is the system damped? */
int reversed=0; /* are we going backwards? */
int mass=1; /* does the wave have mass? (0 or 1 only) */
int blob=0; /* display as blob? */
int drift=0; /* let blob drift */
int slow=0; /* run at slow speed */
int equation=0; /* 0 for light, 1 for mesons, 2 for water */
# define WATER 2
char * equtext[3]={"light","mesons","water"};
char stringbuffer[256]; /* generally useful */
/* to simplify typing */
# define line(x1,y1,x2,y2) XDrawLine(display,window,gc1,x1,y1,x2,y2)
# define line2(x1,y1,x2,y2) XDrawLine(display,window,gc2,x1,y1,x2,y2)
# define color(i) XSetForeground(display,gc1,translate[i])
# define color2(i) XSetForeground(display,gc2,translate[i])
/* various window stuff */
Display *display;
int screen;
static char *progname;
Window window,quitbutton,pausebutton,vacuumbutton,rightbutton,dampbutton,
reversebutton,blobbutton,driftbutton,slowbutton,packetbutton,
equbutton,makebutton();
XColor xcolor,colorcell;
Colormap cmap;
GC gc1,gc2;
int windowwidth,windowheight,playwidth,playheight;
XFontStruct *font=NULL;
int font_height,font_width;
XSizeHints size_hints;
int darkcolor,lightcolor,black,white;
long translate[256]; /* for converting colors */
double norm=1/(1.*NSITES);
int main(int argc, char **argv)
{int i,j,k;
double phipeg,w,theta,theta0;
unsigned int width, height;
XEvent report;
progname=argv[0];
pi=4*atan(1.0);
for (i=0;i<NSITES;i++)
{cn[i]=cos(2*pi*i/(1.*NSITES));
sn[i]=sin(2*pi*i/(1.*NSITES));
}
/* initial configuration */
for (i=0;i<NSITES;i++)
fphi[i]=fcpi[i]=0.;
fphi[15]=NSITES/8.;
fphi[17]=NSITES/8.;
ftox();
makedispersion();
openwindow(argc,argv);
makebuttons();
xcenter=playwidth/2;
ycenter=playheight/2;
vx=0;
vy=0;
/* loop forever, looking for events */
while(1)
{if ((0==XPending(display))&&(0==paused))
update();
else
{XNextEvent(display,&report);
switch (report.type)
{case Expose:
if ((report.xexpose.window)!=window) break; /* cuts down flashing,
but you might remove this line if things aren't being redrawn */
if (report.xexpose.count!=0) break; /* more in queue, wait for them */
repaint();
break;
case ConfigureNotify:
width=report.xconfigure.width;
height=report.xconfigure.height;
if ((width<size_hints.min_width)||(height<size_hints.min_height))
{fprintf(stderr,"%s: window too small to proceed.\n",progname);
cleanup();
}
else if ((width!=windowwidth)||(height!=windowheight))
{windowwidth=width;
windowheight=height;
playwidth=windowwidth;
playheight=windowheight-3*BUTTONHEIGHT-1;
makebuttons();
xcenter=playwidth/2;
ycenter=playheight/2;
}
break;
case ButtonPress:
if (report.xbutton.window==quitbutton)
cleanup();
else if (report.xbutton.window==pausebutton)
{paused=1-paused;
drawbutton(pausebutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"pause",1-2*paused);
}
else if (report.xbutton.window==dampbutton)
{damped=1-damped;
drawbutton(dampbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"damp",1-2*damped);
}
else if (report.xbutton.window==blobbutton)
{blob=1-blob;
drawbutton(blobbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"blob",1-2*blob);
drawit();
}
else if (report.xbutton.window==rightbutton)
{
for (i=NSITES/2;i<NSITES;i++)
{fphi[i]=0;
fcpi[i]=0;
}
fphi[0]=fcpi[0]=0;
ftox();
drawit();
}
else if (report.xbutton.window==vacuumbutton)
{
for (i=0;i<NSITES;i++)
phi[i]=cpi[i]=0.0;
xtof();
drawit();
}
else if (report.xbutton.window==reversebutton)
{reversed=1-reversed;
vx=-vx;
vy=-vy;
for (i=0;i<NSITES;i++)
{cpi[i]*=-1;
}
xtof();
drawbutton(reversebutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"reverse",1-2*reversed);
drawit();
}
else if (report.xbutton.window==driftbutton)
{drift=1-drift;
if (drift)
{vx=1;
vy=2;
}
else
{vx=vy=0;
xcenter=playwidth/2;
ycenter=playheight/2;
}
drawbutton(driftbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"drift",1-2*drift);
drawit();
}
else if (report.xbutton.window==slowbutton)
{slow=1-slow;
drawbutton(slowbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
"slow",1-2*slow);
makedispersion();
}
else if (report.xbutton.window==packetbutton)
{
for (i=0;i<NSITES;i++)
fphi[i]=fcpi[i]=0.;
fphi[15]=NSITES/8.;
fphi[17]=NSITES/8.;
ftox();
drawit();
}
else if (report.xbutton.window==equbutton)
{equation=report.xbutton.x/BUTTONWIDTH;
for (i=0;i<3;i++)
drawbutton(equbutton,i*BUTTONWIDTH,0,BUTTONWIDTH,BUTTONHEIGHT,
equtext[i],1-2*(equation==i));
makedispersion();
}
else /* a mouse click in main window */
{i=report.xbutton.x;
j=report.xbutton.y;
if (j>playheight) {
update(); /* outside of playfield */
}
else /* perturb field */
{if (blob) /* do some trig */
{phipeg=(i-xcenter)*(i-xcenter)
+(j-ycenter)*(j-ycenter);
phipeg=4*sqrt(phipeg)/(1.*playheight)-.60;
if (j!=ycenter)
theta0=atan((xcenter-i)/(1.0*(ycenter-j)));
else
theta0=pi/2.0;
if (j>ycenter) /* put theta0 between 0 and 2 pi */
theta0+=pi;
if (theta0<0.) theta0+=2*pi;
}
else /* not displayed as a blob */
{theta0=2*pi*i/(1.*playwidth);
phipeg=(j-playheight/2)/(0.5*playheight);
}
/* now construct the perturbed field */
for (k=0;k<NSITES;k++)
{theta=2*pi*k/(NSITES-1);
w=sin(.5*(theta0-theta));
w*=(w*100);
w=1/(1+w);
phi[k]=phipeg*w+phi[k]*(1-w);
}
xtof();
drawit();
}
}
break;
default:
break;
} /* end of switch */
} /* end of if XPending */
} /* end of while(1) */
} /* end of main */
void update()
{int i;
double temp;
/* a kludgey replacement for usleep(): */
/* if (slow) */
{timeout.tv_sec=0;
timeout.tv_usec=50000;
select(0,NULL,NULL,NULL,&timeout);
}
if (damped)
for (i=0;i<NSITES;i++)
fcpi[i]*=0.98;
for (i=0;i<NSITES;i++)
{temp =fcpi[i]*comega[i]+fphi[i]*somega[i];
fphi[i]=fphi[i]*comega[i]-fcpi[i]*somega[i];
fcpi[i]=temp;
}
ftox();
drawit();
return;
}
void drawit()
{if (blob) drawblob();
else drawwave();
return;
}
void drawwave()
{int i,iy;
color(CYAN);
color2(MAGENTA);
for (i=0;i<NSITES;i++)
{iy=(playheight/2)*(1+phi[i]);
if (iy>=playheight) iy=playheight-1;
else if (iy<0) iy=0;
xpoints[i].x=(i*playwidth)/(NSITES-1);
xpoints[i].y=iy;
}
xpoints[NSITES].x=playwidth;
xpoints[NSITES].y=playheight-1;
xpoints[NSITES+1].x=0;
xpoints[NSITES+1].y=playheight-1;
XFillPolygon(display,window,gc1,xpoints,NSITES+2,Nonconvex,CoordModeOrigin);
xpoints[NSITES].y=0;
xpoints[NSITES+1].y=0;
XFillPolygon(display,window,gc2,xpoints,NSITES+2,Nonconvex,CoordModeOrigin);
return;
}
void drawblob()
{int i,ix,iy;
double r;
color(CYAN);
color2(MAGENTA);
xcenter+=vx;
ycenter+=vy;
for (i=0;i<NSITES;i++)
{r=.15+.25*phi[i];
iy=ycenter-playheight*r*cn[i];
ix= xcenter-playheight*r*sn[i];
if (iy>(playheight-1)) iy=playheight-1, vy=-abs(vy);
if (iy< 0 ) iy=0 , vy= abs(vy);
if (ix>playwidth ) ix=playwidth , vx=-abs(vx);
if (ix<0 ) ix=0 , vx= abs(vx);
xpoints[i].x=ix;
xpoints[i].y=iy;
}
/* draw inside, allowing for boundary if imposed */
XFillPolygon(display,window,gc2,xpoints,NSITES,Nonconvex,
CoordModeOrigin);
/* draw outside of region */
xpoints[NSITES].x=xpoints[0].x;
xpoints[NSITES].y=xpoints[0].y;
xpoints[NSITES+1].x=xcenter;
xpoints[NSITES+1].y=0;
xpoints[NSITES+2].x=playwidth;
xpoints[NSITES+2].y=0;
xpoints[NSITES+3].x=playwidth;
xpoints[NSITES+3].y=playheight-1;
xpoints[NSITES+4].x=0;
xpoints[NSITES+4].y=playheight-1;
xpoints[NSITES+5].x=0;
xpoints[NSITES+5].y=0;
xpoints[NSITES+6].x=xcenter;
xpoints[NSITES+6].y=0;
XFillPolygon(display,window,gc1,xpoints,NSITES+7,Nonconvex,
CoordModeOrigin);
return;
}
void repaint()
/* this fixes the window up whenever it is uncovered */
{int i;
color(BROWN);
XFillRectangle(display,window,gc1,0,playheight,
playwidth,windowheight-playheight);
drawbutton(pausebutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"pause",1-2*paused);
drawbutton(dampbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"damp",1-2*damped);
drawbutton(quitbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"quit",1);
drawbutton(vacuumbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"vacuum",1);
drawbutton(reversebutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"reverse",1-2*reversed);
drawbutton(blobbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"blob",1-2*blob);
drawbutton(rightbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"right",1);
drawbutton(driftbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"drift",1-2*drift);
drawbutton(slowbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"slow",1-2*slow);
drawbutton(packetbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"packet",1);
for (i=0;i<3;i++)
drawbutton(equbutton,i*BUTTONWIDTH,0,BUTTONWIDTH,BUTTONHEIGHT,
equtext[i],1-2*(equation==i));
color(BLACK);
line(0,playheight,windowwidth,playheight);
drawit();
return;
}
void openwindow(int argc, char **argv)
/* a lot of this is taken from basicwin in the Xlib Programming Manual */
{char *window_name="Wave tank";
char *icon_name="Wave";
long event_mask;
Pixmap icon_pixmap;
char *display_name=NULL;
int i;
# define icon_bitmap_width 16
# define icon_bitmap_height 16
static unsigned char icon_bitmap_bits[] = {
0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
/* open up the display */
if ((display=XOpenDisplay(display_name))==NULL)
{fprintf(stderr,"%s: cannot connect to X server %s\n",
progname,XDisplayName(display_name));
exit(-1);
}
screen=DefaultScreen(display);
cmap=DefaultColormap(display,screen);
darkcolor=black=BlackPixel(display,screen);
lightcolor=white=WhitePixel(display,screen);
translate[WHITE]=white;
translate[BLACK]=black;
translate[2]=lightcolor;
translate[3]=darkcolor;
for(i=4;i<256;i++)
translate[i]=translate[i%4];
if (DefaultDepth(display,screen)>1)
{
if (XAllocNamedColor(display,cmap,"salmon",&colorcell,&xcolor))
darkcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"wheat",&colorcell,&xcolor))
lightcolor=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"red",&colorcell,&xcolor))
translate[RED]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"blue",&colorcell,&xcolor))
translate[BLUE]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"green",&colorcell,&xcolor))
translate[GREEN]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"cyan",&colorcell,&xcolor))
translate[CYAN]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"orange",&colorcell,&xcolor))
translate[ORANGE]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"purple",&colorcell,&xcolor))
translate[PURPLE]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"yellow",&colorcell,&xcolor))
translate[YELLOW]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"pink",&colorcell,&xcolor))
translate[PINK]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"brown",&colorcell,&xcolor))
translate[BROWN]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"grey",&colorcell,&xcolor))
translate[GREY]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"turquoise",&colorcell,&xcolor))
translate[TURQUOISE]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"gold",&colorcell,&xcolor))
translate[GOLD]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"magenta",&colorcell,&xcolor))
translate[MAGENTA]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"navy",&colorcell,&xcolor))
translate[NAVY]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"tan",&colorcell,&xcolor))
translate[TAN]=colorcell.pixel;
if (XAllocNamedColor(display,cmap,"violet",&colorcell,&xcolor))
translate[VIOLET]=colorcell.pixel;
}
/* feel free to type in more colors, I got bored */
for(i=18;i<256;i++)
translate[i]=translate[i%18];
/* make the main window */
windowwidth=6*BUTTONWIDTH;
windowheight=600;
playwidth=windowwidth;
playheight=windowheight-3*BUTTONHEIGHT-1;
window=XCreateSimpleWindow(display,RootWindow(display,screen),
0,0,windowwidth,windowheight,4,translate[14],lightcolor);
/* make the icon */
icon_pixmap=XCreateBitmapFromData(display,window,
icon_bitmap_bits,icon_bitmap_width,icon_bitmap_height);
size_hints.flags=PPosition | PSize | PMinSize;
size_hints.min_width=windowwidth;
size_hints.min_height=100;
#ifdef X11R3
size_hints.x=x;
size_hints.y=y;
size_hints.width=windowwidth;
size_hints.height=windowheight;
XSetStandardProperties(display,window,window_name,icon_name,
icon_pixmap,argv,argc,&size_hints);
#else
{XWMHints wm_hints;
XClassHint class_hints;
XTextProperty windowName, iconName;
if (XStringListToTextProperty(&window_name,1,&windowName)==0)
{fprintf(stderr,"%s: structure allocation for windowName failed.\n"
,progname);
exit(-1);
}
if (XStringListToTextProperty(&icon_name,1,&iconName)==0)
{fprintf(stderr,"%s: structure allocation for iconName failed.\n"
,progname);
exit(-1);
}
wm_hints.initial_state=NormalState;
wm_hints.input=True;
wm_hints.icon_pixmap=icon_pixmap;
wm_hints.flags=StateHint|IconPixmapHint|InputHint;
class_hints.res_name=progname;
class_hints.res_class="Basicwin";
XSetWMProperties(display,window,&windowName,&iconName,
argv,argc,&size_hints,&wm_hints,&class_hints);
}
#endif
/* pick the events to look for */
event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
XSelectInput(display,window,event_mask);
/* pick font: 9x15 is supposed to almost always be there */
if ((font=XLoadQueryFont(display,"9x15"))==NULL)
if ((font=XLoadQueryFont(display,"fixed"))==NULL)
{fprintf(stderr,"%s: Sorry, having font problems.\n",progname);
exit(-1);
}
font_height=font->ascent+font->descent;
font_width=font->max_bounds.width;
/* make graphics context: */
gc1=XCreateGC(display,window,0,NULL);
XSetFont(display,gc1,font->fid);
XSetForeground(display,gc1,black);
XSetBackground(display,gc1,lightcolor);
gc2=XCreateGC(display,window,0,NULL);
XSetFont(display,gc2,font->fid);
XSetForeground(display,gc2,black);
XSetBackground(display,gc2,lightcolor);
/* show the window */
XMapWindow(display,window);
return;
}
void makebuttons()
{
/* first destroy any old buttons */
XDestroySubwindows(display,window);
/* now make the new buttons */
quitbutton=makebutton(0,windowheight-BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);
vacuumbutton=makebutton(3*BUTTONWIDTH,windowheight-BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
pausebutton=makebutton(BUTTONWIDTH,windowheight-BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
dampbutton=makebutton(2*BUTTONWIDTH,windowheight-BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
reversebutton=makebutton(windowwidth-2*BUTTONWIDTH,windowheight-3*BUTTONHEIGHT
,BUTTONWIDTH,BUTTONHEIGHT);
blobbutton=makebutton(windowwidth-BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
rightbutton=makebutton(windowwidth-2*BUTTONWIDTH,windowheight-2*BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
driftbutton=makebutton(windowwidth-BUTTONWIDTH,windowheight-2*BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
slowbutton=makebutton(windowwidth-BUTTONWIDTH,windowheight-BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
equbutton=makebutton(0,windowheight-3*BUTTONHEIGHT,
3*BUTTONWIDTH,BUTTONHEIGHT);
packetbutton=makebutton(windowwidth-2*BUTTONWIDTH,windowheight-BUTTONHEIGHT,
BUTTONWIDTH,BUTTONHEIGHT);
return;
}
void cleanup()
{XUnloadFont(display,font->fid);
XFreeGC(display,gc1);
XCloseDisplay(display);
exit(1);
}
Window makebutton(int xoffset, int yoffset, int xsize, int ysize)
/* Puts button of specified dimensions on window. Nothing is drawn. */
{Window buttonwindow;
long event_mask;
buttonwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
xsize,ysize,0,black,lightcolor);
event_mask=ButtonPressMask|ExposureMask; /* look for mouse-button presses */
XSelectInput(display,buttonwindow,event_mask);
XMapWindow(display,buttonwindow);
return buttonwindow;
}
void drawbutton(Window buttonwindow, int xoffset,
int yoffset, int xsize, int ysize, char *text, int state)
/* Draw a button in buttonwindow of specified dimensions with text
centered. Color of button determined by sign of "state."
size of border determined by magnitude. */
{int textlength,i,j;
int cdark,clight,cup,cdown;
int cleft,cright,cbutton,ctext;
cup=lightcolor;
cdown=darkcolor;
cdark=black;
clight=white;
if (state<0)
{cbutton=cdown;
ctext=clight;
cleft=cdark;
cright=clight;
}
else
{cbutton=cup;
ctext=cdark;
cleft=clight;
cright=cdark;
}
j=abs(state);
XSetForeground(display,gc1,cbutton);
XFillRectangle(display,buttonwindow,gc1,xoffset+j,yoffset+j,
xsize-2*j,ysize-2*j);
XSetForeground(display,gc1,cleft);
XFillRectangle(display,buttonwindow,gc1,xoffset,yoffset,xsize,j);
XFillRectangle(display,buttonwindow,gc1,xoffset,yoffset,j,ysize);
XSetForeground(display,gc1,cright);
for (i=0;i<j;i++)
{XDrawLine(display,buttonwindow,gc1,
xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
XDrawLine(display,buttonwindow,gc1,
xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
}
if (NULL!=text)
{textlength=strlen(text);
XSetForeground(display,gc1,ctext);
XDrawString(display,buttonwindow,gc1,
xoffset+(xsize-textlength*font_width)/2,
yoffset+(ysize+font_height)/2-2,
text,textlength);
}
XSetForeground(display,gc1,black);
return;
}
void xtof()
{fastfourier(phi,cpi,fphi,fcpi,NSITES,1,1,1);
return;
}
void ftox()
{int i;
fastfourier(fphi,fcpi,phi,cpi,NSITES,1,1,-1);
for (i=0;i<NSITES;i++)
{phi[i]*=norm;
cpi[i]*=norm;
}
return;
}
void makedispersion()
/* set up dispersion relations for the models */
{int j;
double k;
for (j=0;j<NSITES;j++)
{if (2*j<NSITES) k=j;
else k=NSITES-j;
if (equation==WATER)
omega[j]=4*sqrt(k);
else
omega[j]=sqrt(k*k+NSITES*NSITES*equation/200);
comega[j]=cos(dt*omega[j]);
somega[j]=sin(dt*omega[j]);
}
return;
}
void fastfourier(double *rex,double *imx, double *ref,double *imf,
int n,int d1,int d2,int direction)
/* recursively calls itself on two halves of the vector */
/* (rex[],imx[]) is the input vector
(ref[],imf[]) is the output
n=vector length
d1=spacing between elements of input vector
d2= same for output
direction= +/-1 for forward and reverse fourier transforms
vector is not normalized */
{int j,Ndn,nd2,i1,i2;
double c,s,temp1,temp2,temp3;
if (n==2) /* do two elements case by hand */
{*ref=*rex+*(rex+d1);
*imf=*imx+*(imx+d1);
*(ref+d2)=*rex-*(rex+d1);
*(imf+d2)=*imx-*(imx+d1);
return;
}
nd2=n>>1; /* half the vector length */
/* recursively call fastfourier */
fastfourier(rex ,imx ,ref ,imf ,nd2,2*d1,d2,direction);
fastfourier(rex+d1,imx+d1,ref+nd2*d2,imf+nd2*d2,nd2,2*d1,d2,direction);
/* combine the results */
Ndn=NSITES/n;
for (j=0;j<nd2;j++)
{c=cn[j*Ndn];
s=direction*sn[j*Ndn];
i1=d2*j;
i2=d2*(j+nd2);
temp1 =ref[i1]+c*ref[i2]-s*imf[i2];
temp2 =imf[i1]+c*imf[i2]+s*ref[i2];
temp3 =ref[i1]-c*ref[i2]+s*imf[i2];
imf[i2]=imf[i1]-c*imf[i2]-s*ref[i2];
ref[i1]=temp1;
imf[i1]=temp2;
ref[i2]=temp3;
}
return;
}
Raw
xwaves.txt
Xwaves
Michael Creutz
[email protected]
This program illustrates the behaviors of three common wave equations,
differing in the relation between the frequency w and the wave number
k. First is the simple wave equation with w proportional to k; this
is the equation is obeyed by free photons. Second is the Klein-Gordon
equation with w proportional to the sqrt(k^2+m^2) with m a constant
representing a particle mass, say that of the electron or a pion.
Finally is the equation for long wavelength water waves, which obey w
proportional to sqrt(k).
This should run on any X platform. To compile, try the xtoys
Makefile. Otherwise try
cc -O -o xwaves xwaves.c -lm -lX11
If this does not work immediately, try adding -I and -L flags to
directories where the X includes and libraries are located. For
example, on a Sun try
cc -O -I/usr/openwin/include -L/usr/openwin/lib -o xwaves xwaves.c -lm -lX11
The program starts with wave packets cycling around the periodic
boundary conditions. The equation is that for light, and the wave
packets move with the same velocity as the short waves of which they
are composed.
If you now click on the "mesons" button, the mass term of the
Klein-Gordon equation is turned on. Note how the short components now
move faster than the packets. This illustrates how the group velocity
of a wave equation can differ from the phase velocity. For the
Klein-Gordon equation, the phase velocity exceeds the speed of light,
but this is not a problem because information is carried only by wave
packets, which move at the group velocity. The latter is always less
than the speed of light.
Selecting the water button shows a similar effect. In this case the
group velocity is exactly one half the phase velocity. This follows
from the frequency going as the square root of the wavelength. This
is my favorite example of an interesting result following from
dimensional analysis.
Mathematically, the phase velocity v_p of a wave equation is given by
the relation
v_p=w/k
while the group velocity v_g is
v_g=dw/dk
Only for the simple wave equation are these equal.
The program has a variety of other buttons to play with. The
``quit,'' ``pause,'' and ``slow'' buttons do the obvious things. A
mouse click in the wave portion of the window will insert a
disturbance, much like throwing a rock into a pond. The ``damp''
button slowly damps out any waves present, while ``vacuum'' settles
the system immediately. The ``reverse'' button reverses the motion of
all waves. The ``right'' button eliminates all left moving waves.
Pressing ``packet'' restores the initial two packet state.
The ``blob'' gadget changes the display so the waves run around a
spherical blob rather than along the periodic line. The ``drift''
button lets that blob drift around and bounce off the window walls.
It has no effect when not in the blob display mode. A highly excited
drifting blob with the slow button seems like a blend of Matisse with
a lava lamp.
For each one of the three equations, try calming the state with the
vacuum button and then exciting a pulse with a mouse click. Observe
how the light equation gives a pulse dividing into oppositely moving
pieces that move without changing shape. With the meson equation, an
oscillating central region emits at first short wavelengths followed
by longer ones. With water, the long waves are the first to leave the
excitation region.
Some technical details: The waves are displayed using the Xlib routine
XPolyFill() with 256 points on the wave surface. To avoid
discretization errors, the evolution is calculated in momentum space
and then a fast Fourier transform gives the position space
coordinates. You can easily obtain other wave equations by modifying
the routine makedispersion().
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment