************************************************************************
C Program: Ising.f
C Purpose: Using Monte Carlo integration to calculate the energy and
C magnetization for the 2D Ising model.  I'm taking J=1.
************************************************************************

      include 'ran2.f'
      implicit none

      integer NITER,NTHERM,NLAT
      parameter(NITER=1000,NTHERM=100)
      parameter(NLAT=30)

      real*8 spins(NLAT,NLAT),temp,energy,magnet
      real*8 sum_energy,sum_magnet,ran2
      real*8 engArray(NITER),magArray(NITER)
      integer iseed,i,j
      
      open(10,file="input.dat",status="old")
      open(80,file="data.dat",status="unknown")

      read(10,*) temp,iseed
*-----------------------------------------------------------------------
C Randomly set spins to start.
*-----------------------------------------------------------------------
      do i=1,NLAT
         do j=1,NLAT
            if(ran2(iseed).lt.0.5) then
               spins(i,j)=1.
            else
               spins(i,j)=-1.
            end if
         end do
      end do
      call obsv(spins,NLAT,energy,magnet)

*-----------------------------------------------------------------------
C Thermalize the configuration.
*-----------------------------------------------------------------------
      do i=1,NTHERM
         call monte(spins,temp,iseed,NLAT)
      end do

*-----------------------------------------------------------------------
C Run the calculation.
*-----------------------------------------------------------------------
      do i=1,NITER
         call monte(spins,temp,iseed,NLAT)
         call obsv(spins,NLAT,energy,magnet)
         sum_energy=sum_energy+energy
         sum_magnet=sum_magnet+magnet
         engArray(i)=energy
         magArray(i)=magnet

      end do
      do i=1,NITER
         write(80,*) engArray(i),magArray(i)
      end do
      print*,sum_energy/float(NITER),sum_magnet/float(NITER)
      
      end

************************************************************************
C Subroutine: monte
C Purpose: Sweeps through the lattice, updating spins with Monte Carlo.
************************************************************************
      subroutine monte(spins,temp,iseed,NLAT)
      implicit none
      real*8 spins(NLAT,NLAT)
      real*8 flipped,energy_old,energy_new,temp,sum_nn,delta,ratio,ran2

      integer i,j,im1,jm1,ip1,jp1,iseed,NLAT
     
      do i=1,NLAT
         im1=i-1
         ip1=i+1
         if(im1.le.0) im1=im1+NLAT
         if(ip1.ge.NLAT+1) ip1=ip1-NLAT
         do j=1,NLAT
            jm1=j-1
            jp1=j+1
            if(jm1.le.0) jm1=jm1+NLAT
            if(jp1.ge.NLAT+1) jp1=jp1-NLAT
            flipped=-spins(i,j)
            sum_nn=spins(ip1,j)+spins(im1,j)+spins(i,jp1)+spins(i,jm1)
            energy_new=flipped*sum_nn
            energy_old=spins(i,j)*sum_nn
            delta=(energy_new-energy_old)/temp
            if(delta.ge.0.0) then
               spins(i,j)=flipped
            else
               ratio=exp(delta)
               if(ratio.ge.ran2(iseed)) spins(i,j)=flipped
            end if
         end do
      end do

      return
      end

************************************************************************
C Subroutine obsv
C Purpose: Calculate the energy and the magnetization
************************************************************************
      subroutine obsv(spins,NLAT,energy,magnet)
      implicit none
      
      real*8 spins(NLAT,NLAT)
      real*8 energy,magnet
      integer nlat,i,ip1,j,jp1
      
      magnet=0.0
      energy=0.0
      do i=1,NLAT
         ip1=i+1
         if(ip1.ge.NLAT+1) ip1=ip1-NLAT
         do j=1,NLAT
            jp1=j+1
            if(jp1.ge.NLAT+1) jp1=jp1-NLAT
            magnet=magnet+spins(i,j)
            energy=energy+spins(i,j)*(spins(ip1,j)+spins(i,jp1))
         end do
      end do
      magnet=magnet/float(NLAT*NLAT)
      energy=energy/float(NLAT*NLAT)
      return
      end
      
      
