%% Problem Set F, Problem 10, Part (e)

% We begin by reproducing the commands that set up the curvature
% computations that we use for the surfaces in this problem, and of course
% for all the others in the chapter and problem set.

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));

% Now for the two surfaces -- Scherk's and Enneper's, both of which are, as we shall see, minimal surfaces.
%
%% (i) Scherk's Minimal Surface

close all
syms u v real
scherk = [u, v, log(cos(v)/cos(u))]
meancurv(scherk, u, v)
figure; fsurf(scherk(1), scherk(2), scherk(3), [(-pi/2)+0.01 (pi/2)-0.01 ...
    (-pi/2)+0.01 (pi/2)-0.01])

%
%% (ii) Enneper's Minimal Surface
% We'll give two views -- the default, and another in which we can see the "hole" better.

enneper = [u - u^3/3 + u*v^2, v - v^3/3 + u^2*v, u^2 - v^2] 
simplify(meancurv(enneper, u, v))
figure; fsurf(enneper(1), enneper(2), enneper(3), [-2 2 -2 2])
figure; fsurf(enneper(1), enneper(2), enneper(3), [-2 2 -2 2])
view([3 0 1])