function dirfield(dx,dy,x,y,s)
%DIRFIELD	Direction field plot.
%	DIRFIELD(Z) draws a graph that displays the angle of the 
%	complex elements of Z as line elements centred on equally
%	spaced grid positions, one element per matrix element.
%
%	DIRFIELD(DX,DY) is equivalent to DIRFIELD(DX+i*DY).  It displays
% the direction field where line elements have slope dy/dy for
% corresponding elements of	matrices DX and DY.
%
%	DIRFIELD(Z,X,Y) and DIRFIELD(DX,DY,X,Y) uses X and Y as the
% position vectors for the x and y axis respectively. X and Y
% should be monotonic and linearlly spaced for best results.
%
%	DIRFIELD(Z,'S'), DIRFIELD(X,Y,'S') and DIRFIELD(DX,DY,X,Y,'S')
% use line style 'S' where 'S' is any legal linetype as described
% under the PLOT command.
%
%	For example, try
%	   [x,y] = meshdom(-2:.2:2, -2:.2:2);
%	   z = x .* exp(-x.^2 - y.^2);
%	   [px,py] = gradient(z,.2,.2);
%	   contour(z),hold on, dirfield(px,py), hold off
%
%	See also QUIVER.

% Modified from MatLab's 'quiver' command to produce 
% direction field plots by Colin Fox, 18 Feb 1991.
% Modifies for MATLAB 5 4 June 1997 by Michelle Wood.

xx = [-0.5 0.5].';
yy = [0 0].';
element = xx + yy.*sqrt(-1);

if nargin == 1      % called as dirfield(z)
   dy = imag(dx); dx = real(dx);
   [m,n] = size(dx); x = 1:n; y = 1:m;
   s = 'r-';
elseif nargin == 2  % dirfield(z,s) or dirfield(dx,dy)
   if isstr(dy)     % dirfield(z,s)
      s = dy;
      dy = imag(dx); dx = real(dx);
   else             % dirfield(dx,dy)
      s = 'r-';
   end
   [m,n] = size(dx); x = 1:n; y = 1:m;
elseif nargin == 3  % dirfield(z,x,y) or dirfield(dx,dy,s)
   if isstr(x)      % dirfield(dx,dy,s)
      s = x;
      [m,n] = size(dx); x = 1:n; y = 1:m;
   else             % dirfield(z,x,y)
disp('I''d say **dirfield(z,x,y)**')
      y = x; x = dy;
      dy = imag(dx); dx = real(dx);
      s = 'r-';
   end
elseif nargin == 4  % dirfield(dx,dy,x,y) or dirfield(z,x,y,s)
   if isstr(y)      % dirfield(z,x,y,s)
      s = y;
      y = x; x = dy;
      dy = imag(dx); dx = real(dx);
   else             % dirfield(dx,dy,x,y)
      s = 'r-';
   end
end

[xx,yy] = meshgrid(x, y);
grid = xx + yy.*sqrt(-1); grid = grid(:);
dx = dx(:); dy = dy(:);
z = (dx + dy.*sqrt(-1)).'; z = z./(abs(z) + eps);
minstep = min([min(diff(x)), min(diff(y))]);
eltlength = 0.5*minstep;
a = eltlength * element * z + ones(size(element)) * grid.';
%keyboard
plot(real(a), imag(a), s);