      PROGRAM ca1
*     one-dimensional Boolean cellular automata
      IMPLICIT NONE
      INTEGER IR_, L, update(0:7), site(0:1601), tmax
      DATA site/1602*0/
      CALL GWopen(IR_, 0)
      CALL setrule(update)
      CALL initial(site,L,tmax)
      CALL iterate(site,L,update,tmax)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE setrule(update)
      IMPLICIT NONE
      CHARACTER Stmp1_*80
      INTEGER IR_, bit0, bit1, bit2, i, rule, update(0:7)
      DATA rule/90/
      CALL GWinput(IR_, 'rule number = ', Stmp1_)
      IF(IR_ .NE. 0) READ(Stmp1_,*) rule
      WRITE(*, *) 'rule number = ', rule
      DO i = 7, 0, -1
         update(i) = rule/2**i              ! find binary representation
         rule = rule - update(i)*2**i
         bit2 = i/4
         bit1 = (i - 4*bit2)/2
         bit0 = i - 4*bit2 - 2*bit1
*        show possible neighborhoods
         WRITE(*, '(3(I1,1X),A,$)') bit2, bit1, bit0, '   '
      END DO
      WRITE(*,*)
      DO i = 7, 0, -1
         WRITE(*, '(2X,I1,$)') update(i)    ! print rules
         WRITE(*, '(A,$)') '      '
      END DO
      WRITE(*,*)
      END
C
      SUBROUTINE initial(site,L,tmax)
      IMPLICIT NONE
      INTEGER IR_, px, py, site(0:1601), L, tmax
      CALL GWcolor(IR_, 0, 2)
      CALL GWsize(IR_, 1, px,py)
      CALL GWindow(IR_, 0.0, REAL(py-1), REAL(px-1), 0.0)
      L = px/4
      IF(L .GE. 1601) L = 1600
      tmax = py/4
      site(L/2) = 1                         ! center site
      CALL GWsrect(IR_, REAL(4*L/2),   0.0,
     +                  REAL(4*L/2+4), 4.0, 15) ! each site 4 x 4 pixels
      END
C
      SUBROUTINE iterate(site,L,update,tmax)
      IMPLICIT NONE
      INTEGER IR_, i, index, sitenew(0:1601), L, update(0:7),
     +       t, site(0:1601), tmax, Itmp1_
*     update lattice
*     need to introduce additional array, sitenew, to temporarily
*     store values of newly updated sites
      DO t = 1, tmax
         DO i = 1, L
            index = 4*site(i-1) + 2*site(i) + site(i+1)
            sitenew(i) = update(index)
            IF(sitenew(i) .EQ. 1) 
     +         CALL GWsrect(IR_, REAL(i*4),   REAL(t*4),
     +                           REAL(i*4+4), REAL(t*4+4), 15)
         END DO
         DO Itmp1_ = 1, L
            site(Itmp1_) = sitenew(Itmp1_)
         END DO
         site(0) = site(L)                  ! periodic boundary conditions
         site(L+1) = site(1)
      END DO
      END
