PROGRAM lattice_gas
! simulation of particle diffusion in a lattice gas
DIM site(40,40),occup$(-1:0)
DIM x(1200),y(1200),x0(1200),y0(1200)
LIBRARY "csgraphics"
CALL initial(L,L_2,N,tmax,#2)
CALL lattice(site(,),x(),y(),x0(),y0(),L,N,occup$(),#1)
FOR t = 1 to tmax
    CALL move(site(,),x(),y(),x0(),y0(),L,N,occup$(),#1)
    CALL data(x(),y(),x0(),y0(),L,L_2,N,t,#2)
NEXT t
END

SUB initial(L,L_2,N,tmax,#2)
    RANDOMIZE
    LET L = 40               ! linear dimension of lattice
    LET L_2 = 0.5*L
    LET N = 500              ! number of particles
    LET tmax = 50            ! # Monte Carlo steps per particle
    SET BACKGROUND COLOR "black"
    OPEN #2: screen 0.52,1.0,0,1.0
    SET WINDOW 0,tmax + 1,0,L_2 + 1
    SET COLOR "white"
    PRINT "Plot of <R^2(t)> versus t"
    BOX LINES 0,tmax,0,L_2
END SUB

SUB lattice(site(,),x(),y(),x0(),y0(),L,N,occup$(),#1)
    OPEN #1: screen 0,0.5,0.2,1
    LET r = 0.5              ! size of box
    CALL compute_aspect_ratio(L+r,xwin,ywin)
    SET WINDOW -0.1*xwin,xwin,-0.1*ywin,ywin
    SET COLOR "blue"
    BOX CIRCLE 2,3,2,3
    FLOOD 2.5,2.5
    BOX KEEP 2,3,2,3 in occup$(-1)
    CLEAR
    SET COLOR "black"
    BOX CIRCLE 2,3,2,3
    FLOOD 2.5,2.5
    BOX KEEP 2,3,2,3 in occup$(0)
    CLEAR
    SET COLOR "white"
    PRINT "density = "; N/(L*L)
    SET COLOR "red"
    BOX LINES 1-r,L+r,1-r,L+r
    FOR iy = 1 to L          ! draw lattice sites 
        FOR ix = 1 to L
            LET site(ix,iy) = 0
            PLOT POINTS: ix,iy
        NEXT ix
    NEXT iy
    LET i = 0
    DO
       LET xadd = int(L*rnd) + 1
       LET yadd = int(L*rnd) + 1
       IF site(xadd,yadd) = 0 then
          LET i = i + 1      ! number of particles added
          LET site(xadd,yadd) = -1     ! site occupied
          BOX SHOW occup$(-1) at xadd - 0.5,yadd - 0.5
          LET x(i) = xadd
          LET y(i) = yadd
          LET x0(i) = x(i)   ! x-coordinate at t = 0
          LET y0(i) = y(i)
       END IF
    LOOP until i = N
END SUB

SUB move(site(,),x(),y(),x0(),y0(),L,N,occup$(),#1)
    WINDOW #1
    FOR particle = 1 to N
        LET i = int(N*rnd) + 1
        LET xtrial = x(i)
        LET ytrial = y(i)
        CALL direction(xtrial,ytrial,L)
        IF site(xtrial,ytrial) >= 0 then    ! unoccupied site
           LET site(x(i),y(i)) = 0
           BOX SHOW occup$(0) at x(i) - 0.5,y(i) - 0.5
           PLOT POINTS: x(i),y(i)
           LET x(i) = xtrial
           LET y(i) = ytrial
           LET site(x(i),y(i)) = -1    ! new site occupied
           BOX SHOW occup$(-1) at x(i) - 0.5,y(i) - 0.5
        END IF
    NEXT particle
    SET COLOR "red"
    BOX LINES 0.5,L+0.5,0.5,L+0.5
END SUB

SUB direction(xtemp,ytemp,L)
    ! choose random direction and use periodic boundary conditions
    LET dir = int(4*rnd) + 1
    SELECT CASE dir
    CASE 1
         LET xtemp = xtemp + 1
         IF xtemp > L then LET xtemp = 1
    CASE 2
         LET xtemp = xtemp - 1
         IF xtemp < 1 then LET xtemp = L
    CASE 3
         LET ytemp = ytemp + 1
         IF ytemp > L then LET ytemp = 1
    CASE 4
         LET ytemp = ytemp - 1
         IF ytemp < 1 then LET ytemp = L
    END SELECT
END SUB

SUB data(x(),y(),x0(),y0(),L,L_2,N,t,#2)
    LET R2 = 0
    FOR i = 1 to N
        LET dx = x(i) - x0(i)
        LET dy = y(i) - y0(i)
        CALL separation(dx,dy,L,L_2)
        LET R2 = R2 + dx*dx + dy*dy
    NEXT i
    LET R2bar = R2/N
    WINDOW #2
    PLOT t,R2bar
END SUB

SUB separation(dx,dy,L,L_2)
    ! use periodic boundary conditions to determine separation
    IF abs(dx) > L_2 then LET dx = dx - sgn(dx)*L
    IF abs(dy) > L_2 then LET dy = dy - sgn(dy)*L
END SUB