function [s,Esave,Msave] = ising(beta,J,B,p,N);
%
%   [s,Esave,Msave] = ising(beta,J,B,p,N);
%
%   beta = inverse temp
%   J = interaction strength
%   B = external field
%   p = initial density of positive spins
%   N = lattice size

%N = 32;
%p = 0.01;
s = (rand(N,N)<p)*2-1;

%beta = 0.41;
%J = 1;
%B = 0;

numIter = 200*N^2;
Esave = zeros(1,numIter+1);
Msave = zeros(1,numIter+1);
E = energy(s,J,B);
M = sum(sum(s));
Esave(1) = E;
Msave(1) = M;

%subplot 121;
%pcolor(s); axis square; title 'IC'; drawnow;

for iter=1:numIter,
    
    i = randi(N);
    j = randi(N);
    ip1 = mod(i,N)+1;
    im1 = mod(i-2,N)+1;
    jp1 = mod(j,N)+1;
    jm1 = mod(j-2,N)+1;
    dE = 2*J*s(i,j)*(s(im1,j)+s(ip1,j)+s(i,jm1)+s(i,jp1)) + 2*B*s(i,j);
    
    if rand(1) < exp(-beta*dE)
        s(i,j) = -s(i,j);
        E = E + dE;
        M = M + 2*s(i,j);
    end
    Esave(iter+1) = E;
    Msave(iter+1) = M;
end
Esave = Esave/N^2;
Msave = Msave/N^2;

%pcolor(s); axis square;
%title (['Ising model, \beta=' num2str(beta) ', J=' num2str(J) ', B=' num2str(B)]);

%subplot(122);
%plot([Esave;Msave]');
%axis square;
%legend('Energy/N^2','Spin/N^2','location','best');
%title (['Ising model, \beta=' num2str(beta) ', J=' num2str(J) ', B=' num2str(B)]);


