PROGRAM hd
*     dynamics of system of hard disks
*     program based in part on Fortran program of Allen and Tildesley
      IMPLICIT NONE
      REAL*8 area, ke, rho, temperature, tij, virial, vsum, Lx, Ly,
     +       collision_time, t, timebig, vx, vy, x, y
      INTEGER IC1_, IC2_, IR_, collisions, i, j, N, partner, ishow
      CHARACTER*80 flag
      LOGICAL Ltmp1_
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      DATA IC1_/1/
      DATA IC2_/2/
      ishow = 1
      CALL GWopen(IR_, 0)
      CALL initial(vsum,rho,area)
      CALL set_up_windows(IC1_,IC2_)
      CALL kinetic_energy(ke,IC1_)
      CALL GWanchor(IR_, 1)
      CALL GWselvp(IR_, IC2_);
      IF(ishow .NE. 0) THEN
         CALL show_positions(flag)
      ELSE
         CALL show_disks(flag)
      ENDIF
      temperature = ke/N
      flag = 'show_output'
      collisions = 0                        ! number of collisions
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL minimum_collision_time(i,j,tij)
*        move particles forward and reduce collision times by tij
         CALL move(tij)
         t = t + tij
         collisions = collisions + 1
         CALL getflag(flag)
         IF(ishow .NE. 0) THEN
            CALL show_positions(flag)
         ELSE
            CALL show_disks(flag)
         ENDIF
         CALL contact(i,j,virial)           ! compute collision dynamics
         vsum = vsum + virial
*        reset collision list for relevant particles
         CALL reset_list(i,j)
         IF(flag .NE. '' .AND. flag .NE. 'no_show') THEN
            CALL GWselvp(IR_, IC1_)
            CALL show_output(t,collisions,temperature,vsum,rho,area)
            CALL GWselvp(IR_, IC2_)
            IF(flag .EQ. 'show_output' .OR. flag .EQ. 'change')
     +         flag = ''
         ENDIF
         Ltmp1_ = flag .NE. 'stop'
      END DO
      CALL save_config
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(vsum,rho,area)
      IMPLICIT NONE
      CHARACTER*80 file, heading, start
      REAL*8 ax, ay, col,
     +       nx, ny, row, vmax, vsum, rho, area, Lx, Ly, collision_time,
     +       t, timebig, vx, vy, x, y
      REAL dmy, rnd
      INTEGER IO1_, IR_, i, partner, N, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      DATA IO1_/11/
      t = 0
      WRITE(*,'(A,$)') 'read file (f) or lattice start (l) = '
      READ(*,'(A)') start
      IF(start .EQ. 'f' .OR. start .EQ. 'F') THEN
         WRITE(*,'(A,$)') 'file name = '
         READ(*,'(A)') file
         OPEN(IO1_, file = file)
         READ(IO1_,*) N
         READ(IO1_,*) Lx
         READ(IO1_,*) Ly
         READ(IO1_,'(A)') heading
         DO i = 1, N
            READ(IO1_,*) x(i),y(i)
         END DO
         READ(IO1_,'(A)') heading
         DO i = 1, N
            READ(IO1_,*) vx(i),vy(i)
         END DO
         CLOSE(IO1_)
      ELSE IF(start .EQ. 'l' .OR. start .EQ. 'L') THEN
         dmy = rnd(-1)
         WRITE(*,'(A,$)') 'N = '
         READ(*,*) N                        ! choose N so that sqr(N) an integer
         WRITE(*,'(A,$)') 'Lx = '
         READ(*,*) Lx
         WRITE(*,'(A,$)') 'Ly = '
         READ(*,*) Ly
         WRITE(*,'(A,$)') 'vmax = '
         READ(*,*) vmax
         nx = sqrt(DBLE(N))
         ny = nx
         IF(nx .GE. Lx .OR. nx .GE. Ly) THEN
            WRITE(*,*) 'box too small'
            CALL GWquit(IR_)
            STOP
         END IF
         ax = Lx/nx                         ! "lattice" spacing
         ay = Ly/ny
         i = 0
         DO col = 1, nx
            DO row = 1, ny
               i = i + 1
               x(i) = (col - 0.5D0)*ax
               y(i) = (row - 0.5D0)*ay
*              choose random positions and velocities
               vx(i) = (2*rnd(0) - 1)*vmax
               vy(i) = (2*rnd(0) - 1)*vmax
            END DO
         END DO
      END IF
      CALL check_overlap                    ! check if two disks overlap
      CALL check_momentum
      area = Lx*Ly
      rho = N/area
      timebig = 1.0D10
      vsum = 0                              ! virial sum
      DO i = 1, N
         partner(i) = N
      END DO
      collision_time(N) = timebig
*     set up initial collision lists
      DO i = 1, N
         CALL uplist(i)
      END DO
      CALL GWsetogn(IR_, 0, N + 1)
      END
C
      SUBROUTINE check_momentum
      IMPLICIT NONE
      REAL*8 vxsum, vysum, Lx, Ly, collision_time, t, timebig, vx,
     +       vy, x, y
      INTEGER i, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      vxsum = 0
      vysum = 0
      DO i = 1, N
         vxsum = vxsum + vx(i)              ! total center of mass velocity (momentum)
         vysum = vysum + vy(i)
      END DO
      vxsum = vxsum/N
      vysum = vysum/N
      DO i = 1, N
         vx(i) = vx(i) - vxsum
         vy(i) = vy(i) - vysum
      END DO
      END
C
      SUBROUTINE minimum_collision_time(i,j,tij)
*     locate minimum collision time
      IMPLICIT NONE
      REAL*8 tij, Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      INTEGER k, i, j, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      tij = timebig
      DO k = 1, N
         IF(collision_time(k) .LT. tij) THEN
            tij = collision_time(k)
            i = k
         END IF
      END DO
      j = partner(i)
      END
C
      SUBROUTINE move(tij)
      IMPLICIT NONE
      REAL*8 tij, Lx, Ly, collision_time, t, timebig, vx, vy, x, y, pbc
      INTEGER k, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      DO k = 1, N
         collision_time(k) = collision_time(k) - tij
         x(k) = x(k) + vx(k)*tij
         y(k) = y(k) + vy(k)*tij
         x(k) = pbc(x(k),Lx)
         y(k) = pbc(y(k),Ly)
      END DO
      END
C
      SUBROUTINE reset_list(i,j)
*     reset collision list for relevant particles
      IMPLICIT NONE
      REAL*8 test, Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      INTEGER k, i, j, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      DO k = 1, N
         test = partner(k)
         IF(k .EQ. i .OR. test .EQ. i .OR. k .EQ. j .OR. test .EQ. j)
     +    THEN
            CALL uplist(k)
         END IF
      END DO
      CALL downlist(i)
      CALL downlist(j)
      END
C
      FUNCTION separation(ds,L)
      IMPLICIT NONE
      REAL*8 separation, ds, L
      IF(ds .GT. 0.5D0*L) THEN
         separation = ds - L
      ELSE IF(ds .LT. -0.5D0*L) THEN
         separation = ds + L
      ELSE
         separation = ds
      END IF
      END
C
      FUNCTION pbc(pos,L)
      IMPLICIT NONE
      REAL*8 pbc, pos, L
      pbc = pos
      DO WHILE(pbc .LT. 0.0D0)
         pbc = pbc + L
      ENDDO
      pbc = mod(pbc,L)
      END
C
      SUBROUTINE check_overlap
      IMPLICIT NONE
      REAL*8 dx, dy, r, r2, tol, Lx, Ly, collision_time, t,
     +       timebig, vx, vy, x, y, separation
      INTEGER IR_, i, j, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      tol = 1.0D-4
      DO i = 1, N - 1
         DO j = i + 1, N
            dx = separation(x(i) - x(j),Lx)
            dy = separation(y(i) - y(j),Ly)
            r2 = dx*dx + dy*dy
            IF(r2 .LT. 1) THEN
               r = sqrt(r2)
               IF((1 - r) .GT. tol) THEN
                  WRITE(*,*) 'particles ', i, ' and ', j, 'overlap'
                  CALL GWquit(IR_)
                  STOP
               END IF
            END IF
         END DO
      END DO
      END
C
      SUBROUTINE kinetic_energy(ke,IC1_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 ke, Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      INTEGER IR_, i, IC1_, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      CALL GWselvp(IR_, IC1_)
      ke = 0
      DO i = 1, N
         ke = ke + vx(i)*vx(i) + vy(i)*vy(i)
      END DO
      ke = 0.5D0*ke
      WRITE(Stmp1_, '(F6.3)') ke/N
      CALL GWputtxt(IR_, 36.0, 2.0, Stmp1_)
      END
C
      SUBROUTINE uplist(i)
*     look for collisions with particles j > i
      IMPLICIT NONE
      REAL*8 Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      INTEGER j, i, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      IF(i .EQ. N) RETURN
      collision_time(i) = timebig
      DO j = i + 1, N
         CALL check_collision(i,j)
      END DO
      END
C
      SUBROUTINE downlist(j)
*     look for collisions with particles i < j
      IMPLICIT NONE
      INTEGER i, j
      IF(j .EQ. 1) RETURN
      DO i = 1, j - 1
         CALL check_collision(i,j)
      END DO
      END
C
      SUBROUTINE check_collision(i,j)
*     consider collisions between i and periodic images of j
      IMPLICIT NONE
      REAL*8 bij, discr, dvx, dvy, dx, dy, r2, tij, v2, xcell, ycell,
     +       Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      INTEGER i, j, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      DO xcell = -1, 1
         DO ycell = -1, 1
            dx = x(i) - x(j) + xcell*Lx
            dy = y(i) - y(j) + ycell*Ly
            dvx = vx(i) - vx(j)
            dvy = vy(i) - vy(j)
            bij = dx*dvx + dy*dvy
            IF(bij .LT. 0) THEN
               r2 = dx*dx + dy*dy
               v2 = dvx*dvx + dvy*dvy
               discr = bij*bij - v2*(r2 - 1)
               IF(discr .GT. 0) THEN
                  tij = (-bij - sqrt(discr))/v2
                  IF(tij .LT. collision_time(i)) THEN
                     collision_time(i) = tij
                     partner(i) = j
                  END IF
               END IF
            END IF
         END DO
      END DO
      END
C
      SUBROUTINE contact(i,j,virial)
*     compute collision dynamics for particles i and j at contact
      IMPLICIT NONE
      REAL*8 delvx, delvy, dvx, dvy, dx, dy, factor, virial, Lx, Ly,
     +       collision_time, t, timebig, vx, vy, x, y, separation
      INTEGER i, j, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      dx = separation(x(i) - x(j),Lx)
      dy = separation(y(i) - y(j),Ly)
      dvx = vx(i) - vx(j)
      dvy = vy(i) - vy(j)
      factor = dx*dvx + dy*dvy
      delvx = - factor*dx
      delvy = - factor*dy
      vx(i) = vx(i) + delvx
      vx(j) = vx(j) - delvx
      vy(i) = vy(i) + delvy
      vy(j) = vy(j) - delvy
      virial = delvx*dx + delvy*dy
      END
C
      SUBROUTINE show_output(t,collisions,temperature,vsum,rho,area)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 mean_pressure, mean_virial, t, temperature, vsum, rho,
     +       area
      INTEGER IR_, collisions, N, partner, ishow
      COMMON /GLBLI/N, partner(100), ishow
      CALL GWsetogn(IR_, 0, -(N + 2))
      mean_virial = vsum/(2*t)
      mean_pressure = rho*temperature + mean_virial/area
      WRITE(Stmp1_, '(I6,$)') collisions
      CALL GWputtxt(IR_, 1.0, 2.0, Stmp1_)
      WRITE(Stmp1_, '(F7.2,$)') t
      CALL GWputtxt(IR_, 12.0, 2.0, Stmp1_)
      WRITE(Stmp1_, '(F7.2,$)') mean_pressure
      CALL GWputtxt(IR_, 24.0, 2.0, Stmp1_)
      CALL GWflush(IR_, -(N + 2))
      CALL GWsetogn(IR_, 0, N + 1)
      END
C
      SUBROUTINE save_config
      IMPLICIT NONE
      CHARACTER*80 file,Stmp1_
      REAL*8 tij, Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      INTEGER IO1_, IR_, i, j, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      DATA IO1_/11/
*     move particles away from collision for final configuration
      CALL minimum_collision_time(i,j,tij)
      CALL move(0.5*tij)
      CALL GWinput(IR_, 'name of saved configuration = ', Stmp1_)
      IF(IR_ .EQ. 0) RETURN
      READ(Stmp1_,'(A)') file
      OPEN(IO1_, file = file)
      WRITE(IO1_,*) N
      WRITE(IO1_,*) Lx
      WRITE(IO1_,*) Ly
      WRITE(IO1_,*) '       x(i)      y(i)'
      DO i = 1, N
         WRITE(IO1_, '(2XF10.5, F10.5)') x(i), y(i)
      END DO
      WRITE(IO1_,*) '      vx(i)     vy(i)'
      DO i = 1, N
         WRITE(IO1_, '(2XF10.5, F10.5)') vx(i), vy(i)
      END DO
      CLOSE(IO1_)
      END
C
      SUBROUTINE set_up_windows(IC1_,IC2_)
      IMPLICIT NONE
      REAL*8 Lx, Ly, collision_time, t,
     +       timebig, vx, vy, x, y
      REAL xwin, ywin, GWaspect
      INTEGER IR_, IC1_, IC2_, N, partner, ishow
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      CALL GWvport(IR_, 0.0, 0.90, GWaspect(1), 1.0)
                                            ! output
      CALL GWindow(IR_, 0.0, 3.0, 80.0, 0.0)
      CALL GWsavevp(IR_, IC1_)
      CALL headings
      CALL GWvport(IR_, 0.02*GWaspect(1), 0.02, GWaspect(1), 0.90)
                                            ! particle trajectories
      CALL compute_aspect_ratio(REAL(Lx),xwin,ywin)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWrect(IR_, 0.0, 0.0, REAL(Lx), REAL(Ly))
      CALL GWsetmrk(IR_, 6, 1.0, 13, -1, 4)
      CALL GWsavevp(IR_, IC2_)
      END
C
      SUBROUTINE headings
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      INTEGER IR_
      WRITE(Stmp1_, '(A)') 'collisions'
      CALL GWputtxt(IR_, 1.0, 1.0, Stmp1_)
      WRITE(Stmp1_, '(A)') 'time'
      CALL GWputtxt(IR_, 12.0, 1.0, Stmp1_)
      WRITE(Stmp1_, '(A)') '<P>'
      CALL GWputtxt(IR_, 24.0, 1.0, Stmp1_)
      WRITE(Stmp1_, '(A)') '   T'
      CALL GWputtxt(IR_, 36.0, 1.0, Stmp1_)
      END
C
      SUBROUTINE show_positions(flag)
      IMPLICIT NONE
      REAL*8 Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      INTEGER IR_, i, N, partner, ishow
      CHARACTER*80 flag
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      IF(flag .NE. 'no_show') THEN
         CALL GWerase(IR_, N + 1, 0)
         DO i = 1, N
            CALL GWsetpxl(IR_, REAL(x(i)), REAL(y(i)), 13)
         END DO
      END IF
      END
C
      SUBROUTINE show_disks(flag)
      IMPLICIT NONE
      CHARACTER*80 flag
      REAL*8 Lx, Ly, collision_time, t, timebig, vx, vy, x, y
      REAL rnd
      INTEGER IR_, i, ncolor, N, partner, ishow, krgb
      SAVE ncolor
      COMMON /GLBLD/Lx, Ly, collision_time(100), t, timebig,
     +             vx(100), vy(100), x(100), y(100)
      COMMON /GLBLI/N, partner(100), ishow
      DATA ncolor/0/
      IF(flag .NE. 'no_show') THEN
         ncolor = mod(ncolor,6) + 1
         IF(ncolor .EQ. 1) CALL GWrefresh(IR_)
         ncolor = 
     +          krgb(int(rnd(0)*256), int(rnd(0)*256), int(rnd(0)*256))
         CALL GWsetpen(IR_, ncolor, -1, -1, 4)
         DO i = 1, N
            CALL GWsetogn(IR_, 0, -i)
            CALL GWputmrk(IR_, REAL(x(i)), REAL(y(i)))
            CALL GWflush(IR_, -i)
         END DO
      END IF
      END
C
      SUBROUTINE getflag(flag)
      IMPLICIT NONE
      CHARACTER*80 flag
      INTEGER IR_, k, N, partner, ishow
      LOGICAL TBkeyinput
      COMMON /GLBLI/N, partner(100), ishow
      IF(TBkeyinput()) THEN
         CALL TBgetkey(k)
         flag = 'show_output'
         IF(k .EQ. ICHAR('r') .OR. k .EQ. ICHAR('R')) THEN
            CALL GWerase(IR_, 0, -1)
            CALL GWrefresh(IR_)
         ELSE IF(k .EQ. ICHAR('s') .OR. k .EQ. ICHAR('S')) THEN
            flag = 'stop'
         ELSE IF(k .EQ. ICHAR('n') .OR. k .EQ. ICHAR('N')) THEN
            flag = 'no_show'
         ELSE IF(k .EQ. ICHAR('c') .OR. k .EQ. ICHAR('C') .OR. 
     +           k .EQ. ICHAR('\t')) THEN
            CALL GWerase(IR_, 0, -1)
            IF(ishow .NE. 0) THEN
                ishow = 0
            ELSE
                ishow = 1
            ENDIF
         ENDIF
      ENDIF
      END