#C14.txt
Help:=proc(): print(` Neis(c,M,N) , NewA(A,c), OSM(A,M,N,x,y), Nick(M,N,x,y,K1)`): end:

with(Statistics):

###stuff from C13.txt
Help13:=proc(): print(`RM(M,N), PlotMat(A), Ene(A) , Mag(A)`):
print(`Prob(A,x,y) , Vecs(N) , Mats(M,N) ,  ParF(M,N,x,y)`):
end:

RM:=proc(M,N) local i,j,ra:
ra:=rand(0..1):
subs(0=-1,[seq([seq(ra(),i=1..M)],j=1..N)]):
end:

with(plots):

#PlotMat(A): plots the matrix A where 1 is black and -1 is white
PlotMat:=proc(A) local i,j,d:

if A[1][1]=1 then
d:=pointplot([[1,1]],axes=none,color=black,symbol=circle):
else
d:=pointplot([[1,1]],axes=none,color=red,symbol=circle):
fi:

for j from 2 to nops(A[1]) do
if A[1][j]=1 then
d:=d,pointplot([[1,j]],color=black,symbol=circle):
else
d:=d,pointplot([[1,j]],color=red,symbol=circle):
fi:
od:

for i from 2 to nops(A) do
for j from 1 to nops(A[i]) do
if A[i][j]=1 then
d:=d,pointplot([[i,j]],color=black,symbol=circle):
else
d:=d,pointplot([[i,j]],color=red,symbol=circle):
fi:
od:
od:



display(d);

end:


#Ene(A): The energy of the matrix A of {1,-1}
Ene:=proc(A) local i,j:
add(add( A[i][j]*A[i][j+1]   , j=1..nops(A[i])-1) ,i=1..nops(A))+
add(add( A[i][j]*A[i+1][j]   , j=1..nops(A[1]) )  ,i=1..nops(A)-1)+
add(A[1][i]*A[-1][i],i=1..nops(A[1]))+
add(A[i][1]*A[i][-1],i=1..nops(A)):
end:

#Mag(A): the magnetization of the spin-config. A
Mag:=proc(A): add(add(A)): end:

#Vecs(N): the set of all 2^N {1,-1} lists of length N
Vecs:=proc(N) local S,v:
option remember:
if N=0 then
 RETURN({[]}):
fi:
S:=Vecs(N-1):

{seq([op(v),1], v in S),seq([op(v),-1], v in S)}:

end:




#Mats(M,N): the set of all M by N {1,-1} matrices
Mats:=proc(M,N) local S,S1,A,v:
if M=0 then
 RETURN({[]}):
fi:

S:=Mats(M-1,N):
S1:=Vecs(N):

{seq(seq([op(A),v],v in S1),A in S)}:

end:

#x=exp(k/T) (T=abs. temp. (in Kelvin)) y=exp(External Magnetication)
#x^Ene*y^Mag
#ProbA(A,x,y):Unnormalied prob.
ProbA:=proc(A,x,y): x^Ene(A)*y^Mag(A):end:

#ParF(M,N,x,y): the Gibbs partition function of M by N Ising model
ParF:=proc(M,N,x,y) local S,A:
S:=Mats(M,N):
add(ProbA(A,x,y),A in S):
end:

###end from C13.txt

#NeisDrZ(c,M,N): Inputs a location c=[i,j] and pos. integers M and N outputs the SET
#of four neighbors of the cell [i,j], {[i \pm 1,j ], [i, j \pm 1] but with toral convention
NeisDrZ:=proc(c,M,N) local i,j,S,S1,s:

i:=c[1]: j:=c[2]:

S:={[i-1 ,j],[i+1,j],[i,j-1],[i,j+1]}:

S1:={}:


for s in S do
 if s[1]=0 then
    S1:=S1 union  {[N,s[2]]}:
  elif s[1]=N+1 then
    S1:=S1 union  {[1,s[2]]}:
   else
    S1:=S1 union {s}:
  fi:


if s[2]=0 then
    S1:=S1 union  {[N,s[2]]}:
  elif s[1]=N+1 then
    S1:=S1 union  {[1,s[2]]}:
   else
    S1:=S1 union {s}:
  fi:
od:
end:

#Neis(c,M,N): Inputs a location c=[i,j] and pos. integers M and N outputs the SET
#of four neighbors of the cell [i,j], {[i \pm 1,j ], [i, j \pm 1] but with toral convention
#CODE by Alison Bu and Yukun Yao, with improvements by Charles Kenney
Neis:=proc(c,M,N) local i,j,S,k,a:

i:=c[1]: j:=c[2]:

S:={[(i-2 mod M)+1,j ],[(i mod M)+1, j ],[i, (j-2 mod N)+1 ] ,[i,(j mod N)+1]}:

#for k from 1 to nops(S)  do
# if S[k][1]=0 then S[k][1]:=M fi:
 #if S[k][2]=0 then S[k][2]:=N fi:
#od:
#S:
end:

#NewA(A,c): The result of flipping the spin at location c in the matrix A
#followed by the change in Energy, and change in Magnetization
NewA:=proc(A,c) local i,j,s, A1,Ne,DiffE, DiffM,a:
i:=c[1]: j:=c[2]:
s:=A[i][j]:

A1:=[op(1..i-1,A),
 [op(1..j-1,A[i]), -s, op(j+1..nops(A[i]), A[i])],  
 op(i+1..nops(A),A)]:

Ne:=Neis(c,nops(A),nops(A[1])):

DiffE:=-2*s*add(A[a[1]][a[2]], a in Ne):
DiffM:=-2*s:
[A1,DiffE,DiffM]:
end:

#OSM(A,M,N,x,y): one step in the Metropolis algorithm inputs an M by N matrix of {1,-1}
#and pos. numbers x,y outputs one step in the Metropolis algorithm according to
#the importance x^Ene(A)*y^Mag(A)
OSM:=proc(A,M,N,x,y) local c,Zidong, A1,Rat,r,U:
c:=[rand(1..M)(), rand(1..N)()]:

Zidong:=NewA(A,c):
A1:=Zidong[1]:

Rat:=x^Zidong[2]*y^Zidong[3]:

r:=Sample(RandomVariable(Uniform(0,1)),1)[1]:

if r<Rat then
  A1:
else
 A:
fi:

end:

#Nick(M,N,x,y,K1): Starts with a random M by N {1,-1} matrix with physical parameters x and y
#and  a positive integer K1 runs Metropolis K1 times
Nick:=proc(M,N,x,y,K1) local A,i:
A:=RM(M,N):

for i from 1 to K1 do
 A:=OSM(A,M,N,x,y):
od:
A:
end: