%% Chapter 5

clear all
close all

%
%% Three-Dimensional Graphs

flowfunction = @(x, y) y*(1 - 1/(x^2 + y^2))
figure; fsurf(flowfunction, [-2 2 -2 2])

%

figure; fsurf(flowfunction, [-1.5 1.5 -1.5 1.5])
axis([-1.5 1.5 -1.5 1.5 -5 5])
view([2 1 2])
xlabel('x'); ylabel('y'); zlabel('z')

%
%% Graphing Level Curves

list1=-2:0.5:2; figure; fcontour(flowfunction, [-2 2 -2 2],...
    'LevelList',list1)

figure; fcontour(flowfunction, [-2 2 -2 2],...
    'LevelList',list1, 'Fill', 'on')

%

syms y c; eqn = (flowfunction(0, y) == c)
solve(eqn, y)
list2=-2:0.1:2; figure; contourplot=fcontour(flowfunction, [-2 2 -2 2],...
    'LevelList',list2)
hold on; plot(0, 0, '.k', 'MarkerSize',500)

%
%% Partial Derivatives and the Gradient

syms x y; diff(x^3 + sin(x*y), x)

gradflowfunction = jacobian(flowfunction(x, y))

%
%% We'll use the output above to feed to the quiver command below.

[X, Y] = meshgrid(-2:0.2:2, -2:0.2:2);
U = (2.*X.*Y)./((X.^2 + Y.^2).^2);
V = (2.*Y.^2)./((X.^2 + Y.^2).^2) - 1./(X.^2 + Y.^2) + 1;
L = sqrt(U.^2 + V.^2);
figure; gradientfield=quiver(X, Y, U./L, V./L, 0.5)
axis equal tight

% Now we'll put the vector field and contour plot together
% Except I can't figure out how to do it with copyobj, Had to run both
% figures using hold on.
hold on
list1=-2:0.5:2; fcontour(flowfunction, [-2 2 -2 2],...
    'LevelList',list1)

%
%% Directional Derivatives

list3=0; figure; fcontour(flowfunction, [-2 2 -2 2],...
    'LevelList',list3,'MeshDensity',200); axis equal

%
%% Functions of Three or More Variables

% There is apparently no implicit 3D plotting outside of MuPad. Ughh!! What
% to do? Jon to the rescue with his function mfile "implicitplot3d. We use
% it here
% NO!! Use fimplicit3
syms x y z; h=z^2-x^2-y^2; 
figure; fimplicit3(h - 1, [-1.1 1.1 -1.1 1.1 -2 2]); axis equal

% According to the descriptive info at the head of the file, "The optional 
% parameter plotpoints (added at the end)
% gives the number of steps in each direction between plotting points.
% I needed to use a big number, like 100, to get a decent plot. It slowed 
% the computation down, but not too terribly.

% Next we superimpose the vector field. We use quiver3 and follow the schem
% laid out in the Matlab help for that command to generate the field.
[X,Y] = meshgrid(-1:0.5:1.,-1:0.5:1.);
Z = sqrt(1+ X.^2 + Y.^2);
[U,V,W] = surfnorm(X,Y,Z);
hold on;quiver3(X,Y,Z,U,V,W,2)
Z1 = - sqrt(1+ X.^2 + Y.^2);
[U,V,W] = surfnorm(X,Y,Z1);
quiver3(X,Y,Z1,U,V,W,2)