function Ising=AL_Ising(dim,b,timesteps)
%% simulate Ising spin model
%features:
% - square lattice 
% - periodic boundary conditions 
% - odd number of elements
% - Moore neighborhood 
%Version 1.0
%Example:
% to simulate the Ising model on a 7x7 lattice with inverse Temperature
% b=1/5 with 1000 time steps enter:
%                           I=AL_Ising(7,1/5,1e3);

% The magnetization as function of time (time steps correspond to array index) can be
% accessed via:
%                           I.Marray

% To visualize the time series type:
%                           plot(I.Marray)

Ising.dim               = dim;  %size of 2d-lattice (dimxdim)
Ising.b                 = b;    %inverse temperature
Ising.Jij               = 1;    %coupling strength
%% initialize system S with uniform distribuion of states {-1,1}
Ising.S   =   rand( [Ising.dim,Ising.dim] );
Ising.S   =   Ising.S - 0.5;
Ising.S( Ising.S < 0 )   =    -1;
Ising.S( Ising.S > 0 )   =    1;

%% calculate lookup table (LUT) of nearest neighbors (Moore neighborhood)
Ising.sizS    = size(Ising.S);
[x,y]=meshgrid(1:Ising.dim,1:Ising.dim);
dim 	= Ising.dim;

Ising.neighbors  =  [   sub2ind( Ising.sizS, x, mod( ( (y - 1) - 1),dim ) + 1),... 
                        sub2ind( Ising.sizS, x, mod( ( (y - 1) + 1),dim ) + 1 ),...  
                        sub2ind( Ising.sizS, mod( ( (x - 1) + 1), dim ) + 1, y ),... 
                        sub2ind( Ising.sizS, mod( ( (x - 1) - 1), dim ) + 1, y ),...     
                        sub2ind( Ising.sizS, mod( ( (x - 1) + 1), dim ) + 1, mod( ( (y - 1) - 1),dim ) + 1 ),...
                        sub2ind( Ising.sizS, mod( ( (x - 1) - 1), dim ) + 1, mod( ( (y - 1) - 1),dim ) + 1),...
                        sub2ind( Ising.sizS, mod( ( (x - 1) + 1), dim ) + 1, mod( ( (y - 1) + 1),dim ) + 1),...
                        sub2ind( Ising.sizS, mod( ( (x - 1) - 1), dim ) + 1, mod( ( (y - 1) + 1),dim ) + 1)    ];
                    
Ising.neighborhood = reshape(Ising.neighbors,dim,dim,8);    
%% some type conversion of LUT for better efficiency in time loop 
Ising.neighbors={};
for i=1:dim
    for j=1:dim
        Ising.neighbors{i,j}     =   squeeze(Ising.neighborhood(i,j,:)); 
    end
end
%% initialize array for storage of magnetization data
Ising.Marray = []; 
%% start simulation
for tstep = 1:timesteps
    %% shuffle indices to pick random matrix elements and still touch each
    %% element only once per time step
    ix      =   randperm( numel( Ising.S ) ); 
    %% loop randomly through matrix S
    for i=ix
            Ising.currentIx     = i;
            Ising.currentS      = Ising.S( Ising.currentIx ); 
            
        % calculate local field at current position ix(i)
            Ising.h     =  sum( Ising.S( Ising.neighbors{ Ising.currentIx } ) ) * Ising.Jij;
        % calc probability for state shift
            Ising.p     = 1/(1+exp( -2*Ising.b*Ising.h ));
        % perform shift of spin
        if rand(1) < Ising.p
            Ising.S( Ising.currentIx )  = 1;
        else
            Ising.S( Ising.currentIx )  = -1;
        end
    end
    %% print status 
    if mod(tstep, round(timesteps/10)) == 0
        disp([num2str(round(tstep/timesteps * 100)),'% done'])
    end
    %% store resulting total magentization (normalized to be between 0 and 1) at the end of array Marray
    Ising.Marray(end + 1)   =   sum(sum(Ising.S))/numel(Ising.S);
end
%%  save results to hard disk 
save(['dim',num2str(dim),'b',num2str(b),'.mat'],'Ising')
end