//Metropolis pour le modèle d'Ising


N=50;
beta=0.1;
//beta=0.5*log(1+sqrt(2))+0.2;
//beta=7;
T=5000;
pas=100;



clf();
B=1:N; C=ones(1,N);
S=[];
for i=1:N
  I=[B;i*C;C;C];
  S=[S I];
end,
A=grand(1,N*N,'uin',0,1);
A(1:N)=0;
A(N*N-N+1:N*N)=1;
for i=1:N
    A((i-1)*N+1)=0;
    A((i-1)*N+N)=1;
end,
B=A*2+3*ones(1,N*N); //couleurs 3 5
xrects(S,B);
plot(1,0);
plot(N+1,N);
  
function k=coor(s,t);
  k=(s-1)*N+t;
endfunction

  
couleur=[];
for i=1:T
    if (modulo(i,100)==0) then
      disp(i);
   end,
   a=grand(1,2,'uin',2,N-1);
  x1=pmodulo(a(1,1)-2,N)+1;
  x2=pmodulo(a(1,1),N)+1;
  y1=pmodulo(a(1,2)-2,N)+1;
  y2=pmodulo(a(1,2),N)+1;
  a0=1-A(coor(a(1,1),a(1,2)));
  b0=A(coor(a(1,1),a(1,2)));
  a1=A(coor(x1,a(1,2)));
  a2=A(coor(a(1,1),y1));
  a3=A(coor(x2,a(1,2)));
  a4=A(coor(a(1,1),y2));
  r=exp(-beta*((a1-a0)^2+(a2-a0)^2+(a2-a0)^2+(a3-a0)^2+(a4-a0)^2)+beta*((a1-b0)^2+(a2-b0)^2+(a2-b0)^2+(a3-b0)^2+(a4-b0)^2));
  u=grand(1,1,'unf',0,1);
  p=min(1,r);
  if (u<p) then
    A(coor(a(1,1),a(1,2)))=1-A(coor(a(1,1),a(1,2)));
  end,
  if (pmodulo(i,pas)==0) then
    B=A*2+3*ones(1,N*N); //couleurs 3 5
    xrects(S,B);
  end,
  //end, end,
end,

  
  
