% vector_fields generates a "quiver" plot of a two-dimensional
%  vector field with components Wx and Wy, and then approximately
%  calculates the (z-component of the) curl.
%
% The sample functions are taken in part from Problem 7.6.2
%  in "Basic Training in Mathematics" by R. Shankar.
% 

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

close all; clear;       % this 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) );

% This could be the velocity field of water in a bathtub: It has a curl!
 Wx = -Y;
 Wy = X;

% Problem 7.6.2 (i)
%$ Wx = Y.*sin(X);
%$ Wy = cos(Y);

% Problem 7.6.2 (ii)
%$ Wx = 2*X.*Y.^3;
%$ Wy = 3*Y.^2.*X.^2;

% Problem 7.6.2 (iii)
%$ Wx = sin(X);
%$ Wy = cos(Y);

% Problem 7.6.2 (iv)
%$ Wx = cosh(X)*cosh(Y);
%$ Wy = sinh(X)*sinh(Y);

% A radial function
%$ R = (X.^2 + Y.^2).^(1/2);   % R is the radial distance
%$ f = exp(-R);      % define some radial function (exponential here)
%$ Wx = f.*(X./R);   % define the matrix of x-components of the vector field
%$ Wy = f.*(Y./R);   % define the matrix of y-components of the vector field

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 
                                         % (we'll use the y-derivative dWxdy)
[dWydx dWydy] = gradient(Wy,dx,dy);      % find the gradient of Wy 
                                         % (we'll use the x-derivative dWydx)
curlW_z = dWydx - dWxdy;                 % z-component of the curl

% Draw the figures 
figure(1);                               % figure 1 plots the field itself
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');

figure(2);                               % figure 2 is an approximate curl
surf(X,Y,curlW_z); colorbar;             % make the plot and add a color bar
xlabel('x-axis');  ylabel('y-axis');     % add labels
title('Curl of W');