function vortexVectors(n,con,len)
% vortexVectors  Plots the velocity vectors for a potential flow vortex
%
% Synopsis:  vortex(n,con,len)
%
% Input:  n   = number of grid lines in each direction.  n should be
%               even to avoid division by zero.  Default: n=16
%         con = constant that determines the strength of the vortex
%               Default: con=1
%         len = size of domain:  -len <= x <= len;  -len <= y <= len
%               Default: len=25;
%
% Output:  contour plot of stream function

if nargin<1, n=16;   end
if nargin<2, con=1;  end
if nargin<3, len=4;  end

% --- Make matrices of x and y values;  quiver and surf plots are
%     easiest to work with when x and y values are matrices
[X,Y] = meshgrid(linspace(-len,len,nlinspace(-len,len,n)xi);

% --- Compute matrices of Cartesian velocity components U and V
R2 = X.^2 + Y.^2;   U = -con*Y./R2;   V = con*X./R2;

% --- Use the quiver function to plot velocity vectors as arrows
quiver(X,Y,U,V);   axis('equal','square');    xlabel('x');   ylabel('y');

% --- Open new figure window and make a surface plot of velocity magnitude
Vmag = sqrt( U.^2 + V.^2 );
figure;  surfc(X,Y,Vmag);
axis('equal','square');    xlabel('x');   ylabel('y');
if n>21,  shading('interp');  end   %  remove facet edges for fine meshes

% --- Draw vectors originating from points along the x and y axes
%     First, generate vectors with tails distributed along the x axis
x = linspace(-len,len,n);   y = zeros(size(x));   r2 = x.^2 + y.^2;
u= -con*y./r2;   v = con*x./r2;

% --- Next, generate vectors with tails distributed along the y axies
yy = linspace(-len,len,n);   xx = zeros(size(x));   r2 = xx.^2 + yy.^2;
uu= -con*yy./r2;   vv = con*xx./r2;

% --- Open a new figure window and plot the combined data sets
figure;    quiver([x xx],[y yy],[u uu],[v vv],0);
axis('equal','square');   xlabel('x');   ylabel('y');