function im = gibbs(im,beta)

nRows = size(im, 1);
nCols = size(im, 2);

for i = 1:nRows
    for j = 1:nCols
        
        % voisins V4
        i1=mod(i-2,nRows)+1;
        i2=mod(i,nRows)+1;
        j1=mod(j-2,nCols)+1;
        j2=mod(j,nCols)+1;
       
	% evaluate energy 
        p0 = exp( energy( 0, im(i1, j), im(i2, j), im(i, j1), im(i, j2) ) * beta );
        p1 = exp( energy( 1, im(i1, j), im(i2, j), im(i, j1), im(i, j2) ) * beta );
        
	% conditional transition
        if( rand() < p0/(p0+p1) )
            im(i,j) = -1;
        else
            im(i,j) = 1;
        end

    end
end

end

function val = energy( me, up, down, left, right )
	% val = -1 * me * left - me * right; 
	val = -1 * me * left - me * right + me * up + me * down; 
end
