/* program latgas */
#include <stdio.h>

void initial(int lat[],int *lx,int *ly,int *nstep,int *nplot);
/*    void initial_graphics(int lx,int ly);  */
void update(int lat[],int rule[],int nn[][6],int lx,int ly);
void nntable(int nn[][6],int lx,int ly);
void ruletable(int rule[]);
/*    void plot(int lat[],int lx,int ly);   */

/* note how arrays are initialized */
int nv[] = {3,4,5,0,1,2};       /* global variables available */
int mask[] = {1,2,4,8,16,32,64,128}; /* to all functions */

main()
/* simulates lattice gas fluid model
  graphics routines not listed */
{
  int lat[10000],nn[10000][6];
  int lx,ly,nstep,nplot,istep;
  int rule[1024];
  initial(lat,&lx,&ly,&nstep,&nplot);
/*    initial_graphics(lx,ly); */
  nntable(nn,lx,ly);
  ruletable(rule);
/*    plot(lat,lx,ly);    */
  for (istep = 1; istep <= nstep; istep++)
    {
     update(lat,rule,nn,lx,ly);
/*        if (istep % nplot==0)
             plot(lat,lx,ly);  */
    }
}

void initial(int lat[],int *lx,int *ly,int *nstep,int *nplot)
{
  int i,j,n;
  /* system size */
  *lx = 50;
  *ly = 20;
  /* simulation parameters */
  *nstep = 100;
  *nplot = 1;
  /* begin with no particles */
  for (n = 0; n < (*lx)*(*ly); n++)
    lat[n] = 0;
  /* fill block in center of lattice with 6 particles per site */
  for (j = -1; j <= 1; j++)
    for (i = -1; i <= 1; i++)
     {
        n = (j+(*ly)/2)*(*lx) + i + (*lx)/2;
        lat[n] = 63;
     }
}

void update(int lat[],int rule[],int nn[][6],int lx,int ly)
{
  int i,j,n,dir;
  int latn[10000];
  /* bounce back boundary conditions v goes to -v */
  extern int nv[];
  extern int mask[];
  /* move particles */
  for (j=0; j<ly; j++)
    for (i=0; i<lx; i++)
      {
         n = j*lx + i;
         for (dir=0; dir<6 ;dir++)
         if ((lat[nn[n][dir]] & mask[dir]) != 0)
           {
             if ((lat[n] & mask[7]) != 0)
             /* reflection */
             latn[nn[n][dir]] = latn[nn[n][dir]] | mask[nv[dir]];
         else
            /* particle moves from nearest neighbor */
            latn[n] = latn[n] | mask[dir];
          }
      }/* end for */

  /* collisions */
  for (j=0; j<ly; j++)
    for (i=0; i<lx; i++)
      {
        n = j*lx + i;
        if ((lat[n] & mask[7]) == 0)
         {
           lat[n] = rule[latn[n]];
           latn[n] = 0;
         }
      } /* end for */
}

void ruletable(int rule[])
/* 6-bit saturated deterministic rule */
{
  int i;
  for (i=0; i<256; i++)
    rule[i] = i;
  rule[21] = 42;
  rule[42] = 21;
  rule[9] = 36;
  rule[18] = 9;
  rule[36] = 18;
  rule[27] = 45;
  rule[45] = 54;
  rule[54] = 27;
  rule[19] = 37;
  rule[37] = 19;
  rule[50] = 41;
  rule[41] = 50;
  rule[22] = 13;
  rule[13] = 22;
  rule[26] = 44;
  rule[44] = 26;
  rule[11] = 38;
  rule[38] = 11;
  rule[25] = 52;
  rule[52] = 25;
}

void nntable(int nn[][6],int lx,int ly)
     /* nn[n][dir] gives the neighbor whose particle moving in
        direction dir will move to site n. */
{
  int i,j,ip,im,jp,jm,n;
  for (j=0; j < ly; j = j+2)
    for (i=0; i < lx; i++)
      {
        n = j*lx + i;
        ip = i + 1;
        if (ip > lx-1) ip = 0;
          im = i-1;
        if (im < 0) im = lx - 1;
          jp = j + 1;
        if (jp > ly-1) jp = 0;
          jm = j - 1;
        if (jm < 0) jm = ly - 1;
        nn[n][0] = j*lx + im;
        nn[n][1] = jp*lx + im;
        nn[n][2] = jp*lx + i;
        nn[n][3] = j*lx + ip;
        nn[n][4] = jm*lx + i;
        nn[n][5] = jm*lx + im;
      }
  for (j=1; j < ly;j = j+2)
    for (i=0; i < lx;i++)
      {
       n = j*lx + i;
       ip = i + 1;
       if (ip > lx-1) ip = 0;
       im = i - 1;
       if (im < 0) im = lx - 1;
       jp = j + 1;
       if (jp > ly-1) jp = 0;
       jm = j - 1;
       if (jm < 0) jm = ly - 1;
       nn[n][0] = j*lx + im;
       nn[n][1] = jp*lx + i;
       nn[n][2] = jp*lx + ip;
       nn[n][3] = j*lx + ip;
       nn[n][4] = jm*lx + ip;
       nn[n][5] = jm*lx + i;
      }
}