%                  Function t(riangle) Area 
%  This function calculates the area of a surface consisting of triangles
%  over a mesh with spacing dx and dy. The call is area(z,dx,dy,Z).
%  The matrix Z, (m+1) by (n+1), contains the heights of the surface
%  over the meshpoints.  The vector z contains all the "interior" 
%  elements of Z,arranged  as a row vector,
%  Z(i,j), 2 <= i <= m-1 and 2 <= j <= n-1.
%  The resulting function of z can be minimized by fmins.

function A = tarea(z,dx,dy,Z)
m = size(Z,1)-1;
n = size(Z,2)-1;
nn = n-1;
for i = 2:m
   Z(i, 2:n) = z((i-2)*nn +1: (i-1)*nn); 
end

  Z1 = diff(Z);
  Zvert = Z1(:, 1:n);
  Z2 = (diff(Z'))';
  Zhort = Z2(1:m, :);

  dx2 = dx^2; dy2 = dy^2;

  dAlower = .5*sqrt(dx2*dy2 + dx2*Zvert.^2 + dy2*Zhort.^2);
  Alower = sum(sum(dAlower));


  Zvert = Z1(:,2:n+1);
  Zhort = Z2(2:m+1,:);
  dAupper = .5*sqrt(dx2*dy2 + dx2*Zvert.^2 + dy2*Zhort.^2);
  Aupper = sum(sum(dAupper));

  A = Alower + Aupper;