function spin = metropolis_mc_ising(spin,T)

% Metropolis Monte Carlo update (single sweep) of 2D Ising ferromagnet:
% It takes matrix of spins and temperature as an input and returns 
% matrix of updated spins via Metropolis algorithm; such conversion of 
% input spins into out spins counts as one Monte Carlo time step

% INPUT:
% spin       Input matrix of Ising classical spins on the square lattice
% T          temperature in units in which J=1, k_B=1

% OUTPUT:
% spin       Output matrix of Ising classical spins on the square lattice


beta=1/T; % beta=J/(k_B*T)
size = length(spin); % find the dimension of the input matrix (which is equal to the size of the lattice)

for i=1:size %Loop over sites:
    for j=1:size
        
        % Find coordinates of neighboring spins for spin at position i,j 
        % while using periodic boundary conditions:
        ipos = mod(i,size) + 1;
        jpos = mod(j,size) + 1;
        ineg = mod(i+size-2,size) + 1;
        jneg = mod(j+size-2,size) + 1;
        
        old_spin = spin(i,j); % spin at position i,j of the input matrix
        new_spin = -spin(i,j); % flip the spin and test below if this is allowed
        
        % sum of four spins surrounding spin at position i,j
	    spin_sum = spin(i,jpos) + spin(ipos,j) + spin(i,jneg) + spin(ineg,j);
        
        % energy difference between old and new configuration of five spins 
        % consisting of spin at i,j plus its four nearest neighbors:
        old_energy=-old_spin*spin_sum;
        new_energy=-new_spin*spin_sum;
        energy_diff = beta * (new_energy - old_energy); 
        
        % update spin at position i,j by checking if new_spin is allowed 
        % in accord with the laws of statistical physics
        
        p = rand; %generate a uniform random number between 0 and 1
        
        if( (energy_diff <= 0) || (exp(-energy_diff) > p) ) % Metropolis criterion
            spin(i,j) = new_spin;  % accept the change so that spin at position i,j is now flipped
        end

    end
end %Loop over sites
end %function metropolis_mc_ising