! Some utility functions
module utils

  integer, parameter :: Lsystem = 50, Nsystem = Lsystem * Lsystem 
  integer, parameter :: Nstep = 3000, Navg = 10000
  real, parameter :: J = 1.0, kb = 1.0

  real :: h

  contains

    ! Computes the integer square root of n
    integer function isqrt(n)
      integer :: i
      i = 0
      do while (i*i <= n)
        i = i+1
      enddo
      isqrt = i-1
    end function isqrt

end module utils

! Outputs a random walk of n steps into t
subroutine random_walk(n, t)
  implicit none
  integer, intent(in) :: n
  integer, dimension(n), intent(inout) :: t

  integer :: i

  t(1) = 0
  do i = 2, n
    if (rand() < 0.5) then
      t(i) = t(i-1) - 1
    else
      t(i) = t(i-1) + 1
    endif
  enddo
end subroutine random_walk

! Computes the integral of f from a to b with n steps using the MC method
real function mc_integration(f, n, a, b)
  implicit none
  real, intent(in) :: a, b
  integer, intent(in) :: n
  real, external :: f

  integer :: i
  real :: tmp

  mc_integration = 0
  do i = 1, n
    tmp = rand() * (b-a) + a
    mc_integration = mc_integration + f(tmp)
  enddo
  mc_integration = mc_integration * (b-a) / real(n) 
end function mc_integration

! Inverse of x, to test the previous function
real function invx(x)
  implicit none
  real, intent(in) :: x

  invx = 1.0 / x
end function invx

! Translate i to be in the segment [1;l]
! Note: only works if i >= 1-l
integer function period(i, l)
  implicit none
  integer, intent(in) :: i, l

  period = mod(i-1+l,l) + 1 ! Modulo hack
end function period

! Returns the 1D-index associated to the coordinates (ix,iy) 
! in a grid of size l*l
integer function index_of(ix, iy, l)
  implicit none
  integer, intent(in) :: ix, iy, l

  index_of = (iy - 1) * l + ix
end function index_of

! Returns the 2D-indices associated to the 1D-index i in
! a grid of size l*l
subroutine reverse_index(i, l, ix, iy)
  implicit none
  integer, intent(in) :: i, l
  integer, intent(out) :: ix, iy

  iy = (i-1)/l + 1
  ix = mod(i-1,l)+1
end subroutine reverse_index

! Returns the list l(4,n) of the coordinates of the 4
! neighbours of every point in a grid of size n
subroutine mklist(list, n)
  use utils
  implicit none
  integer, dimension(4,n), intent(out) :: list
  integer, intent(in) :: n
  integer, external :: index_of, period

  integer :: l, ix, iy

  l = isqrt(n)
  do ix = 1, l
    do iy = 1, l
      list(1, index_of(ix,iy,l)) = index_of(ix,period(iy-1,l),l)
      list(2, index_of(ix,iy,l)) = index_of(period(ix+1,l),iy,l)
      list(3, index_of(ix,iy,l)) = index_of(ix,period(iy+1,l),l)
      list(4, index_of(ix,iy,l)) = index_of(period(ix-1,l),iy,l)
    enddo
  enddo
end subroutine mklist

! Computes the potential of a grid of size n given the spins s and
! the list returned by mklist
real function potential(s, n, list)
  use utils
  implicit none
  integer, dimension(4,n), intent(in) :: list
  integer, intent(in) :: n
  integer, dimension(n), intent(in) :: s

  integer :: i

  potential = 0.0
  do i = 1, n
    potential = potential - J * s(i) * s(list(2,i)) - J * s(i) * s(list(3,i)) - h * s(i)
  enddo
end function potential

! Improved Monte Carlo algorithm 
subroutine monte_carlo(s,n,list,e,t)
  use utils
  implicit none
  integer, dimension(n), intent(inout) :: s
  integer, intent(in) :: n
  integer, dimension(4,n), intent(in) :: list
  real, intent(in) :: t
  real, intent(inout) :: e

  real :: de
  integer :: i, ix, iy, iter, step

  do step = 1, n 
    i = int(rand()*real(n))+1
    de = 2.0 * h * s(i)

    do iter = 1,4
      de = de + 2.0 * J * s(i) * s(list(iter,i)) 
    enddo

    if ((de <= 0.0) .or. rand() <= exp(-de / (kb * t))) then
      s(i) = -s(i)
      e = e + de
    endif
  enddo
end subroutine monte_carlo

! Initializes the spins (at 1)
subroutine init_spins(s,n)
  use utils
  implicit none
  integer, dimension(n), intent(inout) :: s
  integer, intent(in) :: n

  integer :: i

  do i = 1,n
    s(i) = 1
  enddo
end subroutine init_spins

! Computes the magnetization of the grid of spins s
real function magnetism(s, n)
  use utils
  implicit none
  integer, intent(in) :: n
  integer, dimension(n), intent(in) :: s

  integer :: i

  magnetism = 0.0
  do i = 1, n
    magnetism = magnetism + real(s(i))
  enddo
  magnetism = magnetism / real(n)
end function magnetism

! Main program
program main
  use utils
  implicit none

  integer, dimension(4,Nsystem) :: list
  integer, dimension(Nsystem) :: spins
  real, external :: potential, magnetism 
  integer, external :: index_of

  real :: e, T, avge, avge2, avgm, avgm2
  integer :: i, i2, Ti, hi

  call srand(-34445)

  call mklist(list, Nsystem)

  h = 0.0
  T = 2.0

  open(1, file='results/energy_histogram_T1.res')
    
  call init_spins(spins, Nsystem)
    
  e = potential(spins, Nsystem, list)

  do i = 1, Nstep
    call monte_carlo(spins, Nsystem, list, e, T * J / kb)
  enddo

  do i = 1, Nstep
    call monte_carlo(spins, Nsystem, list, e, T * J / kb)
    write(1,*) e
  enddo

  close(1)
end program main
