function H = Mcurvb(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
% returns mean curvature
%     H = Mcurvb(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%
%  use bicubicm to get k array's

[m,n] = size(k1);
onemat = ones(m,n);
fr = bicubicfr(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
fr2 = fr.*fr;
fc = bicubicfc(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
fc2 = fc.*fc;
frr = bicubicfrr(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
frc = bicubicfrc(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
fcc = bicubicfcc(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);

H = ((onemat+fc2).*frr + (onemat+fr2).*fcc - 2*fr.*fc.*frc)./  ...
    (2*sqrt(onemat+fr2+fc2).*sqrt(onemat+fr2+fc2).*sqrt(onemat+fr2+fc2));
