% divergence generates a "quiver" plot of a two-dimensional
%  vector field with components Wx and Wy, and then approximately
%  calculates the divergence.

% Programmer: Dick Furnstahl  furnstahl.1@osu.edu
%
% Revision history:
%   13-May-2005 -- original version, based on vector_fields.m

close all; clear;       % this closes plots and sets variables to zero

% Set a grid X Y with the desired range (and gridsize points on each axis)
xmin = -1;  xmax = 1;   % set the minimum and maximum values of x
ymin = -1;  ymax = 1;   % set the minimum and maximum values of y
gridsize = 10;
[X Y] = meshgrid( linspace(xmin,xmax,gridsize), linspace(ymin,ymax,gridsize) );

% Vector Field 1
Wx = X;  Wy = Y;

% Vector Field 2
%$ Wx = X.*Y.^2;  Wy = Y.*X.^2;

% Vector Field 3  [Problem 7.6.2 (iii)]
%$ Wx = sin(X);  Wy = cos(Y);

% Vector Field 4:  This has a curl, but does it have a divergence?
%$ Wx = -Y;  Wy = X;

% Vector Field 5
%$ Wx = X.*Y;  Wy = X.*Y;

% Vector Field 6  [Problem 7.6.2 (ii)]
%$ Wx = 2*X.*Y.^3;  Wy = 3*Y.^2.*X.^2;

dx = (xmax-xmin)/(gridsize-1);           % spacing between X points
dy = (ymax-ymin)/(gridsize-1);           % spacing between Y points
[dWxdx dWxdy] = gradient(Wx,dx,dy);      % find the gradient of Wx 
[dWydx dWydy] = gradient(Wy,dx,dy);      % find the gradient of Wy 

%curlW_z = dWydx - dWxdy;                % z-component of the curl
divW = dWxdx + dWydy;                    % the divergence!

% Draw the figures 
figure(1);                               % figure 1 plots the vector field 
scale = 2;                               % scaling of vector lengths
quiver(X,Y,Wx,Wy,scale);                 % plot the vector field W
xlabel('x-axis');  ylabel('y-axis');     % add labels
title('Vector Field W');                 % add a title

pause;                                   % pause between plots 

figure(2);                                   % contour plot of divergence
num_contours = 10;                           % set the number of contours
contourf(X,Y,divW,num_contours); colorbar;   % make the plot and add a color bar
figure(3);                                   % surface plot of divergence
surf(X,Y,divW); colorbar;                    % make the plot and add a color bar
xlabel('x-axis');  ylabel('y-axis');         % add labels
title('Divergence of W');                    % add a title