function hh = quiver_jvm(x, y, u, v, rescale, length_legend, tail_point, ...
    axis_val, arrow_head, alpha, beta, col_arrows, col_background)
% QUIVER_JVM: like quiver but with extra features.
%   quiver_jvm rescales vectors, has a legend and maintains the magnitude and
%   direction of the original vectors. Additionally, any cases where the
%   magnitude of the vector is zero or a NaN will not be printed - this
%   avoids the printing of a dot on the paper.
%
% hh = quiver_jvm(x, y, u, v, rescale, length_legend, tail_point, axis_val, ...
%            arrow_head, alpha, beta, col_arrows, col_background)
%
%  INPUT ARGUMENTS.
%
% The first 4 arguments are the same as in a call to quiver.
% RESCALE: A vector (RESCALE, 0) in (u, v) space will be plotted as length 1
%          in the x direction. If AXIS_VAL (see below) is also passed
%          scaling will be carried out so that all vectors of the same
%          magnitude in (u, v) space will have the same length on the screen
%          or paper. Note that (x, y) space is not the same as space on the
%          paper and if AXIS_VAL is not passed then all vectors of the same
%          magnitude in (u, v) space will have the same length in (x, y)
%          space - which is not what you usually want. If it is [] or a NaN
%          then a guess is made for rescale. (Actually this is
%             rescale = 10*sqrt(max(u(:).^2 + v(:).^2))/ ...
%             max(abs([max(x(:)) - min(x(:)); max(y(:)) - min(y(:))]));
% LENGTH_LEGEND: If +ve this will be the length in (u, v) space for the arrow
%                in the legend. If 0 or -ve then no legend will be drawn. If
%                it is [] or NaN or no argument is passed then the maximum
%                vector will be used.
% TAIL_POINT: point in (x, y) space for the tail of the arrow in the
%             legend. If no argument is passed or any element of tail_point
%             is a NaN then the user is prompted to use the mouse to mark
%             the tail of the arrow.
% AXIS_VAL: [xmin xmax ymin ymax] as passed to the axis command to specify
%           the figure axis. Passing this vector also causes the vector
%           lengths to be scaled so that vectors of the same magnitude in
%           (u, v) space will have the same length on the screen or paper.
%           If AXIS_VAL is not passed or an empty matrix or NaN is passed
%           then all vectors of the same magnitude in (u, v) space will have
%           the same length in (x, y) space - which is often not what you
%           want.
% ARROW_HEAD: == 1 then the usual open arrow heads will be used (like QUIVER)
%             == 2 then filled arrow heads will be used. The default is 1.
% ALPHA: Size of arrow head relative to the length of the vector (default
%        0.33).
% BETA: Width of the base of the arrow head relative to the length (default
%        0.33).
% COL_ARROWS: The colour of the arrows (default black).
% COL_BACKGROUND: The colour of the background to the legend (default white)
%
%  OUTPUT ARGUMENTS.
%
% HH: a vector of line handles. Thus you can do "set(h, 'color', [0 0 0])" to
% change the vectors and legend from the default blue to black.
%
%       Notes on the magnitude and direction of the vectors:
%
% Consider a common problem - plotting a series of wind vectors on a lat/lon
% grid. Assume that the wind at point A is a south-westerly. Furthermore,
% consider a point B which is 10 degrees north and 10 degrees east of A. Using
% quiver, the standard matlab routine for plotting vectors, the wind vector
% with its tail at A will point directly at B. This will happen even if the
% plot is stretched in one direction, i.e., delta(lon) ~= delta(lat) and does
% not allow for the non-linear mapping between lon/lat space and physical space
% (we are missing a cos(lat)). Accordingly, the vector need not be at 45
% degrees to the paper or screen. Furthermore, winds of the same magnitude but
% different direction will appear on the plot with different magnitudes (and
% the wrong directions) due to differential stretching in the lon and lat
% direction. By using the AXIS_VAL argument quiver_jvm can show a plot in which
% the magnitudes and angles of each vector are "correct" (relative to the
% paper). Of course, this means that the south-westerly vector described above
% will not point from A to B. Although this is usually what is required it
% can be misleading when plotting something like pressure contours with
% geostrophic winds overlaid. In this case it is best to just avoid the
% problem by fiddling things so that the delta(lon) == delta(lat). The
% easiest way to do this is by using 'axis equal'. Alternatively you can
% choose various settings to get the same effect. This can be confusing,
% especially because crucial axes lengths may be stored in normalised
% co-ordinates. Hence the relevant bit of code is:
%
% >>  gca_units_old = get(gca, 'units');
% >>  set(gca, 'units', 'inches');
% >>  Position = get(gca, 'Position')
% >>  new_scale = (axis_val(4) - axis_val(3))/(axis_val(2) - axis_val(1))* ...
% >>      Position(3)/Position(4);
% >>  set(gca, 'units', gca_units_old) 
% >>  v = v*new_scale;
%
% new_scale is the value used in quiver_jvm for the rescaling of the vector
% components and so the aim is to choose the settings so that new_scale == 1.
% The default settings have:
% >> Position(3)/Position(4) == 1.2679
% and so 'axis equal' simply sets 
% (axis_val(4) - axis_val(3))/(axis_val(2) - axis_val(1)) == 1/1.2679
%
% All of this implies that if you are going to set the sizes of the figures
% and axes (by set(gcf, 'Position'), set(gcf, 'PaperPosition'),
% set(gca, 'Position')) you should do these before calling quiver_jvm.

% $Id: quiver_jvm.m,v 1.15 2000/09/28 07:15:44 mansbrid Exp $
% Copyright J. V. Mansbridge, CSIRO, Mon Mar  9 18:24:39 EST 1998

% Arrow head parameters
sym = '';

% Set default values

if nargin == 4
  rescale = NaN;
  length_legend = NaN;
  tail_point = NaN;
  scale_y = 0;
  arrow_head = 1; % Plot arrows
  alpha = 0.33; % Size of arrow head relative to the length of the vector
  beta = 0.33;  % Width of the base of the arrow head relative to the length
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 5
  length_legend = NaN;
  tail_point = NaN;
  scale_y = 0;
  arrow_head = 1; % Plot arrows
  alpha = 0.33; % Size of arrow head relative to the length of the vector
  beta = 0.33;  % Width of the base of the arrow head relative to the length
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 6
  tail_point = NaN;
  scale_y = 0;
  arrow_head = 1; % Plot arrows
  alpha = 0.33; % Size of arrow head relative to the length of the vector
  beta = 0.33;  % Width of the base of the arrow head relative to the length
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 7
  scale_y = 0;
  arrow_head = 1; % Plot arrows
  alpha = 0.33; % Size of arrow head relative to the length of the vector
  beta = 0.33;  % Width of the base of the arrow head relative to the length
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 8
  arrow_head = 1; % Plot arrows
  alpha = 0.33; % Size of arrow head relative to the length of the vector
  beta = 0.33;  % Width of the base of the arrow head relative to the length
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 9
  alpha = 0.33; % Size of arrow head relative to the length of the vector
  beta = 0.33;  % Width of the base of the arrow head relative to the length
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 10
  beta = 0.33;  % Width of the base of the arrow head relative to the length
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 11
  col_arrows = [0 0 0];
  col_background = [1 1 1];
elseif nargin == 12
  col_background = [1 1 1];
elseif nargin ~= 13
  error('quiver_jvm needs between 4 and 13 arguments')
end

% set scale_y according to axis_val

if nargin >= 8
  if isempty(axis_val) | isnan(axis_val)
    axis_val = axis;
    scale_y = 0;
  else
    scale_y = 1;
  end
end
  
% Set arrow style stuff

if isempty(arrow_head)
  arrow_head = 1;
elseif isnan(arrow_head)
  arrow_head = 1;
elseif (arrow_head ~= 1) & (arrow_head ~= 2)
  error('arrow_head must equal 1 or 2')
end
if isempty(alpha)
  alpha = 0.33; % Size of arrow head relative to the length of the vector
elseif isnan(alpha)
  alpha = 0.33; % Size of arrow head relative to the length of the vector
end
if isempty(beta)
  if arrow_head == 1
    beta = alpha;
  else
    beta = 0.3*alpha;
  end
elseif isnan(beta)
  if arrow_head == 1
    beta = alpha;
  else
    beta = 0.3*alpha;
  end
end

% Set colors for arrows and background in legend

if isempty(col_arrows)
  col_arrows = [0 0 0];
elseif isnan(col_arrows)
  col_arrows = [0 0 0];
end
col_legend = col_arrows;
if isempty(col_background)
  col_background = [1 1 1];
elseif isnan(col_background)
  col_background = [1 1 1];
end

nextpl_gcf = get(gcf, 'NextPlot');
ax = gca;
nextpl_gca = get(ax, 'NextPlot');
% xlim = get(gca, 'XLim');
% ylim = get(gca, 'YLim');
han = [];

vel_max = sqrt(max(u(:).^2 + v(:).^2));
if isempty(length_legend) | isnan(length_legend)
  length_legend = vel_max;
end

if scale_y == 1
  % Temporarily change units to absolute ones to enable the calculation of
  % new_scale, the stretching value in the y direction.
  
  % Position = get(gca, 'Position')
  % PaperPosition = get(gcf, 'PaperPosition')
  % new_scale = (axis_val(4) - axis_val(3))/(axis_val(2) - axis_val(1))* ...
  %     Position(3)*PaperPosition(3)/(Position(4)*PaperPosition(4));
   
  gca_units_old = get(gca, 'units');
  set(gca, 'units', 'inches');
  Position = get(gca, 'Position');
  new_scale = (axis_val(4) - axis_val(3))/(axis_val(2) - axis_val(1))* ...
      Position(3)/Position(4);
  set(gca, 'units', gca_units_old) 
  v = v*new_scale;
else
  new_scale = 1;
end
if isempty(rescale) | isnan(rescale)
  rescale = 10*vel_max/ ...
	    max(abs([max(x(:)) - min(x(:)); max(y(:)) - min(y(:))]));
end
extra = [alpha beta arrow_head scale_y new_scale];
han = quiver_old(x, y, u/rescale, v/rescale, 0, extra);

% Change the colour of the arrows.

hxx = findobj(han, 'Type', 'line');
set(hxx, 'Color', col_arrows)
hxx = findobj(han, 'Type', 'patch');
set(hxx, 'FaceColor', col_arrows)
set(hxx, 'EdgeColor', col_arrows)

hold on
ax_lims = axis;
if exist('axis_val') == 1
  vert_ratio = (axis_val(4) - axis_val(3))/(ax_lims(4) - ax_lims(3));
else
  vert_ratio = 1;
end

% Do legend stuff

if length_legend > 0
  length_arrow = length_legend/rescale;
  draw_legend = 1;
else
  draw_legend = 0;
end

if draw_legend
  if isempty(tail_point) | any(isnan(tail_point))
    disp('click mouse at desired tail point of the arrow in the legend'); 
    tail_point = ginput(1);
  end

  % Write the text to find out the height of the letters; this enables us to
  % draw a patch of the appropriate size. We then delete the text since the
  % patch will be drawn over it anyway and then draw the text we again.
  
  h_text = text(tail_point(1)+0.5*length_arrow, ...
		tail_point(2), ['  ' num2str(length_legend)]);
  set(h_text, 'HorizontalAlignment', 'Center');
  set(h_text, 'VerticalAlignment', 'bottom');
  set(h_text, 'Color', col_legend)
  extent = get(h_text, 'Extent');
  height_text = extent(4);
  delete(h_text)
  delta_x = length_arrow*0.2;
  delta_y = max([length_arrow*new_scale*0.2 height_text*vert_ratio]);
  h_p = patch(tail_point(1) + length_arrow*[0 1 1 0] + delta_x*[-1 1 1 -1], ...
	      tail_point(2) + delta_y*[-1 -1 1 1], col_background);
  set(h_p, 'EdgeColor', col_legend)
  h_text = text(tail_point(1)+0.5*length_arrow, ...
		tail_point(2), ['  ' num2str(length_legend)]);
  set(h_text, 'HorizontalAlignment', 'Center');
  set(h_text, 'VerticalAlignment', 'bottom');
  set(h_text, 'Color', col_legend)
  if arrow_head == 1 % open arrowheads
    xx = tail_point(1) + length_arrow*[0 1 NaN 1 1-alpha NaN 1 1-alpha];
    yy = tail_point(2) + new_scale*length_arrow* ...
	 [0 0 NaN 0 -alpha*beta NaN 0 alpha*beta];
    h_pl = plot(xx, yy);
    set(h_pl, 'Color', col_legend)
  elseif arrow_head == 2 % filled arrowheads
    xx = tail_point(1) + length_arrow*[0 1];
    yy = tail_point(2)*[1 1];
    h_pl = plot(xx, yy);
    set(h_pl, 'Color', col_legend)
    xx = tail_point(1) + length_arrow*[1 1-alpha 1-alpha];
    yy = tail_point(2) + new_scale*length_arrow*[0 -alpha*beta alpha*beta];
    h_p = patch(xx, yy, col_legend);
    set(h_p, 'EdgeColor', col_legend)
  end
end

hold off
if scale_y == 1
  axis(axis_val)
end
set(gcf, 'NextPlot', nextpl_gcf);
set(ax, 'NextPlot', nextpl_gca);

if nargout > 0
  hh=han;
end

% The following is a copy of the matlab routine quiver except that I have
% added code to avoid the printing of any vector whose magnitude is zero or a
% NaN. There is also an extra input argument. The new
% 6th argument is EXTRA. It contains values used to specify the
% characteristics of the arrow heads.

function hh = quiver_old(varargin)
%QUIVER Quiver plot.
%   QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
%   at the points (x,y).  The matrices X,Y,U,V must all be the same size
%   and contain corresponding position and vecocity components (X and Y
%   can also be vectors to specify a uniform grid).  QUIVER automatically
%   scales the arrows to fit within the grid.
%
%   QUIVER(U,V) plots velocity vectors at equally spaced points in
%   the x-y plane.
%
%   QUIVER(U,V,S) or QUIVER(X,Y,U,V,S) automatically scales the 
%   arrows to fit within the grid and then stretches them by S.  Use
%   S=0 to plot the arrows without the automatic scaling.
%
%   QUIVER(...,LINESPEC) uses the plot linestyle specified for
%   the velocity vectors.  Any marker in LINESPEC is drawn at the base
%   instead of an arrow on the tip.  Use a marker of '.' to specify
%   no marker at all.  See PLOT for other possibilities.
%
%   QUIVER(...,'filled') fills any markers specified.
%
%   H = QUIVER(...) returns a vector of line handles.
%
%   Example:
%      [x,y] = meshgrid(-2:.2:2,-1:.15:1);
%      z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.15);
%      contour(x,y,z), hold on
%      quiver(x,y,px,py), hold off, axis image
%
%   See also FEATHER, QUIVER3, PLOT.

%   Clay M. Thompson 3-3-94
%   Copyright (c) 1984-97 by The MathWorks, Inc.
%   $Revision: 1.15 $  $Date: 2000/09/28 07:15:44 $

% Arrow head parameters
alpha = 0.33; % Size of arrow head relative to the length of the vector
beta = 0.33;  % Width of the base of the arrow head relative to the length
autoscale = 1; % Autoscale if ~= 0 then scale by this.
arrow_head = 1; % Plot arrows
sym = '';

filled = 0;
ls = '-';
ms = '';
col = '';

nin = nargin;
% Parse the string inputs
while isstr(varargin{nin}),
  vv = varargin{nin};
  if ~isempty(vv) & strcmp(lower(vv(1)),'f')
    filled = 1;
    nin = nin-1;
  else
    [l,c,m,msg] = colstyle(vv);
    if ~isempty(msg), 
      error(sprintf('Unknown option "%s".',vv));
    end
    if ~isempty(l), ls = l; end
    if ~isempty(c), col = c; end
    if ~isempty(m), ms = m; arrow_head = 0; end
    if isequal(m,'.'), ms = ''; end % Don't plot '.'
    nin = nin-1;
  end
end

error(nargchk(2,6,nin));

% Check numeric input arguments
if nin<4, % quiver(u,v) or quiver(u,v,s)
  [msg,x,y,u,v] = xyzchk(varargin{1:2});
else
  [msg,x,y,u,v] = xyzchk(varargin{1:4});
end
if ~isempty(msg), error(msg); end

if nin==3 | nin==5, % quiver(u,v,s) or quiver(x,y,u,v,s)
  autoscale = varargin{nin};
elseif nin == 6
  autoscale = varargin{5};
  extra = varargin{6};
  alpha = extra(1); % Size of arrow head relative to the length of the vector
  beta = extra(2); % Width of the base of the arrow head relative to the length
  arrow_head = extra(3); % Plot arrows
  scale_y = extra(4);
  new_scale = extra(5);
end

% Scalar expand u,v
if prod(size(u))==1, u = u(ones(size(x))); end
if prod(size(v))==1, v = v(ones(size(u))); end

if autoscale,
  % Base autoscale value on average spacing in the x and y
  % directions.  Estimate number of points in each direction as
  % either the size of the input arrays or the effective square
  % spacing if x and y are vectors.
  if min(size(x))==1, n=sqrt(prod(size(x))); m=n; else [m,n]=size(x); end
  delx = diff([min(x(:)) max(x(:))])/n;
  dely = diff([min(y(:)) max(y(:))])/m;
  len = sqrt((u.^2 + v.^2)/(delx.^2 + dely.^2));
  autoscale = autoscale*0.9 / max(len(:));
  u = u*autoscale; v = v*autoscale;
end

ax = newplot;
next = lower(get(ax,'NextPlot'));
hold_state = ishold;

% Make velocity vectors
x = x(:).'; y = y(:).';
u = u(:).'; v = v(:).';

% Extra code to eliminate those vectors whose magnitude is zero or a NaN.

mag = u.*u + v.*v;
ff = find(~isnan(mag) & mag > 0);
if length(ff) > 0
  x = x(ff);
  y = y(ff);
  u = u(ff);
  v = v(ff);
end
uu = [x;x+u;repmat(NaN,size(u))];
vv = [y;y+v;repmat(NaN,size(u))];

h1 = plot(uu(:),vv(:),[col ls]);

if arrow_head == 1 % open arrowheads
  % Make arrow heads and plot them
  hu = [x+u-alpha*(u+beta*(v+eps)/new_scale);x+u; ...
        x+u-alpha*(u-beta*(v+eps)/new_scale);repmat(NaN,size(u))];
  hv = [y+v-alpha*(v-beta*(u+eps)*new_scale);y+v; ...
        y+v-alpha*(v+beta*(u+eps)*new_scale);repmat(NaN,size(v))];
  hold on
  h2 = plot(hu(:),hv(:),[col ls]);
elseif arrow_head == 2 % filled arrowheads
  % Make arrow heads and plot them
  hu = [x+u-alpha*(u+beta*(v+eps)/new_scale);x+u; ...
        x+u-alpha*(u-beta*(v+eps)/new_scale)];
  hv = [y+v-alpha*(v-beta*(u+eps)*new_scale);y+v; ...
        y+v-alpha*(v+beta*(u+eps)*new_scale)];
  hold on
  h2 = patch(hu, hv, 'b');
else
  h2 = [];
end

if ~isempty(ms), % Plot marker on base
  hu = x; hv = y;
  hold on
  h3 = plot(hu(:),hv(:),[col ms]);
  if filled, set(h3,'markerfacecolor',get(h1,'color')); end
else
  h3 = [];
end

if ~hold_state, hold off, view(2); set(ax,'NextPlot',next); end
if nargout>0, hh = [h1;h2;h3]; end
