%% Surface Invariants
% These are the commands for computing the geometric invariants that
% characterize a parametrized surface.

% In the following, surf = [x(u, v), y(u, v), z(u, v)] is a parametrized 
% surface. For each such surface, the following program computes its 
% geometric invariants.

veclength = @(vec) sqrt(dot(vec, vec));
normalvector = @(surf, u, v) simplify(cross(diff(surf, u), diff(surf, v)));
tangentplane = @(surf, u, v, x, y, z) dot([x, y, z] - ...
    surf, normalvector(surf, u, v)== 0);
unitnorm = @(surf, u, v) normalvector(surf, u, v)/... 
    veclength(normalvector(surf, u, v));

E = @(surf, u, v) dot(diff(surf, u), diff(surf, u));
F = @(surf, u, v) dot(diff(surf, u), diff(surf, v));
G = @(surf, u, v) dot(diff(surf, v), diff(surf, v));
e = @(surf, u, v) dot(unitnorm(surf, u, v), diff(surf, u, 2));
f = @(surf, u, v) dot(unitnorm(surf, u, v), diff(surf, u, v));
g = @(surf, u, v) dot(unitnorm(surf, u, v), diff(surf, v, 2));

shapeoperator = @(surf, u, v) ((E(surf,u,v)*G(surf,u,v)-F(surf,u,v)^2)^(-1)) ...
   *[[e(surf,u,v)*G(surf,u,v)-f(surf,u,v)*F(surf,u,v), f(surf,u,v)*E(surf,u,v) - e(surf,u,v)*F(surf,u,v)]; ...
   [f(surf,u,v)*G(surf,u,v) - g(surf,u,v)*F(surf,u,v), g(surf,u,v)*E(surf,u,v) - f(surf,u,v)*F(surf,u,v)]];

principalcurves = @(surf, u, v) eig(shapeoperator(surf, u, v));
gausscurv = @(surf, u, v) det(shapeoperator(surf, u, v));
meancurv = @(surf, u, v) (1/2)*trace(shapeoperator(surf, u, v));