PROGRAM ising
! Monte Carlo simulation of the Ising model on the square lattice
! using the Metropolis algorithm
DIM spin(32,32),w(-8 to 8),accum(10)
DIM Ce(0:20),Cm(0:20),esave(100),msave(100)
LIBRARY "csgraphics"
CALL initial(N,L,T,spin(,),mcs,nequil,w(),E,M)
CALL plotspins(L,spin(,))
FOR i = 1 to nequil               ! equilibrate system
    CALL Metropolis(N,L,spin(,),E,M,w(),accept)
NEXT i
CALL initialize(accum(),accept)
CALL initialize_correl(Ce(),Cm(),esave(),msave(),nsave)
FOR pass = 1 to mcs               ! accumulate data while updating spins 
    CALL Metropolis(N,L,spin(,),E,M,w(),accept)
    CALL data(E,M,accum())
    CALL correl(Ce(),Cm(),E,M,esave(),msave(),pass,nsave)
NEXT pass
CALL output(T,N,mcs,accum(),accept)
CALL c_output(N,Ce(),Cm(),accum(),mcs,nsave)
END

SUB initial(N,L,T,spin(,),mcs,nequil,w(),E,M)
    RANDOMIZE
    INPUT prompt "linear dimension of lattice = ": L
    LET N = L*L                   ! number of spins
    ! temperature measured in units of J/k
    INPUT prompt "temperature = ": T
    INPUT prompt "# MC steps per spin for equilibrium = ": nequil
    INPUT prompt "# MC steps per spin for data = ": mcs
    LET M = 0
    FOR y = 1 to L                ! random initial configuration
        FOR x = 1 to L
            IF rnd < 0.5 then
               LET spin(x,y) = 1  ! spin up
            ELSE
               LET spin(x,y) = -1
            END IF
            LET M = M + spin(x,y)      ! total magnetization
        NEXT x
    NEXT y
    LET E = 0
    FOR y = 1 to L                ! compute initial energy
        IF y = L then
           LET up = 1             ! periodic boundary conditions
        ELSE
           LET up = y + 1
        END IF
        FOR x = 1 to L
            IF x = L then
               LET right = 1
            ELSE
               LET right = x + 1
            END IF
            LET sum = spin(x,up) + spin(right,y)
            LET E = E - spin(x,y)*sum  ! total energy
        NEXT x
    NEXT y
    ! compute Boltzmann probability ratios
    FOR dE = -8 to 8 step 4
        LET w(dE) = exp(-dE/T)
    NEXT dE
END SUB

SUB plotspins(L,spin(,))
    CLEAR
    CALL compute_aspect_ratio(L+1,xwin,ywin)
    SET WINDOW 0,xwin,0,ywin
    FOR i = 1 to L
        FOR j = 1 to L
            IF spin(i,j) = 1 then
               SET COLOR "red"
            ELSE
               SET COLOR "green"
            END IF
            BOX AREA i,i+1,j,j+1
        NEXT j
    NEXT i
END SUB

SUB initialize_correl(Ce(),Cm(),esave(),msave(),nsave)
    LET nsave = 10
    FOR i = 1 to nsave
        LET esave(i) = 0
        LET msave(i) = 0
    NEXT i
    FOR i = 1 to nsave
        LET Ce(i) = 0
        LET Cm(i) = 0
    NEXT i
END SUB

SUB initialize(accum(),accept)
    ! use array to save accumulated values of magnetization and
    ! energy. Array used to make it easier to add other quantities      
    FOR i = 1 to 5
        LET accum(i) = 0
    NEXT i
    LET accept = 0
END SUB

SUB Metropolis(N,L,spin(,),E,M,w(),accept)
    DECLARE DEF DeltaE
    ! one Monte Carlo step per spin
    FOR ispin = 1 to N
        LET x = int(L*rnd + 1)    ! random x coordinate
        LET y = int(L*rnd + 1)    ! random y coordinate
        LET dE = DeltaE(x,y,L,spin(,))      ! compute change in energy
        IF rnd <= w(dE) then
           LET spin(x,y) = -spin(x,y)  ! flip spin
           LET accept = accept + 1
           LET M = M + 2*spin(x,y)
           LET E = E + dE
           LET dir = spin(x,y)
           CALL show_change(x,y,dir)
        END IF
    NEXT ispin
END SUB

FUNCTION DeltaE(x,y,L,spin(,))
    ! periodic boundary conditions
    IF x = 1 then
       LET left = spin(L,y)
    ELSE
       LET left = spin(x - 1,y)
    END IF
    IF x = L then
       LET right = spin(1,y)
    ELSE
       LET right = spin(x + 1,y)
    END IF
    IF y = 1 then
       LET down = spin(x,L)
    ELSE
       LET down = spin(x,y - 1)
    END IF
    IF y = L then
       LET up = spin(x,1)
    ELSE
       LET up = spin(x,y + 1)
    END IF
    LET DeltaE = 2*spin(x,y)*(left + right + up + down)
END DEF

SUB show_change(i,j,dir)
    IF dir = 1 then
       SET COLOR "red"
    ELSE
       SET COLOR "green"
    END IF
    BOX AREA i,i+1,j,j+1
END SUB

SUB data(E,M,accum())
    ! accumulate data after every Monte Carlo step per spin
    LET accum(1) = accum(1) + E
    LET accum(2) = accum(2) + E*E
    LET accum(3) = accum(3) + M
    LET accum(4) = accum(4) + M*M
    LET accum(5) = accum(5) + abs(M)
END SUB

SUB output(T,N,mcs,accum(),accept)
    LET norm = 1/(mcs*N)          ! averages per spin
    LET accept = accept*norm
    LET eave = accum(1)*norm
    LET e2ave = accum(2)*norm
    LET mave = accum(3)*norm
    LET m2ave = accum(4)*norm
    LET abs_mave = accum(5)*norm
    CLEAR
    SET BACKGROUND COLOR "black"
    SET COLOR "yellow"
    PRINT "temperature = "; T
    PRINT "acceptance probability = "; accept
    PRINT "mean energy per spin = "; eave
    PRINT "mean squared energy per spin = "; e2ave
    PRINT "mean magnetization per spin = "; mave
    PRINT "mean of absolute magnetization per spin = "; abs_mave
    PRINT "mean squared magnetization per spin = "; m2ave
END SUB

SUB save_config(N,L,T,spin(,))
    INPUT prompt "name of file for last configuration = ": file$
    OPEN #2: name file$, access output, create new
    PRINT #2: T
    FOR y = 1 to L
        FOR x = 1 to L
            PRINT #2: spin(x,y)
        NEXT x
    NEXT y
    CLOSE #2
END SUB

SUB read_config(N,L,T,spin(,))
    INPUT prompt "filename?": file$
    OPEN #1: name file$, access input
    INPUT #1: T
    FOR y = 1 to L
        FOR x = 1 to L
            INPUT #1: spin(x,y)
        NEXT x
    NEXT y
    CLOSE #1
END SUB

SUB correl(Ce(),Cm(),E,M,esave(),msave(),pass,nsave)
    ! accumulate data for time correlation functions
    ! save last nsave values of M and E
    ! index0 = array index for earliest saved time
    LET index0 = mod(pass-1,nsave) + 1
    IF pass > nsave then
       ! compute Ce and Cm after nsave values are saved
       LET index = index0
       FOR tdiff = nsave to 1 step -1
           LET Ce(tdiff) = Ce(tdiff) + E*esave(index)
           LET Cm(tdiff) = Cm(tdiff) + M*msave(index)
           LET index = index + 1
           IF index > nsave then LET index = 1
       NEXT tdiff
    END IF
    ! save latest value in array position of earliest value
    LET esave(index0) = E
    LET msave(index0) = M
END SUB

SUB c_output(N,Ce(),Cm(),accum(),mcs,nsave)
    ! compute time correlation functions
    LET ebar = accum(1)/mcs
    LET e2bar = accum(2)/mcs
    LET Ce(0) = e2bar - ebar*ebar
    LET mbar = accum(3)/mcs
    LET m2bar = accum(4)/mcs
    LET Cm(0) = m2bar - mbar*mbar
    LET norm = 1/(mcs - nsave)
    PRINT
    PRINT "t","Ce(t)","Cm(t)"
    PRINT
    FOR tdiff = 1 to nsave
        ! correlation functions defined so that C(t=0) = 1
        ! and C(infinity) = 0
        LET Ce(tdiff) = (Ce(tdiff)*norm - ebar*ebar)/Ce(0)
        LET Cm(tdiff) = (Cm(tdiff)*norm - mbar*mbar)/Cm(0)
        PRINT tdiff,Ce(tdiff),Cm(tdiff)
    NEXT tdiff
END SUB