% Function tsurf (triangular surface)
% [X,Y] is a meshgrid of points in rectangle R. The heights of
% the surface at the meshpoints are given in a matrix Z.
% Frequently Z will be computed by a function Z = f(X,Y).
% The call tsurf(X,Y,Z) calculates the area of the surface
% consisting of triangles spanning the points (xj,yi,f(xj,yi)).
% The surface consisting of triangles is also graphed.
function out = tsurf(X,Y,Z)
n = size(X,2)-1; m = size(Y,1)-1;
a = X(1,1); b = X(1,n+1); c = Y(1,1); d = Y(m+1,1);
dx = (b-a)/n; dy = (d-c)/m;
A = 0;
for i = 1:m
for j = 1:n
x1 = X(i,j); x2 = x1+dx;
y1 = Y(i,j); y2 = y1+dy;
z1 = Z(i,j); z2 = Z(i,j+1);
z3 = Z(i+1, j+1); z4 = Z(i+1,j);
xedge = [x1 x2 x1];
yedge = [y1 y1 y2];
zedge = [z1 z2 z4];
fill3(xedge,yedge, zedge, 'b')
hold on
xedge = [x2 x2 x1];
yedge = [y1 y2 y2];
zedge = [z2 z3 z4];
fill3(xedge, yedge, zedge, 'c')
dAlower = .5*sqrt(dx^2*dy^2 + (z2-z1)^2*dy^2 + (z4-z1)^2*dx^2);
dAupper = .5*sqrt(dx^2*dy^2 + (z4-z3)^2*dy^2 + (z2-z3)^2*dx^2);
A = A + dAlower + dAupper;
end
end
out= A;
hold off