% Simulazione del modello di Ising con l'algoritmo di Metropolis
% per eseguire dalla command window: [x,mag]=ising;
% per compilare e creare un file eseguibile dal prompt: mcc -B sgl ising.m
function [x,mag]=ising(n,nsweep)

if nargin<2
    nsweep=500;
    if nargin<1
        n=20;
    end
end

ti=cputime;
bj=0.6;
for i=1:2
    pdelta(i)=exp(-4*bj*i);
end
n2=n^2;
% inizializza il reticolo in modo casuale
s=2*randint(n)-1;
nstep=nsweep*n2+1;
count=0;
%ciclo sui passi di MC
for k=1:nstep
% seleziona casualmente gli indici di uno spin 
    ss=randint(2,1,[1,n]);
    i=ss(1);
    j=ss(2);
    % tiene conto delle condizioni periodiche al contorno
    im=mod(i+n-2,n)+1;
    ip=mod(i+n,n)+1;
    jm=mod(j+n-2,n)+1;
    jp=mod(j+n,n)+1;
    deltaE=2*s(i,j)*(s(i,jm)+s(i,jp)+s(im,j)+s(ip,j));
    if deltaE<=0
        s(i,j)=-s(i,j);
    else
     % genera un numero casuale per decidere sul rigetto
        if rand<pdelta(deltaE/4)
            s(i,j)=-s(i,j);
        end
    end
    if mod(k,n2)==1
        count=count+1;
        mag(count)=sum(sum(s));
    end
end
mag=mag./n^2;
x=0:nsweep;
plot(x,mag);
dt=cputime-ti