...

function [A,Qx,Qy,Ix,Iy,Ixy,xbar,ybar]= ...
         prop(x,y)
%
% [A,Qx,Qy,Ix,Iy,Ixy,xbar,ybar]=prop(x,y)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes the geometrical 
% properties for a polygon. 
%
% x      - vector of x coordinates
% y      - vector of x coordinates
%
% A      - area of polygon
% Qx     - first moment of area about 
%          y-axis
% Qy     - first moment of area about 
%          x-axis
% Ix     - inertia about y axis
% Iy     - inertia about x axis
% Ixy    - product of inertia
% xbar   - x distance to centroid
% ybar   - y distance to centroid
%
% User m functions called:  none
%----------------------------------------------

i=1:length(x); j=i+1;
A=0; Qx=0; Qy=0; Ix=0; Iy=0; Ixy=0;
x=[x(:);x(1)]; y=[y(:);y(1)];
xi=x(i); yi=y(i); xj=x(j); yj=y(j);
Ai=xi.*yj-xj.*yi; A=sum(Ai)/2;
Qx=sum((yi+yj).*Ai)/6; Qy=sum((xi+xj).*Ai)/6;
Ix=sum((yi.^2+yi.*yj+yj.^2).*Ai)/12;
Iy=sum((xi.^2+xi.*xj+xj.^2).*Ai)/12;
Ixy=sum((xi.*yi+0.5*(xi.*yj+xj.*yi)+xj.*yj) ...
         .*Ai)/12;
xbar=Qy/A; ybar=Qx/A;
%...Reverse signs if necessary (when nodes
%   traversed clockwise)
if A < 0
  A=-A; Qx=-Qx; Qy=-Qy; Ix=-Ix; Iy=-Iy; Ixy=-Ixy;
end