function jquiver(x,y,uM,vM, varargin)
% jquiver(x,y,uM,vM, varargin)
%
% plots like squiver, but with arrowheads: you scale it.

%x,y are mx1 and nx1 vectors and uM, vM are mxn matrices 
% -- or --
%x,y,u,v are all nx1 vectors
%
%varargin is a coolection of name, value pairs:
%
%
% name                      what                   default 
% ------                   ------                 ---------
% scale            vector length as a fraction          1
%                     of one unit in x or y
% lthick           vector line thickness                1
% msize            marker size at (x,y)                 6
% head             arrowhead length is this 
%                     fraction of the longest 
%                     vector                            .05
% mcolor           marker color                         'g'
% vcolor           vector color                         'r'
% lcolor           line color (see 'addline')           'k'
% addline          plot (x,y) as a line                 0 
%                     (i.e. false; true is 1)
%
%
% NOTE: x and y are assumed to have the same scale.
%
% usage:
% jquiver(x, y, u, v, 'scale', 2, 'head', .01);

%% JH: 2001/09/14 - finally rewriting a_head from 1993

cfg.lthick = 1;
cfg.scale  = 1;
cfg.msize  = 6;
cfg.head   = .05;
cfg.mcolor = 'g';
cfg.vcolor = 'r';
cfg.lcolor = 'k';
cfg.addline = 0;
       
cfg = fillstruct(cfg, varargin);


[urows,ucols]=size(uM);      
[vrows,vcols]=size(vM);
x=x(:);y=y(:);

n=length(x);
m=length(y);


%if it's a cruisetrack, add it
if (cfg.addline)
   plot(x,y,'color',cfg.lcolor)
   hold on
end


%first make sure whether it's vectors or matrices:
if ((n~=m) | ((ucols>1) & (urows>1)) | ((vcols>1) & (vrows>1)) )
	%then it's matrices:

	[xM,yM]=meshgrid(x,y);
	lx=reshape((flipud(xM))',m*n,1);
	ly=reshape((flipud(yM))',m*n,1);
	lx=reshape(xM,m*n,1);
	ly=reshape(yM,m*n,1);
	lu=reshape(uM,m*n,1);
	lv=reshape(vM,m*n,1);
	
else %assume vectors
   
   lu=uM(:);lv=vM(:);lx=x(:);ly=y(:);	
   n1 = find(~isnan(uM) & ~isnan(vM));
   lu = lu(n1); lv = lv(n1); lx = lx(n1); ly = ly(n1);
	clear pi
	x=[lx lx+cfg.scale*lu]';y=[ly ly+cfg.scale*lv]';
   
   %plot the emination points
   h=plot(x,y,cfg.vcolor);set(h,'LineWidth',cfg.lthick)
	hold on
	h=plot(lx,ly,[cfg.mcolor,'.']);set(h,'MarkerSize',cfg.msize);
	

	headx=lx+cfg.scale*lu;
	heady=ly+cfg.scale*lv;
	uvangle=angle([lu + sqrt(-1)*lv]);% in radians
	
	d_angle=35*pi/180; %+- radians from uvangle to arrowhead angles
	a=sqrt(lu.^2+lv.^2);
   arrowlength=cfg.head*max(a);
	angle1=uvangle + pi - d_angle ;
	angle2=uvangle + pi + d_angle;
   
   u1=arrowlength.*cos(angle1);v1=arrowlength.*sin(angle1);
	u2=arrowlength.*cos(angle2);v2=arrowlength.*sin(angle2);
	
	x=[headx  headx+cfg.scale*u1]';y=[heady heady+cfg.scale*v1]';
	h=plot(x,y,cfg.vcolor);set(h,'LineWidth',cfg.lthick)
	x=[headx  headx+cfg.scale*u2]';y=[heady heady+cfg.scale*v2]';
   h=plot(x,y,cfg.vcolor);set(h,'LineWidth',cfg.lthick)
end
	
