function [gm samc] = mcurvature_vec(x,y,z)
% Description: The function calculates mean curvature and
% Surface Avg Mean Curvature (SAMC) of a surface formed by x, y & z.
% The input are the coordinate matrices x, y & z. The
% matrices can be formed using meshgrid or similar functions.
% This code is a vectorized form of original code (mcurvature) posted
% on Matlab File Exchange.
% The mean curvature is calculated according to the formula:
% If x:U->R^3 is a regular patch, then the mean curvature is given by
%
% H = (eG-2fF+gE)/(2(EG-F^2)),
%
% where E, F, and G are coefficients of the first fundamental form and
% e, f, and g are coefficients of the second fundamental form
% Reference: Gray, A. "The Gaussian and Mean Curvatures." 16.5 in
% Modern Differential Geometry of Curves and Surfaces with Mathematica,
% 2nd ed. Boca Raton, FL: CRC Press, pp. 373-380, 1997 (p. 377).
% Inspired by:
% Title: Mean Curvature
% Author: Ahmed Elnaggar
% Summary: Calculate the Mean curvature of a given surface (x,y,z).
gm = zeros(size(z));
[xu,xv] = gradient(x);
[xuu,xuv] = gradient(xu);
[xvu,xvv] = gradient(xv);
[yu,yv] = gradient(y);
[yuu,yuv] = gradient(yu);
[yvu,yvv] = gradient(yv);
[zu,zv] = gradient(z);
[zuu,zuv] = gradient(zu);
[zvu,zvv] = gradient(zv);
Xu(:,:,1) = xu;
Xu(:,:,2) = yu;
Xu(:,:,3) = zu;
Xv(:,:,1) = xv;
Xv(:,:,2) = yv;
Xv(:,:,3) = zv;
Xuu(:,:,1) = xuu;
Xuu(:,:,2) = yuu;
Xuu(:,:,3) = zuu;
Xuv(:,:,1) = xuv;
Xuv(:,:,2) = yuv;
Xuv(:,:,3) = zuv;
Xvv(:,:,1) = xvv;
Xvv(:,:,2) = yvv;
Xvv(:,:,3) = zvv;
E = dot(Xu,Xu,3);
F = dot(Xu,Xv,3);
G = dot(Xv,Xv,3);
m = cross(Xu,Xv,3);
temp(:,:,1) = sqrt(sum(m.*m,3));
temp(:,:,2) = temp(:,:,1);
temp(:,:,3) = temp(:,:,1);
n = m./temp;
L = dot(Xuu,n,3);
M = dot(Xuv,n,3);
N = dot(Xvv,n,3);
gm = ((E.*N)+(G.*L)-(2.*F.*M))./(2.*(E.*G - F.^2));
dim = size(z);
samc = 1/ (dim(1) * dim(2)) * sum (gm(:).^2);
