!============================================================================
! Name        : IsingModelF.f90
! Author      : Georg Bergner
! Version     :
! Copyright   :
! Description : Simple Fortran program to run Ising model simulations.
!               Warning: several features (like random number generator)
!               not optimal for precise investigations.
!============================================================================

! Functions to find the neighbouring points in up and down direction
! using periodic boundary conditions.
integer function shiftup(num,NS)
 implicit none
 integer num
 integer NS
 shiftup=num+1
 if(num==NS) then
  shiftup=1
 endif
end function shiftup

integer function shiftdown(num,NS)
 implicit none
 integer num
 integer NS
 shiftdown=num-1
 if(num==1) then
  shiftdown=NS
 endif
end function shiftdown


program IsingModelF
    implicit none
    !  The grid size and the parameter for the scan are fixed at compile time.
    integer, parameter :: XS=40
    integer, parameter :: YS=40
    ! Minimum and maximum value for the scan.
    double precision, parameter :: Jmin=0.2
    double precision, parameter :: Jmax=0.8
    integer, parameter :: nsteps=20
    ! The thermalization has to be checked / adapted.
    integer, parameter :: nthermal=200
    integer, parameter :: numconfigs=1000
    ! Fulloutput provides the information about the thermalization.
    integer, parameter :: fulloutput=0
    integer configuration(1:XS,1:YS)
    integer x,y,n,step
    double precision J,tmpsum, average, average2,averageabs,r
    integer shiftup, shiftdown
    double precision correlator(1:XS)
    double precision correlator2(1:XS)
    double precision correlatoroutput(2*nsteps+1,1:XS)
    print *,"#Parameter Nconfigs=",numconfigs," nthermal=",nthermal," XS=",XS," YS=",YS,"\n"
    ! Scan of parameter range of J.
    do step=1,nsteps
        J=Jmin+(step-1)*(Jmax-Jmin)/dble(nsteps-1)
        ! Initial configuration.
        do x=1,XS
            do y=1,YS
                configuration(x,y)=1
            end do
        end do

        average=0
        average2=0
        averageabs=0
        ! start the updater.
        do n=1,numconfigs
            ! measure magnetization
            tmpsum=0.0
            do x=1,XS
                do y=1,YS
                    tmpsum=tmpsum+configuration(x,y)
                end do
            end do
            tmpsum=tmpsum/dble(XS*YS)

            if(fulloutput.EQ.1) then
                print*, n,tmpsum
            end if

            if(n>=nthermal) then
                average=average+tmpsum
                average2=average2+tmpsum*tmpsum
                averageabs=averageabs+dabs(tmpsum)
            end if

            ! metropolis update
            ! This is simple typewriter ordering
            ! which is not the best for all applications.
            do x=1,XS
                do y=1,YS
                    !change of local action
                    tmpsum=configuration(shiftup(x,XS),y)+configuration(shiftdown(x,XS),y)+&
                        configuration(x,shiftup(y,YS))+configuration(x,shiftdown(y,YS))
                    tmpsum=2.0*J*tmpsum*configuration(x,y)
                    ! Simple random number generator.
                    ! Not very precise results can be expected.
                    call RANDOM_NUMBER(r)
                    if (r < dexp(-tmpsum)) then
                        configuration(x,y) = -1.0*configuration(x,y)
                    end if
                end do
            end do

        end do

        ! Print out of average value for each J
        ! Limited number of observables, more can be added.
        ! Note that the error estimate is drastically underestimated
        ! since the autocorrelation time is not included.
        ! The absolute value of the order parameter provides a
        ! clearer identification of the phase transition.
        average=average/dble(numconfigs-nthermal+1)
        average2=average2/dble(numconfigs-nthermal+1)
        averageabs=averageabs/dble(numconfigs-nthermal+1)
        print *,J,average,dsqrt((average2-average*average)/dble(numconfigs-nthermal)),averageabs,&
            (average2-average*average)/dble(XS*YS)
    end do
end program IsingModelF
