%% Chapter 6

% The Sphere
syms u v real; figure; fsurf(cos(u)*sin(v), sin(u)*sin(v), cos(v),...
    [0 2*pi 0 pi])
view([1, 1, 1])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;

%
%% Sinusoidal Cylinder

syms u v real; figure; fsurf(u, sin(u), v,...
    [-2*pi 2*pi -1 1])
view([1, 1, 1])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;

%
%% Cone

syms u v real; figure; fsurf((1 + v)*cos(u), (1 + v)*sin(u), 1 + v,...
    [0 2*pi -2 0])
view([1, 1, 1])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;

%
%% Monkey Saddle

syms u v real; figure; fsurf(u^3 - 3*u*v^2, [-2 2 -2 2])
view([1, 1, 1])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis off;

%
%% Torus

syms u v real; figure; 
fsurf((2 + cos(u))*cos(v), (2 + cos(u))*sin(v), sin(u), [0 2*pi 0 2*pi])
view([4, 4, 2])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;

%
%% Helicoid

syms u v real; figure; 
fsurf(v*cos(u), v*sin(u), u, [0 3*pi 0 3])
view([4, 0, 2])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;

%
%% Hyperboloid of One Sheet

syms u v real; figure
hyperboloid = [cos(u), sin(u), 0] + v*[-sin(u), cos(u), 1];
simplify(hyperboloid(1)^2 + hyperboloid(2)^2 - hyperboloid(3)^2)
fsurf(cos(u) - v*sin(u), sin(u) + v*cos(u), v, [0 2*pi, -2 2])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;

%
%% Figure Eight Curve
% Here's the code to generate the figure eight graph.

syms x y real
list=0; figure; fcontour(x^4-x^2+y^2, 'LevelList', list, 'LineColor','k')
axis([-1.1 1.1 -0.6 0.6])
hold on; syms x y real; fplot(x,sym('0'), [-1.1 1.1], 'Color', 'k')
fplot(sym('0'), y, [-0.5 0.5], 'Color', 'k')

%
%% Surface Graphs with Principal Directions and Normal
% First, negative curvature

syms x y real; figure; fsurf(x^2 - y^2, [-1 1 -1 1])
view([1 -2 1.5]);
hold on; syms t; fplot3(t, sym(0), t^2, [-1, 1], 'LineWidth', 5,'Color', 'k')
fplot3(sym(0), t, -t^2, [-1, 1], 'LineWidth', 5,'Color', 'k')
fplot3(sym(0), t, sym(0), [-1/2 1/2], 'LineWidth', 3,'Color', 'k')
fplot3(t, sym(0), sym(0), [-2 1.2], 'LineWidth', 3,'Color', 'k')
%fplot3(sym(0), sym(0), t, [0 3/4], 'LineWidth', 3, 'Color', 'k')
text(0,0,1.2, 'U'), text(-0.2, 3/5, 0, 'T_1'), text(1.5, 0,0, 'T_2')
quiver3(0,0,0,0,0,1, 'Color', 'k', 'LineWidth', 2)
text(0,-0.1,.95, '\Lambda')
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off; 

% Then positive curvature
syms x y real; figure; fsurf(-x^2 - y^2, [-1 1 -1 1])
view([1 -2 1.5]);
hold on; syms t; fplot3(t, sym(0), -t^2, [-1, 1], 'LineWidth', 5,'Color', 'k')
fplot3(sym(0), t, -t^2, [-1, 1], 'LineWidth', 5,'Color', 'k')
fplot3(sym(0), t, sym(0), [-1/2 1/2], 'LineWidth', 3,'Color', 'k')
fplot3(t, sym(0), sym(0), [-1 1.2], 'LineWidth', 3,'Color', 'k')
%fplot3(sym(0), sym(0), t, [0 3/4], 'LineWidth', 3, 'Color', 'k')
text(0,0,1.2, 'U'), text(-0.2, 3/5, 0, 'T_1'), text(1.5, 0,0, 'T_2')
quiver3(0,0,0,0,0,1, 'Color', 'k', 'LineWidth', 2)
text(0,-0.1,.95, '\Lambda')
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;

%
%% Geodesic Triangle

syms u v real; figure; fsurf(cos(u)*sin(v), sin(u)*sin(v), cos(v),...
    [0 pi/2 0 pi/2])
view([1, 1, 1])
title(''); xlabel(''); ylabel(''); zlabel('')
grid off; axis equal; axis off;
text(0,0,1.05, 'P'), text(1.05,0,0, 'Q'), text(0,1.05,0, 'R')

%
%% Surfaces with Tangent Planes

syms x y real; figure; fsurf(x^2 - y^2, [-1 1 -1 1])
view([1 -2 1.5]); hold on;;
fsurf(0, [ -1/2 1/2 -1/2 1/2])

syms x y real; figure; fsurf(-x^2 - y^2, [-1 1 -1 1])
view([1 -2 1.5]); hold on;;
fsurf(0, [ -1/2 1/2 -1/2 1/2])

%
%% Curvature Calculations with Matlab

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 specific computations of various curvatures for several surfaces

syms u v real; sphere = [cos(u)*sin(v), sin(u)*sin(v), cos(v)];
principalcurves(sphere, u, v)
simplify(ans, 'Steps', 30)

gausscurv(sphere, u, v)
simplify(ans, 'Steps', 30)

meancurv(sphere, u, v)
simplify(ans, 'Steps', 30)

hyperboloid = [cos(u), sin(u), 0] + v*[-sin(u), cos(u),1];
principalcurves(hyperboloid, u, v)
simplify(ans, 'Steps', 50)

gausscurv(hyperboloid, u, v)
simplify(ans, 'Steps', 30)

helicoid = [v*cos(u), v*sin(u), u];
principalcurves(helicoid, u, v)
simplify(ans, 'Steps', 30)

meancurv(helicoid, u, v)
% Note: didn't have to simplify this one.