      PROGRAM R2
*     compute mean square displacement of one particle
*     by averaging over all possible origins
      IMPLICIT NONE
      REAL*8 Lx, Ly, R2cum(20), x(1000), y(1000)
      INTEGER ndata, ndiff, norm(20)
      DATA norm/20*0/
      CALL initial(Lx,Ly,R2cum,ndiff)
      CALL read_data(x,y,ndata)
      CALL displacements(x,y,Lx,Ly,R2cum,norm,ndata,ndiff)
      CALL normalization(ndiff,R2cum,norm)
      END
C
      SUBROUTINE initial(Lx,Ly,R2cum,ndiff)
      IMPLICIT NONE
      REAL*8 Lx, Ly, R2cum(20)
      INTEGER idiff, ndiff
      Lx = 5
      Ly = 5
      ndiff = 10                            ! maximum time difference
      DO idiff = 1, ndiff
         R2cum(idiff) = 0
      END DO
      END
C
      SUBROUTINE read_data(x,y,ndata)
*     read file for position of particle at regular intervals 
      IMPLICIT NONE
      REAL*8 x(1000), y(1000)
      INTEGER IO1_, t, ndata
      DATA IO1_/11/
      OPEN(IO1_, file = 'xy.dat')
      t = 0
      DO WHILE(.TRUE.)
         t = t + 1
         READ(IO1_, *, END = 999) x(t),y(t)
      END DO
  999 CLOSE(IO1_)
      ndata = t - 1                         ! # of data points
      END
C
      SUBROUTINE displacements(x,y,Lx,Ly,R2cum,norm,ndata,ndiff)
      IMPLICIT NONE
      REAL*8 dx, dy, x(1000), y(1000), Lx, Ly, R2cum(20), separation
      INTEGER i, idiff, norm(20), ndata, ndiff
      DO idiff = 1, ndiff
         DO i = 1, ndata - idiff
            dx = separation(x(i+idiff) - x(i),Lx)
            dy = separation(y(i+idiff) - y(i),Ly)
            R2cum(idiff) = R2cum(idiff) + dx*dx + dy*dy
            norm(idiff) = norm(idiff) + 1
         END DO
      END DO
      END
C
      SUBROUTINE normalization(ndiff,R2cum,norm)
      IMPLICIT NONE
      REAL*8 R2bar, R2cum(20)
      INTEGER idiff, ndiff, norm(20)
      WRITE(*,*) 'time difference     <R2>'
      WRITE(*,*)
      DO idiff = 1, ndiff
         IF(R2cum(idiff) .GT. 0) THEN
            R2bar = R2cum(idiff)/norm(idiff)
            WRITE(*,'(I12,1X,F12.6)') idiff, R2bar
         END IF
      END DO
      END
C
      FUNCTION separation(ds,L)             ! function same as in Program md
      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
