PROGRAM demon
! demon algorithm for the d = 1 Ising model in zero magnetic field 
DIM spin(1000)
LIBRARY "mygraphics"
CALL initial(N,spin(),E,Ed,M,mcs,Ecum,Edcum,Mcum,M2cum,accept)
CALL setupscreen(N,spin(),up$,down$)
FOR imcs = 1 to mcs
    CALL change(N,spin(),E,Ed,M,accept,up$,down$)
    CALL data(E,Ed,M,Ecum,Edcum,Mcum,M2cum)
NEXT imcs
CALL averages(N,Ecum,Edcum,Mcum,M2cum,mcs,accept)
END

SUB initial(N,spin(),E,Ed,M,mcs,Ecum,Edcum,Mcum,M2cum,accept)
    RANDOMIZE
    INPUT prompt "number of spins = ": N
    ! choose total energy to be multiple of 4J
    ! coupling constant J is unity
    INPUT prompt "desired total energy = ": Etot
    LET Etot = 4*int(Etot/4)
    INPUT prompt "number of Monte Carlo steps per spin = ": mcs
    ! initial configuration of spins in minimum energy state
    FOR isite = 1 to N
        LET spin(isite) = 1
    NEXT isite
    LET M = N                     ! net magnetization
    ! compute initial system energy
    LET E = -N                    ! periodic boundary conditions
    LET Ed = (Etot - E)
    PRINT "total energy = "; E + Ed
    ! initialize sums
    LET Ecum = 0
    LET Edcum = 0
    LET Mcum = 0
    LET M2cum = 0
END SUB

SUB setupscreen(N,spin(),up$,down$)
    LET r = N/(2*pi)
    CALL compute_aspect_ratio(r+2,xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
    LET dtheta = 2*pi/N
    SET COLOR "red"
    BOX AREA 1,1+0.5,1,1+0.5
    BOX KEEP 1,1+0.5,1,1+0.5 in up$
    CLEAR
    SET COLOR "blue"
    BOX AREA 1,1+0.5,1,1+0.5
    BOX KEEP 1,1+0.5,1,1+0.5 in down$
    CLEAR
    FOR i = 1 to N
        CALL showspin(N,spin(i),i,up$,down$)
    NEXT i
END SUB

SUB change(N,spin(),E,Ed,M,accept,up$,down$)
    FOR i = 1 to N
        LET isite = int(rnd*N + 1)     ! random spin
        ! determine neighboring spin values
        IF isite = 1 then
           LET left = spin(N)
        ELSE
           LET left = spin(isite - 1)
        END IF
        IF isite = N then
           LET right = spin(1)
        ELSE
           LET right = spin(isite + 1)
        END IF
        ! trial energy change 
        LET de =  2*spin(isite)*(left + right)
        IF de <= Ed then
           ! spin flip dynamics
           LET spin(isite) = -spin(isite)
           LET M = M + 2*spin(isite)
           LET accept = accept + 1     ! number of changes accepted
           LET Ed = Ed - de
           LET E = E + de
        END IF
        CALL showspin(N,spin(isite),isite,up$,down$)
    NEXT i
END SUB

SUB data(E,Ed,M,Ecum,Edcum,Mcum,M2cum)
    ! accumulate data
    LET Ecum = Ecum + E
    LET Edcum = Edcum + Ed
    LET Mcum = Mcum + M
    LET M2cum = M2cum + M*M
END SUB

SUB averages(N,Ecum,Edcum,Mcum,M2cum,mcs,accept)
    SET COLOR "black/white"
    LET norm = 1/mcs              ! collected data after every attempt
    LET Edave = Edcum*norm
    PRINT "mean demon energy ="; Edave
    LET T = 4/log(1 + 4/Edave)
    PRINT "T ="; T
    LET Eave = Ecum*norm
    PRINT "<E> = "; Eave
    LET Mave = Mcum*norm
    PRINT "<M> ="; Mave
    LET M2ave = M2cum*norm
    PRINT "<M*M> ="; M2ave
    LET accept_prob = accept*norm/N
    PRINT "acceptance probability ="; accept_prob
END SUB

SUB showspin(N,dir,isite,up$,down$)
    LET r = N/(2*pi)
    LET theta = isite/r
    LET x = r*cos(theta)
    LET y = r*sin(theta)
    IF dir = 1 then
       BOX SHOW up$ at x,y
    ELSE
       BOX SHOW down$ at x,y
    END IF
END SUB