In reply to Steven Chipps.

The area is computed as a local variable inside the function. You could return it as part of the output of _Centroid, but it probably makes more sense to create a new function:

/* Area: Return the area of a simple polygon.
P is an N x 2 matrix of (x,y) values for consectutive vertices. N > 2. */
start Area(P);
lagIdx = 2:nrow(P) || 1; /* vertices 2, 3, ..., N, 1 */
xi = P[ ,1]; yi = P[ ,2]; /* x_i and y_i for i=1..N */
xip1 = xi[lagIdx]; yip1 = yi[lagIdx]; /* x_{i+1} and y_{i+1}, x_{N+1}=x_1 */
d = xi#yip1 - xip1#yi; /* vector of difference of products */
A = 0.5 * sum(d); /* signed area of polygon */
return( A );
finish;

Area = Area(L);
print Area;