PROGRAM R2
! compute mean square displacement of one particle
! by averaging over all possible origins
DIM x(1000),y(1000),R2cum(20),norm(20)
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

SUB initial(Lx,Ly,R2cum(),ndiff)
    LET Lx = 5
    LET Ly = 5
    LET ndiff = 10                ! maximum time difference
    FOR idiff = 1 to ndiff
        LET R2cum(idiff) = 0
    NEXT idiff
END SUB

SUB read_data(x(),y(),ndata)
    ! read file for position of particle at regular intervals 
    OPEN #1: name "xy.dat",access input
    LET t = 0
    DO while more #1
       LET t = t + 1
       INPUT #1: x(t),y(t)
    LOOP
    CLOSE #1
    LET ndata = t                 ! # of data points
END SUB

SUB displacements(x(),y(),Lx,Ly,R2cum(),norm(),ndata,ndiff)
    DECLARE DEF separation        ! function same as in Program md
    FOR idiff = 1 to ndiff
        FOR i = 1 to ndata - idiff
            LET dx = separation(x(i+idiff) - x(i),Lx)
            LET dy = separation(y(i+idiff) - y(i),Lx)
            LET R2cum(idiff) = R2cum(idiff) + dx*dx + dy*dy
            LET norm(idiff) = norm(idiff) + 1
        NEXT i
    NEXT idiff
END SUB

SUB normalization(ndiff,R2cum(),norm())
    PRINT "time difference","   <R2>"
    PRINT
    FOR idiff = 1 to ndiff
        IF R2cum(idiff) > 0 then
           LET R2bar = R2cum(idiff)/norm(idiff)
           PRINT idiff,R2bar
        END IF
    NEXT idiff
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