function [ north, south, east, west ] = ising_neighbors ( spin, i, j )
%
%  function [ north, south, east, west ] = ising_neighbors ( spin, i, j )
%
%  returns the values of spin at the four neighbors of atom (I,J).
%
%  SPIN is an M by N array of spins;
%  I, J are the row and column indices of the atom we are interested in.
%
%  NORTH, SOUTH, EAST and WEST are the values of SPIN at the neighbor
%    atom.  Neighbors "wrap around", so the north neighbor of atom (1,5)
%    is atom (M,5).
%
[m, n] = size ( spin );

if ( i < m ) 
  k = i + 1;
else
  k = 1;
end
north = spin(k,j);

if ( 1 < i ) then
  k = i - 1;
else
  k = m;
end
south = spin(k,j);

if ( j < n ) then
  k = j + 1;
else
  k = 1;
end
east = spin(i,k);

if ( 1 < j ) then
  k = j - 1;
else
  k = n;
end
west = spin(i,k);