% gradient_hemisphere generates a contour and "quiver" plot of the 
%  hemisphere height function h(x,y) and its gradient, following 
%  problem 7.5.1 in "Basic Training in Mathematics" by R. Shankar.
% 

% Programmer: Dick Furnstahl  furnstahl.1@osu.edu
%
% Revision history:
%   01-May-2005 -- original version for Physics 263
%

clear;    % this sets variables to zero, so we start with a clean slate

R = 1;    % set the radius of the hemisphere

% Set a grid X Y with the desired range and number of points
xmin = -1.1;  xmax = 1.1;   % set the minimum and maximum values of x
ymin = -1.1;  ymax = 1.1;   % set the minimum and maximum values of y
num_x = 25; num_y = 25;     % set the number of points in each direction
delta_x = (xmax-xmin)/(num_x-1);  % calculate the spacing in x
delta_y = (ymax-ymin)/(num_y-1);  % calculate the spacing in y

for i = 1:num_x                 % step through each of the x points
    x = xmin + (i-1)*delta_x;   % find the current value of x 
    for j = 1:num_y             % step through each of the y points
        y = ymin + (j-1)*delta_y;  % find the current value of y
        X(i,j) = x;                % set the (i,j) element of the X matrix
        Y(i,j) = y;                % set the (i,j) element of the Y matrix
        if ((x^2+y^2) < R^2)       % check if we're inside the circle 
           h(i,j) = (R^2 - x^2 - y^2)^(1/2);      % height function 
           hx(i,j) = -x/(R^2 - x^2 - y^2)^(1/2);  % x-component of gradient 
           hy(i,j) = -y/(R^2 - x^2 - y^2)^(1/2);  % y-component of gradient
        else
           % everything is zero outside
           h(i,j) = 0;
           hx(i,j) = 0;
           hy(i,j) = 0;
        end 
        hz(i,j) = 0;      % set z-component of gradient to zero for plotting
    end        
end

% Draw the two-D figure (we use "hold" to let us draw two plots on one figure)
figure(1);                               % give the figure a label
num_contours = 10;                       % use num_contours contour lines
contourf(X,Y,h,num_contours); colorbar;  % make the plot and add a color bar
hold on;                                 % hold the current plot
axis equal;                              % make the axes the same scale
scale = 5;                               % scaling of vector lengths
quiver(X,Y,hx,hy,scale);                 % plot the gradient vectors
xlabel('x-axis');  ylabel('y-axis');     % add labels
hold off;                                % release the hold

% Draw the three-D figure
figure(2);                               % give the figure a new label
surf(X,Y,h); colorbar;                   % make a surface plot with colorbar
hold on;                                 % hold the current plot
grid off;                                % turn the grid off
shading interp;                          % interpolate the colors 
xlabel('x-axis');  ylabel('y-axis'); zlabel('h(x,y)');     % add labels
axis equal;                              % make the axes the same scale
scale = 5;                               % scaling of vector lengths
quiver3(X,Y,h,-hx,-hy,hz,scale);         % plot minus the gradient vectors 
                                         %  horizontally at height h
hold off;                                % release the hole