function e = ising_energy ( spin )
%
%  There is a linear algebra way of doing this:
%  e = S'*T*S
%
[ m, n ] = size ( spin );
  
e = 0;
  
for i = 1 : m-1
  for j = 1 : n
    e = e - spin(i,j) * spin(i+1,j);
  end
end

i = m;
for j = 1 : n
  e = e - spin(i,j) * spin(1,j);
end

for j = 1 : n-1
  for i = 1 : m
    e = e - spin(i,j) * spin(i,j+1);
  end
end

j = n;
for i = 1 : m
  e = e - spin(i,j) * spin(i,1);
end
