PROGRAM boltzmann
! Metropolis algorithm for a particle in one dimension
DIM P(-100 to 100),accum(3)
CALL initial(v,E,beta,mcs,nequil,delta,nbin,delv)
FOR imcs = 1 to nequil            ! equilibrate system
    CALL metropolis(v,E,beta,delta,accept)
NEXT imcs
CALL initialize_sums(P(),accum(),accept,nbin)
FOR imcs = 1 to mcs
    CALL metropolis(v,E,beta,delta,accept)
    ! accumulate data after each trial change
    CALL data(P(),accum(),v,E,nbin,delv)
NEXT imcs
CALL output(P(),accum(),mcs,accept,nbin,delv)
END

SUB initial(v0,E0,beta,mcs,nequil,delta,nbin,delv)
    RANDOMIZE
    INPUT prompt "temperature = ": T
    LET beta = 1/T
    INPUT prompt "number of Monte Carlo steps = ": mcs
    LET nequil = int(0.1*mcs)
    INPUT prompt "initial speed = ": v0
    LET E0 = 0.5*v0*v0            ! initial kinetic energy
    INPUT prompt "maximum change in velocity = ": delta
    LET vmax = 10*sqr(T)
    LET nbin = 20                 ! number of bins
    LET delv = vmax/nbin          ! velocity interval
END SUB

SUB initialize_sums(P(),accum(),accept,nbin)
    FOR ibin = -nbin to nbin
        LET P(ibin) = 0
    NEXT ibin
    FOR i = 1 to 3
        LET accum(i) = 0
    NEXT i
    LET accept = 0
END SUB

SUB metropolis(v,E,beta,delta,accept)
    LET dv = (2*rnd - 1)*delta    ! trial velocity change
    LET vtrial = v + dv           ! trial velocity
    LET dE = 0.5*(vtrial*vtrial - v*v)      ! trial energy change 
    IF dE > 0 then
       IF exp(-beta*dE) < rnd then
          EXIT SUB                ! step not accepted
       END IF
    END IF
    LET v = vtrial
    LET accept = accept + 1
    LET E = E + dE
END SUB

SUB data(P(),accum(),v,E,nbin,delv)
    LET accum(1) = accum(1) + E
    LET accum(2) = accum(2) + E*E
    LET accum(3) = accum(3) + v
    LET ibin = round(v/delv)
    LET P(ibin) = P(ibin) + 1
END SUB

SUB output(P(),accum(),mcs,accept,nbin,delv)
    LET accept = accept/mcs
    PRINT "acceptance probability ="; accept
    LET vave = accum(3)/mcs
    PRINT "mean velocity ="; vave
    LET Eave = accum(1)/mcs
    PRINT "mean energy ="; Eave
    LET E2ave = accum(2)/mcs
    LET sigma2 = E2ave - Eave*Eave
    PRINT "sigma_E = "; sqr(sigma2)
    PRINT
    PRINT " v ", "P(v)"
    PRINT
    LET v = -nbin*delv
    FOR ibin = -nbin to nbin
        IF p(ibin) > 0 then
           LET prob = p(ibin)/mcs
           PRINT v,
           PRINT using "--.###": prob
        END IF
        LET v = v + delv
    NEXT ibin
END SUB