PROGRAM hd
! dynamics of system of hard disks
! program based in part on Fortran program of Allen and Tildesley
PUBLIC x(100),y(100),vx(100),vy(100)
PUBLIC collision_time(100),partner(100)
PUBLIC N,Lx,Ly,t,timebig
LIBRARY "csgraphics"
CALL initial(vsum,rho,area)
CALL set_up_windows(ncolor,#1,#2)
CALL kinetic_energy(ke,#1)
CALL show_positions(flag$,#2)
LET temperature = ke/N
LET flag$ = ""
LET collisions = 0                ! number of collisions
DO
   CALL minimum_collision_time(i,j,tij)
   ! move particles forward and reduce collision times by tij
   CALL move(tij)
   LET t = t + tij
   LET collisions = collisions + 1
   CALL show_positions(flag$,#2)
   !CALL show_disks(ncolor,flag$,#2)
   CALL contact(i,j,virial)       ! compute collision dynamics
   LET vsum = vsum + virial
   CALL show_output(t,collisions,temperature,vsum,rho,area,#1)
   ! reset collision list for relevant particles
   CALL reset_list(i,j)
LOOP until flag$ = "stop"
CALL save_config(#2)
END

SUB initial(vsum,rho,area)
    DECLARE PUBLIC x(),y(),vx(),vy()
    DECLARE PUBLIC N,Lx,Ly,t
    DECLARE PUBLIC collision_time(),partner(),timebig
    LET t = 0
    INPUT prompt "read file (f) or lattice start (l) = ": start$
    IF start$ = "f" or start$ = "F" then
       INPUT prompt "file name = ": file$
       OPEN #1: name file$, access input
       INPUT #1: N
       INPUT #1: Lx
       INPUT #1: Ly
       INPUT #1: heading$
       FOR i = 1 to N
           INPUT #1: x(i),y(i)
       NEXT i
       INPUT #1: heading$
       FOR i = 1 to N
           INPUT #1: vx(i),vy(i)
       NEXT i
       CLOSE #1
    ELSE IF start$ = "l" or start$ = "L" then
       RANDOMIZE
       INPUT prompt "N = ": N     ! choose N so that sqr(N) an integer
       INPUT prompt "Lx = ": Lx
       INPUT prompt "Ly = ": Ly
       INPUT prompt "vmax = ": vmax
       LET nx = sqr(N)
       LET ny = nx
       IF nx >= Lx or nx >= Ly then
          PRINT "box too small"
          STOP
       END IF
       LET ax = Lx/nx             ! "lattice" spacing
       LET ay = Ly/ny
       LET i = 0
       FOR col = 1 to nx
           FOR row = 1 to ny
               LET i = i + 1
               LET x(i) = (col - 0.5)*ax
               LET y(i) = (row - 0.5)*ay
               ! choose random positions and velocities
               LET vx(i) = (2*rnd - 1)*vmax
               LET vy(i) = (2*rnd - 1)*vmax
           NEXT row
       NEXT col
    END IF
    CLEAR
    CALL check_overlap            ! check if two disks overlap
    CALL check_momentum
    LET area = Lx*Ly
    LET rho = N/area
    LET timebig = 1.0e10
    LET vsum = 0                  ! virial sum
    FOR i = 1 to N
        LET partner(i) = N
    NEXT i
    LET collision_time(N) = timebig
    ! set up initial collision lists
    FOR i = 1 to N
        CALL uplist(i)
    NEXT i
END SUB

SUB check_momentum
    DECLARE PUBLIC vx(),vy(),N
    LET vxsum = 0
    LET vysum = 0
    FOR i = 1 to N
        LET vxsum = vxsum + vx(i)      ! total center of mass velocity (momentum)
        LET vysum = vysum + vy(i)
    NEXT i
    LET vxsum = vxsum/N
    LET vysum = vysum/N
    FOR i = 1 to N
        LET vx(i) = vx(i) - vxsum
        LET vy(i) = vy(i) - vysum
    NEXT i
END SUB

SUB minimum_collision_time(i,j,tij)
    DECLARE PUBLIC N,collision_time(),partner(),timebig
    ! locate minimum collision time
    LET tij = timebig
    FOR k = 1 to N
        IF collision_time(k) < tij then
           LET tij = collision_time(k)
           LET i = k
        END IF
    NEXT k
    LET j = partner(i)
END SUB

SUB move(tij)
    DECLARE PUBLIC x(),y(),vx(),vy()
    DECLARE PUBLIC N,collision_time()
    DECLARE PUBLIC Lx,Ly
    DECLARE DEF pbc
    FOR k = 1 to N
        LET collision_time(k) = collision_time(k) - tij
        LET x(k) = x(k) + vx(k)*tij
        LET y(k) = y(k) + vy(k)*tij
        LET x(k) = pbc(x(k),Lx)
        LET y(k) = pbc(y(k),Ly)
    NEXT k
END SUB

SUB reset_list(i,j)
    DECLARE PUBLIC N,partner()
    ! reset collision list for relevant particles
    FOR k = 1 to N
        LET test = partner(k)
        IF k = i or test = i or k = j or test = j then
           CALL uplist(k)
        END IF
    NEXT k
    CALL downlist(i)
    CALL downlist(j)
END SUB

DEF separation(ds,L)
    IF ds > 0.5*L then
       LET separation = ds - L
    ELSE IF ds < -0.5*L then
       LET separation = ds + L
    ELSE
       LET separation = ds
    END IF
END DEF

DEF pbc(pos,L)
    LET pbc = mod(pos,L)
END DEF

SUB check_overlap
    DECLARE PUBLIC x(),y()
    DECLARE PUBLIC N,Lx,Ly
    DECLARE DEF separation
    LET tol = 1.0e-4
    FOR i = 1 to N - 1
        FOR j = i + 1 to N
            LET dx = separation(x(i) - x(j),Lx)
            LET dy = separation(y(i) - y(j),Ly)
            LET r2 = dx*dx + dy*dy
            IF r2 < 1 then
               LET r = sqr(r2)
               IF (1 - r) > tol then
                  PRINT "particles ";i;" and ";j;"overlap"
                  STOP
               END IF
            END IF
        NEXT j
    NEXT i
END SUB

SUB kinetic_energy(ke,#1)
    DECLARE PUBLIC vx(),vy()
    DECLARE PUBLIC N
    WINDOW #1
    LET ke = 0
    FOR i = 1 to N
        LET ke = ke + vx(i)*vx(i) + vy(i)*vy(i)
    NEXT i
    LET ke = 0.5*ke
    SET CURSOR 2,1
    PRINT,,,
    PRINT using "##.###": ke/N;
    PRINT,
END SUB

SUB uplist(i)
    DECLARE PUBLIC N,collision_time(),timebig
    ! look for collisions with particles j > i
    IF i = N then EXIT SUB
    LET collision_time(i) = timebig
    FOR j = i + 1 to N
        CALL check_collision(i,j)
    NEXT j
END SUB

SUB downlist(j)
    DECLARE PUBLIC collision_time()
    ! look for collisions with particles i < j
    IF j = 1 then EXIT SUB
    FOR i = 1 to j - 1
        CALL check_collision(i,j)
    NEXT i
END SUB

SUB check_collision(i,j)
    DECLARE PUBLIC x(),y(),vx(),vy()
    DECLARE PUBLIC Lx,Ly,collision_time(),partner()
    ! consider collisions between i and periodic images of j
    FOR xcell = -1 to 1
        FOR ycell = -1 to 1
            LET dx = x(i) - x(j) + xcell*Lx
            LET dy = y(i) - y(j) + ycell*Ly
            LET dvx = vx(i) - vx(j)
            LET dvy = vy(i) - vy(j)
            LET bij = dx*dvx + dy*dvy
            IF bij < 0 then
               LET r2 = dx*dx + dy*dy
               LET v2 = dvx*dvx + dvy*dvy
               LET discr = bij*bij - v2*(r2 - 1)
               IF discr > 0 then
                  LET tij = (-bij - sqr(discr))/v2
                  IF tij < collision_time(i) then
                     LET collision_time(i) = tij
                     LET partner(i) = j
                  END IF
               END IF
            END IF
        NEXT ycell
    NEXT xcell
END SUB

SUB contact(i,j,virial)
    DECLARE PUBLIC x(),y(),vx(),vy(),Lx,Ly
    DECLARE DEF separation
    ! compute collision dynamics for particles i and j at contact
    LET dx = separation(x(i) - x(j),Lx)
    LET dy = separation(y(i) - y(j),Ly)
    LET dvx = vx(i) - vx(j)
    LET dvy = vy(i) - vy(j)
    LET factor = dx*dvx + dy*dvy
    LET delvx = - factor*dx
    LET delvy = - factor*dy
    LET vx(i) = vx(i) + delvx
    LET vx(j) = vx(j) - delvx
    LET vy(i) = vy(i) + delvy
    LET vy(j) = vy(j) - delvy
    LET virial = delvx*dx + delvy*dy
END SUB

SUB show_output(t,collisions,temperature,vsum,rho,area,#1)
    WINDOW #1
    SET CURSOR 2,1
    SET COLOR "black/white"
    PRINT using "######": collisions;
    PRINT,
    PRINT using "####.##": t;
    PRINT,
    LET mean_virial = vsum/(2*t)
    LET mean_pressure = rho*temperature + mean_virial/area
    PRINT using "####.##": mean_pressure;
END SUB

SUB save_config(#2)
    DECLARE PUBLIC x(),y(),vx(),vy()
    DECLARE PUBLIC N,Lx,Ly
    WINDOW #2
    SET COLOR "black"
    ! move particles away from collision for final configuration
    CALL minimum_collision_time(i,j,tij)
    CALL move(0.5*tij)
    SET CURSOR 1,1
    INPUT prompt "name of saved configuration = ": file$
    OPEN #1: name file$, access output, create new
    PRINT #1: N
    PRINT #1: Lx
    PRINT #1: Ly
    PRINT #1: "x(i)","y(i)"
    FOR i = 1 to N
        PRINT #1, using "----.#####, ----.#####": x(i),y(i)
    NEXT i
    PRINT #1: "vx(i)","vy(i)"
    FOR i = 1 to N
        PRINT #1, using "----.#####, ----.#####": vx(i),vy(i)
    NEXT i
    CLOSE #1
END SUB

SUB set_up_windows(ncolor,#1,#2)
    DECLARE PUBLIC Lx,Ly
    OPEN #1: screen 0,1,0.90,1.0  ! output
    OPEN #2: screen 0.02,1,0.02,0.90   ! particle trajectories
    CALL compute_aspect_ratio(Lx,xwin,ywin)
    SET WINDOW 0,xwin,0,ywin
    BOX LINES 0,Lx,0,Ly
    CALL headings(#1)
    LET ncolor = 0
END SUB

SUB headings(#1)
    WINDOW #1
    SET CURSOR 1,1
    PRINT using "##########": "collisions";
    PRINT,
    PRINT using "##########": "time";
    PRINT,
    PRINT using "########.##": "<P>";
    PRINT,
    PRINT using "####.##": "T";
    PRINT,
END SUB

SUB show_positions(flag$,#2)
    DECLARE PUBLIC x(),y(),N,Lx,Ly
    IF key input then
       GET KEY k
       IF k = ord("r") then
          WINDOW #2
          CLEAR
          SET COLOR "black"
          BOX LINES 0,Lx,0,Ly
          LET flag$ = ""
       ELSEIF k = ord("s") then
          LET flag$ = "stop"
       ELSEIF k = ord("n") then
          LET flag$ = "no_show"
       END IF
    END IF
    IF flag$ <> "no_show" then
       SET COLOR "red"
       WINDOW #2
       FOR i = 1 to N
           PLOT x(i),y(i)
       NEXT i
    END IF
END SUB

SUB show_disks(ncolor,flag$,#2)
    DECLARE PUBLIC x(),y(),N,Lx,Ly
    WINDOW #2
    IF key input then
       GET KEY k
       IF k = ord("r") then
          CLEAR
          BOX LINES 0,Lx,0,Ly
          LET flag$ = ""
       ELSEIF k = ord("s") then
          LET flag$ = "stop"
       ELSEIF k = ord("n") then
          LET flag$ = "no_show"
       END IF
    END IF
    IF flag$ <> "no_show" then
       LET ncolor = mod(ncolor,6) + 1
       IF ncolor = 1 then
          CLEAR
          BOX LINES 0,Lx,0,Ly
       END IF
       SET COLOR MIX(ncolor) rnd,rnd,rnd
       SET COLOR ncolor
       FOR i = 1 to N
           BOX CIRCLE x(i)-0.5,x(i)+0.5,y(i)-0.5,y(i)+0.5
       NEXT i
    END IF
END SUB