program ising
! program to simulate the 2Dising model generating random
! this version moves Hamiltonian calc to user def fn
! it uses four 2-D arrrays:
! s() is the inital spin values
! ju() is J "up" -  the vertical interaction values
! jl() is J "left" -  the horizontal interaction values
! sn() is the new spin values
! uses /dev/urandom so needs some modernish *nix to run

integer, parameter :: n = 8
character (len=2), parameter  :: anx=' X', ano=' O', anp=' +',anm=' -',dsp='  '
integer :: i,j,k,l,delt,up,dn,lf,rt, un, rsd
integer, allocatable :: seed(:)
integer, dimension(n,n) :: s, ju, jl, sn
real, dimension(n,n) :: ss, si, id
id = 1
!s  = reshape([-1,-1,-1, 1, 1,-1, 1, 1, 1],[n,n])
!ju = reshape([ 1, 1, 1,-1, 1, 1, 1,-1,-1],[n,n])
!jl = reshape([-1, 1, 1, 1, 1,-1,-1,-1,-1],[n,n])

call random_seed(size=rsd)
allocate(seed(rsd))
open(newunit=un, file="/dev/urandom", access="stream",form="unformatted", action="read", status="old",iostat=istat)
read(un) seed
close(un)
call random_seed(put=seed)

!set random s
call random_number(ss)
si = sign(id,ss - 0.5)
s = int(si)

!set random s
call random_number(ss)
si = sign(id,ss - 0.5)
ju = int(si)

!set random s
call random_number(ss)
si = sign(id,ss - 0.5)
jl = int(si)


sn = s
!print *, 's='
!do i = 1,n
!   print *, (s(i,j), j=1,n)
!end do
!print *, 'ju='
!do i = 1,n
!   print *, (ju(i,j), j=1,n)
!end do
!print *, 'jl='
!do i = 1,n
!   print *, (jl(i,j), j=1,n)
!end do


do 20 k=1,n
do 10 l=1,n

if (k-1 < 1) then
   up = s(n,l)*ju(k,l)
else
   up = s(k-1,l)*ju(k,l)
end if

if (k+1 > n) then
   dn = s(1,l)*ju(1,l)
else
   dn = s(k+1,l)*ju(k+1,l)
end if

if (l-1 < 1) then
   lf = s(k,n)*jl(k,l)
else
   lf = s(k,l-1)*jl(k,l)
end if

if (l+1 > n) then
   rt = s(k,1)*jl(k,1)
else
   rt = s(k,l+1)*jl(k,l+1)
end if

delt = up + dn + lf +rt

if (delt > 0) then
   sn(k,l) = 1
else if (delt < 0) then
   sn(k,l) = -1
else
   sn(k,l) = s(k,l)
end if

!print *, 'up= ', up
!print *, 'dn= ', dn
!print *, 'lf= ', lf
!print *, 'rt= ', rt
!print *, 'delt= ', delt
!print *, 'sn',k,l,'=', sn(k,l)


10 continue
20 continue




!print *, 'sn='
!do i = 1,n
!   print *, (sn(i,j), j=1,n)
!end do

do i = 1,n
   print *, (pnm(jl(i,j))//xno(s(i,j)), j=1,n)
   print *, (dsp//pnm(ju(i,j)), j=1,n)
end do

print *, 'function first H= ', hamiltonian(s,ju,jl)

do i = 1,n
   print *, (pnm(jl(i,j))//xno(sn(i,j)), j=1,n)
   print *, (dsp//pnm(ju(i,j)), j=1,n)
end do

print *, 'function new H= ', hamiltonian(sn,ju,jl)

contains

!calculate H ----------------------------------

integer function hamiltonian(hs,hu,hl)
   implicit none
   integer, dimension(n,n), intent(in) :: hs,hu,hl
   integer :: p,q, hh


hh = 0

!the up Js
do 120 p=1,n
do 110 q=1,n

if (p==1) then
   hh = hh - hs(n,q)*hu(p,q)*hs(p,q)
else
   hh = hh - hs(p-1,q)*hu(p,q)*hs(p,q)
end if


110 continue
120 continue

!the left Js
do 140 p=1,n
do 130 q=1,n

if (q==1) then
   hh = hh - hs(p,n)*hl(p,q)*hs(p,q)
else
   hh = hh - hs(p,q-1)*hl(p,q)*hs(p,q)
end if

130 continue
140 continue

hamiltonian = hh

end function hamiltonian

!end calculate H

character (len=2) function xno(a)
   implicit none
   integer, intent(in) :: a
   character (len=2) :: xnov

if (a == 1) then
   xnov = ano
else
   xnov = anx
end if

xno = xnov

end function xno

character (len=2) function pnm(a)
   implicit none
   integer, intent(in) :: a
   character (len=2) :: pnmv

if (a == 1) then
   pnmv = anp
else
   pnmv = anm
end if

pnm = pnmv

end function pnm


end program