function z = area_poly(x, y)
%This returns the area of a closed polygon.  The area is calculated by
%applying Green's Theorem in the special case N = x
%=> int(x*dy) = int(int(dx*dy)).
%For a polygon this gives area = sum((x2 + x1)*(y2 - y1))/2 over all of
%the vertices of the polygon.  Note that since we are de
%   INPUT:
%x: vector of the x positions of the polygon.
%y: vector of the y positions of the polygon.
%   OUTPUT:
%z: the area of the polygon.
%
%Note that since we are dealing with a closed polygon the first point
%specified should also be the last one specified, i.e., x(1) == x(len)
%and y(1) == y(len).  If this is not true of the input vectors then they
%are extended suitably.

%     Copyright J. V. Mansbridge, CSIRO, Mon Aug 22 15:37:08 EST 1994

l_x = length(x);
l_y = length(y);
if l_x ~= l_y
  error('x and y must be the same length')
end

len = l_x;
x = reshape(x, 1, len);
y = reshape(y, 1, len);

%If the polynomial is not closed then do so.

if ( (x(1) ~= x(len)) | (y(1) ~= y(len)) )
  x = [x x(1)];
  y = [y y(1)];
  len = len + 1;
end

xmean = [ (x(2:len) + x(1:(len-1))) (x(1) + x(len)) ];
ydiff = [ (y(2:len) - y(1:(len-1))) (y(1) - y(len)) ];
z = abs(xmean*ydiff'/2);
