real L,dE,TM,TE,E,M,M1,M2,X,E1,C1,C2,Esq,TEsq,C
real AvgE,AvgM,TMsq,Msq
Integer S(0:N+1,0:N+1)
R=rand(23123)
do while (.true.)
read(*,*,end=10)T
do i=0,N+1
do j=0,N+1
S(i,j)=1
enddo
enddo
AJ=1.0
BJ=0.1
TM=0.0
TE=0.0
TEsq=0.0
TMsq=0.0
do kmc=1,1100
M=0.0
E=0.0
Esq=0.0
Msq=0.0
do i=1, N
do j=1, N
If (i==1) S(N+1,j) = S(i,j)
If (i==N) S(0,j) = S(i,j)
If (j==1) S(i,N+1) = S(i,j)
If (j==N) S(i,0) = S(i,j)
Bq=0.0
do i1=1,N
do j1=1,N
ix = iabs(i1 - i)
iy = iabs(j1 - j)
If (N-ix < ix) ix = N-ix
If (N-iy < iy) iy = N-iy
B = sqrt(1.*ix**2 + 1.*iy**2)
If(i1.eq.i.and.j1.eq.j) Then
Bq = Bq
Else
Bq=Bq+S(i1,j1)/B
Endif
enddo
enddo
dE=2*AJ*S(i,j)*(S(i+1,j)+S(i-1,j)+S(i,j+1)+S(i,j-1))-2.*BJ*S(i,j)*Bq
if (dE<=0) then
S(i,j)=(-1)*S(i,j)
else
L=exp(-dE/T)
R=rand()
if (R < L) then
S(i,j)=(-1)*S(i,j)
else
endif
endif
enddo
enddo
if (kmc>1000) then
do i=1, N
do j=1, N
M=M+S(i,j)
Bq=0.0
do i1=1,N
do j1=1,N
ix = iabs(i1 - i)
iy = iabs(j1 - j)
If (N-ix < ix) ix = N-ix
If (N-iy < iy) iy = N-iy
B=sqrt(1.*ix**2+1.*iy**2)
If(i1.eq.i.and.j1.eq.j) Then
Bq = Bq
Else
Bq=Bq+S(i1,j1)/B
Endif
enddo
enddo
E1=-AJ*S(i,j)*(S(i+1,j)+S(i-1,j)+S(i,j+1)+S(i,j-1))+BJ*S(i,j)*Bq
E=E+E1
enddo
enddo
E = E/2.0
Esq= E*E
Msq= M*M
TM=TM+M
TE=TE+E
TESq=TEsq+Esq
TMsq=TMsq+Msq
else
endif
enddo
Vinv = 1./(1.*N*N)
C2=TEsq*Vinv*Vinv/100.
AvgE=TE*Vinv/100.
C1=AvgE*AvgE
C=(1./T*T)*((1.*C2)-(1.*C1))
AvgM=Vinv*TM/100.
M2=AvgM*AvgM
M1=TMsq*Vinv*Vinv/100.
X=(1./T)*(M1-M2)
write(11,*)T,AvgE
write(12,*)T,AvgM
write(13,*)T,C
write(14,*)T,X
write(*,*)T,AvgM,AvgE,X
open
(2,file='lattice.dat',status='unknown',access='append')
Do i1=1,N
write(2,*) ((1+S(i1,j1))/2,j1=1,N)
Enddo
write(2,*)'------------------------------------------'
close(2)
enddo
10 continue
stop
end