program ising
implicit none
integer,parameter :: sizeofbox=10
real, dimension (sizeofbox,sizeofbox) :: grid
integer,dimension(sizeofbox):: NI,PI
integer, parameter ::seed=42
integer :: i,j,k,l=1,m=1,t
real ::temperature,dE,boltzman,w,R
real, dimension(2)::wlist
real, dimension(4200)::magnetisation,mcs, temp_energy
real, dimension(sizeofbox**2):: hamilton
real::average_magnetisation, average_magnetisation_2, average_magnetisation_4, energy, energy_2
CHARACTER(LEN=15) :: str
!INICJACJA TABLICY NAJBLIŻSZYCH SĄSIADÓW
do i=1,sizeofbox
    NI(i)=i+1
    PI(i)=i-1
end do
NI(sizeofbox)=1
PI(1)=sizeofbox
!DONE


!inicjacja grida
call random_number(grid)
do i=1,sizeofbox
    do j=1,sizeofbox
        IF (grid(i,j)>0.5) THEN
        grid(i,j)=1
        ELSE
        grid(i,j)=-1
    end IF
end do
end do
!DONE

open(unit=1,file='average10x10')
write(1,*)'Średniamagnetyzacja', 'Średniamagnetyzacja^2', 'sredniamagnetyzacja4', 'e', 'e2'

!metropolis
do t=1,100
    temperature=1.5+t*0.02
    write(str, *) t
    do k=1,500000
        do i=1,sizeofbox
            do j=1,sizeofbox
                dE=-2*grid(i,j)*(grid(NI(i),j)+grid(PI(i),j)+grid(i,NI(j))+grid(i,PI(j)))
                wlist(1)=1
                boltzman=exp(dE/temperature)
                wlist(2)=boltzman
                w=minval(wlist)
                R=rand()
                if (R<=w) then
                    grid(i,j)=grid(i,j)*(-1)
                end if
            end do
        end do
        if (k>80000 .AND. (MOD(k,100) == 0)) then
            magnetisation(l)=sum(grid)/sizeofbox**2
            mcs(l)=k
            !print*, sum(grid)/sizeofbox**2
            m=1
            do i=1,sizeofbox
                do j=1,sizeofbox
                    hamilton(m)=-2*grid(i,j)*(grid(NI(i),j)+grid(PI(i),j)+grid(i,NI(j))+grid(i,PI(j)))
                    m=m+1
                end do
            end do
            temp_energy(l)=sum(hamilton)
            l=l+1
        end if
    end do
l=1

average_magnetisation=sum(abs(magnetisation))/4200.0
average_magnetisation_2=sum(abs(magnetisation**2))/4200.0
average_magnetisation_4=sum(abs(magnetisation**4))/4200.0
energy=sum(temp_energy)/4200.0
energy_2=sum(temp_energy**2)/4200.0
!print*, 'SREDNIA MAGNETYZACJA', average_magnetisation
print*, average_magnetisation
write(1,*) temperature, average_magnetisation, average_magnetisation_2,average_magnetisation_4,energy, energy_2
!do i=1,4200
!    write(1,*) mcs(i), magnetisation(i)
!end do
end do

close(1)
!print*, SIZE(average_magnetisation)
end program ising