%-----------------------------------------------------------------------
% File: ising1.m                                        (P335: 3/05) 
%
% Monte Carlo Simulation for the two-dimensional Ising model on a square
% lattice with zero magnetic field and periodic boundary condition.
%
% Spin configurations are displayed through the run.
% Magnetization vs MC step are plotted at end of run.
%
% Please set the temperature (T), system size (L), and number of
% MC cycles (NCYCLE)to your liking.
% Larger systems and more MC cycles give better results
%           ... but they also cost a lot more CPU time.
%
% Try temperatures in the range 1.0 -> 10.0
%
% Try running T=2.0 several times in a row.
%    Do you get different final M values?
%    How many cycles are required for equlibration?
%
% PRESS F5-key to run this program
%
%--------------------------------------------------------------------------
clear;  %clear variable names

% PLEASE ADJUST THESE PARAMETERS:
T=3.00;     %temperature
L=25;       % number of rows in the lattice
NCYCLE=500; % number of Monte Carlo steps


%---------------------------------------------------------%
% Main Program
%---------------------------------------------------------%
%Define some other parameters and start RNG:
A=L*L; %total number of lattice sites and spins
j=1;   % interaction constant: j=+1 for ferromagnetic 
	   % systems, j=-1 for antiferromagnetic systems
h=0;   % applied magnetic field (zero for this program)
rand('state',sum(100*clock)); % set state for random number gen.

% set up the initial configuration
spin=zeros(L,L);  % initialize variables to zero
for y=1:L
   for x=1:L
      if rand>0.5   % spins are randomly set to 1 or 0
         spin(x,y)=1; 
      else
         spin(x,y)=-1;
      end 
   end
end
figure(1) % show the initial configurations
imagesc(spin,[0.2,1]); colormap(gray);
title(['T=',num2str(T) ,'   (Starting Configuration)'])
axis image
pause(1)

% identify nearest neighbors (periodic boundary conditions)
for i=1:L
   if i==1
      ip(i)=L;
   else
      ip(i)=i-1;
   end
   if i==L
      is(i)=1;
   else
      is(i)=i+1;
   end
end
Mst=zeros(1,NCYCLE);  % initialize the magnetization to zero

for k=1:NCYCLE  % Monte Carlo steps   
   for x=1:L  % visit the sites in turn
      for y=1:L
         %calculate the interaction energy (n.n. interaction)
         spinsum=spin(ip(x),y)+spin(is(x),y)...
            +spin(x,ip(y))+spin(x,is(y)); 
         Espin = -j*spin(x,y)*spinsum - h*spin(x,y);
         %if Espin>0 then energy is gained by flipping a spin
         %therefore flip the spin. 
         %if Espin<0 then energy is lost by a spin flip, 
         %therefore flip only with Boltzmann probability
         if Espin>0 
            spin(x,y)=-spin(x,y);
         else
            boltz=exp(2*Espin/T);
            if rand<boltz 
               spin(x,y)=-spin(x,y);
            end
         end
      end
   end
   Mst(k)=sum(sum(spin))/(L*L); %determine magnetization per site
   % show the configuration after every 10 MC steps
   if mod(k,10)==0
   	imagesc(spin,[0.2,1]); colormap(gray);
	title(['T=',num2str(T) ,'   (',num2str(k), '  MC sweeps)'])
   	axis image
   	pause(0.05)
   end
end  % end of Monte Carlo steps

% show the magnetization as a function of MC steps
figure(2)
plot(Mst)
axis([0,NCYCLE,-1,1])
xlabel('MC step')
ylabel('Magnetization per Site')
title(['2D Ising Model:  T=',num2str(T),',  L=',num2str(L),' square lattice'])