! Ising model in two dimensions on a square lattice
! calculate the magnetization and energy as functions of temperature
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program ising_2d
   option nolet
   library "sgfunc*","sglib*"
   dim mag(100),temp(100),energy(100),spin(0,0)
     ! mag() and energy() contain running averages for magnetization and energy
     ! at temperature temp()   spin(i,j) = value of spin at site i,j
   randomize
   call initialize(spin,nmax,t_init,t_final,n_points,n_passes,current_m,current_e)
     ! set up for a temperature sweep
   dt = (t_final - t_init) / (n_points - 1)
   t = t_init
   call init_spin(spin,nmax,current_m,current_e,1)
   for i = 1 to n_points    ! sweep temperature
      call calculate(spin,nmax,t,n_passes,current_m,current_e,mag(i),energy(i))
      temp(i) = t   ! NOT "time"
      print t,mag(i),energy(i)
      t = t + dt
   next i
   mat redim mag(n_points),temp(n_points),energy(n_points)
   call display(temp,mag)     ! display final results
   get key z
   call display(temp,energy)
end
! initialize variables
! nmax = size of lattice  t_init,t_final = starting and ending temperature
! n_points = number of temperatures to consider
! n_passes = number of passes through the lattice at each temperature
sub initialize(spin(,),nmax,t_init,t_final,n_points,n_passes,current_m,current_e)
   nmax = 10   ! lattice size is nmax by nmax
   mat redim spin(nmax,nmax)
   t_init = 1
   t_final = 5
   n_points = 9      ! change this to change the number of temperatures considered
   n_passes = 100    ! make this larger to get more accurate averages
   call init_spin(spin,nmax,current_m,current_e,1)
end sub
! initialize spin directions
! also initialize magnetization and energy
sub init_spin(spin(,),nmax,m,e,flag)
   for i = 1 to nmax
      for j = 1 to nmax
         if flag = 1 then     ! if flag = 1 start with an ordered array of spins
            spin(i,j) = 1
         else                 ! else start with a random array
            if(rnd >= 0.5) then
               spin(i,j) = 1
            else 
               spin(i,j) = -1
            end if
         end if
      next j
   next i
   if flag = 1 then           ! initialize mag and energy
      m = nmax^2
      e = -2 * nmax^2
   else 
      m = 0
      e = 0
      for i = 1 to nmax
         for j = 1 to nmax
            m = m + spin(i,j)
            call neighbors(spin,nmax,i,j,sum)
            e = - spin(i,j) * sum / 2
         next j
      next i
   end if
end sub
! do the real work here at temperature t
sub calculate(spin(,),nmax,t,n_passes,current_mag,current_energy,ave_mag,ave_energy)
   dim m(0),e(0),ts(0)
   mat redim m(n_passes+1),e(n_passes+1),ts(n_passes+1)
   m_tmp = 0           ! set up for running averages of mag and energy
   e_tmp = 0
   n = 0
   m(1) = current_mag
   e(1) = current_energy
   ts(1) = 0
   for i = 1 to n_passes   ! make this many sweeps through the lattice
      for j = 1 to nmax    ! sweep along each row and column
         for k = 1 to nmax
            call neighbors(spin,nmax,j,k,sum) !check each spin for possible flip
            call test_for_flip(spin,j,k,t,sum,current_mag,current_energy)
         next k
      next j
      m(i+1) = current_mag
      e(i+1) = current_energy
      ts(i+1) = i
      if i > n_passes / 2 then
         m_tmp = m_tmp + current_mag
         e_tmp = e_tmp + current_energy
         n = n + 1
      end if
   next i
   ave_mag = m_tmp / n
   ave_energy = e_tmp / n
end sub
! find the net alignment of the nearest neighbors of spin j,k
! note periodic boundary conditions
sub neighbors(spin(,),nmax,j,k,sum)
   x1 = j
   y1 = k + 1
   call test_for_edge(y1,nmax)   ! have to deal with edges carefully
   x2 = j
   y2 = k - 1
   call test_for_edge(y2,nmax)
   x3 = j + 1
   y3 = k
   call test_for_edge(x3,nmax)
   x4 = j - 1
   y4 = k 
   call test_for_edge(x4,nmax)
   sum = spin(x1,y1) + spin(x2,y2) + spin(x3,y3) + spin(x4,y4)
end sub
sub test_for_edge(x,nmax)
   if x < 1    then x = nmax
   if x > nmax then x = 1
end sub
! apply the Monte-Carlo flipping rules here
sub test_for_flip(spin(,),j,k,t,sum,current_mag,current_energy)
   denergy = 2 * spin(j,k) * sum
   if denergy < 0 then        ! always flip if it lowers the energy
      spin(j,k) = - spin(j,k)
      current_mag = current_mag + 2 * spin(j,k)
      current_energy = current_energy + denergy
   else                       ! if energy would increase, flip with probability
      if exp(-denergy / t) > rnd then   ! equal to Boltzmann factor
         spin(j,k) = - spin(j,k)
         current_mag = current_mag + 2 * spin(j,k)
         current_energy = current_energy + denergy
      end if
   end if
end sub
! plot final results
sub display(x(),y())
   set background color "white"
   set color "black"
   clear
   call datagraph(x,y,4,1,"black")
end sub