function h = swan_vectorplot(xax,yax,mag,wdir,dx,dy,S,rot)
% SWAN_VECTORPLOT - Plots a field of vectors from SWAN
% h = swan_vectorplot(xax,yax,mag,wdir,dx,dy,S,rot)
% xax, yax  - x- and y-axis vectors
% mag, wdir - magnitude and direction of vector field
% dx, dy    - number of points to skip
% S         - scale factor for arrow lengths (start with 1)
% rot       - degrees to rotate vector field (geographic convention)
% (needs xycoord.m)

% Chris Sherwood, USGS
% Last revised June 20, 2000

if(nargin<7),rot=90;,end;
[m,n]=size(wdir);
wdir=reshape(wdir,m*n,1);
% this line takes SWAN Cartesian coord and converts to arrows pointing
% in the direction that wind/waves blow/propagate towards, in geographic conv.
wdir = -1 * wdir + rot;	% rotate so 0=North, 90=E, 180=S, 270=W
ilo = find(wdir < 0.);
if length(ilo) > 0, wdir(ilo) = wdir(ilo) + 360; end
wmag =  ones(m*n,1).*reshape(mag,m*n,1);
[u,v]=xycoord( wmag, wdir );
baddir = find( (abs(v)./abs(u))> 2 );
u(baddir)=zeros(size(baddir));
v(baddir)=zeros(size(baddir));
u=reshape(u,m,n);
v=reshape(v,m,n);
[x,y]=meshgrid(xax,yax);
% subsample vector field
if(dy > 1 | dx >1),
  x=x(fix(dy/2):dy:m,fix(dx/2):dx:n);
  y=y(fix(dy/2):dy:m,fix(dx/2):dx:n);
  u=u(fix(dy/2):dy:m,fix(dx/2):dx:n);
  v=v(fix(dy/2):dy:m,fix(dx/2):dx:n);
end
h=quiver(x,y,u,v,S,'-k');
