program krysztal
implicit none
integer,parameter :: sizeofbox=20
integer, dimension (sizeofbox,sizeofbox) :: grid
integer,dimension(sizeofbox):: NI,PI
integer, parameter ::seed=42
integer :: i,j,k,l=1,m=1,E, fi_new
real ::electric_field,dE,boltzman,w,R, dEold, dEnew, pie=3.141, n1=1.5, n2=1.7,neff_node
real, dimension(2)::wlist
real, dimension(4200)::neff
real, dimension(360) :: P2(-180:180)
!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 tablicy P2
do i=-180,180
    P2(i)=0.5*(3.0*cos(real(i)*pie/180.0)*cos(real(i)*pie/180)-1)
end do

!inicjacja grida
do i=1,sizeofbox
    do j=1,sizeofbox
        grid(i,j)=0
    end do
end do


do j=1,sizeofbox
    grid(1,j)=0
    grid(sizeofbox,j)=0
end do
open(unit=1,file='average5x5')

print*, grid
do E=0,10
    print*, E
    electric_field=0.5+E*0.1
    do k=1,500000
        do i=2,sizeofbox-1
            do j=1,sizeofbox
                fi_new=grid(i,j)+int((rand()-0.5)*3)
                if (fi_new>90) then
                    fi_new=fi_new-180
                end if

                if (fi_new<-90) then
                    fi_new=fi_new+180
                end if
                dEold=-20.0*(P2(grid(i,j)-grid(NI(i),j))+P2(grid(i,j)-grid(i,NI(j)))&
                +P2(grid(i,j)-grid(i,PI(j)))+P2(grid(i,j)-grid(PI(i),j)))&
                +(-electric_field*electric_field*P2(grid(i,j)-90))

                dEnew=-20.0*(P2(fi_new-grid(NI(i),j))+P2(fi_new-grid(i,NI(j)))&
                +P2(fi_new-grid(i,PI(j)))+P2(fi_new-grid(PI(i),j)))&
                +(-electric_field*electric_field*P2(fi_new-90))
                dE=dEnew-dEold
                wlist(1)=1
                boltzman=exp(-dE)
                wlist(2)=boltzman

                w=minval(wlist)
                R=rand()
                if (R<=w) then
                    grid(i,j)=fi_new
                end if
            end do
        end do



        if (k>80000 .AND. (MOD(k,100) == 0)) then
            m=1
            neff_node=0
            do i=1,sizeofbox
                do j=1,sizeofbox
                    neff_node=neff_node+(n1*n2)/sqrt(n1**2*cos(grid(i,j)*pie/180)**2&
                    +n2**2*sin(grid(i,j)*pie/180)**2)
                end do
            end do
            neff(l)=neff_node/(sizeofbox*sizeofbox)
            l=l+1
        end if
    end do
    l=1
    print*, sum(neff)/4200
    write(1,*) electric_field , sum(neff)/4200

end do

write(1,*) grid
close(1)
end