function [hp, ht] = m_vec(s, x, y, varargin)  % and color, etc
% M_VEC Draws fancy arrows/quiverplots on a map.
%    [HP, HT]=M_VEC(S,LONG,LAT,VARARGIN) draws arrows as a patch 
%    object on a map created by the M_Map package.
%      HP: handle to the patch object
%      HT: handle to text object if this is a key; see below.
%      S:  scale factor, arrow data units per inch.
%          It is based on the figure PaperPosition, so be
%          sure to set this (e.g. via "landscape" or "portrait")
%          before calling m_vec.
%      LONG,LAT: longitude, latitude of the arrows
%          LONG and LAT must have the same dimensions, but can be
%          scalars or vectors; if scalars, multiple arrows can
%          be plotted at a single location.
%      VARARGIN can consist of any of the following:
%        U,V   U,V,C   Z,U,V,C
%      followed by optional arrow parameters,
%      followed by optional patch parameters.
%        U,V are vectors containing the east and north components
%                 of the arrows.
%        C is an optional colorspec for all arrows, or an
%                 array of CData, one value per arrow.
%                 Defaults to black.
%        Z is a height in axes data units: this is subject
%                 to future modification or omission, and
%                 is probably not useful now as-is.
%      optional arrow parameters: keyword-value pairs, shown
%                 here with default values:
%         'headangle',60     degrees: angle of arrow tip
%         'headwidth',NaN    points: direct specification of
%                                width, instead of headangle
%         'headlength',5     points: length of tip; set to 0
%                                to omit arrowhead entirely
%         'shaftwidth',1     points: width of arrow shaft
%         'centered', 'no'   'yes' to make x,y the arrow
%                                center instead of its tail
%         'key', ''          make a labelled horizontal arrow
%                                if the string is not empty;
%                                then the string labels the
%                                arrow, and the second argument
%                                returned, ht, is the handle of
%                                the string.
%          'edgeclip', 'off'  If 'on' then arrows IN the axes
%                             are clipped if their heads are
%                             OUT of the axes.
%
%      optional patch parameters: any valid patch properties
%         may be specified here; they are passed directly
%         to the patch function.
%
%    M_VEC called without any parameters generates a demonstration plot.

%
%  Mon  98/02/16 Eric Firing, efiring@soest.hawaii.edu
%
% 6/Nov/00 - eliminate returned stuff if ';' neglected (thx to D Byrne)
% 2/May/01 - small bug fix (thx to Pierre Jaccard)
% 7/jun/06 - arrows near boundaries were not done correctly - now fixed.
% 12/Dec/12 - added clipping for arrows at boundaries ('edgeclip')

global MAP_PROJECTION

% Have to have initialized a map first


if nargin==0,  % demo

    % demvec.m
    % This is a demonstration of m_vec.
    % Sun  98/02/22 Eric Firing

    % Set up the figure and axes at the start, so that the vector lengths
    % will come out right when the figure is printed:

  %  orient tall
    % Main axes, for the map.
    ha1 = axes;
    pos = get(ha1,'position');
    pos1 = pos;
    pos1(4) = pos1(4) - 0.1;
    pos1(2) = pos1(2) + 0.1;
    set(ha1,'position',pos1)

    % Colormap axes.  Make them smaller than for the default colormap.
    pos2 = [0.3  0  0.4  0];
    pos2(2) = pos(2);
    pos2(4) = 0.02;
    % This positioning of the colormap is actually not very good as
    % printed, but I am not going to fiddle with it any more now.
    % This is just a typical difficulty with normalized coordinates
    % combined with fixed DataAspectRatio; if
    % something looks reasonable on the screen, it is likely to look
    % wrong when printed.

    ha2 = axes('position', pos2);

    axes(ha1); % Back to the main axes.

    m_proj('ortho','lat',48','long',-123', 'rad', 10, ...
        	'rec', 'off' );
    m_coast('patch',[0.9 0.95 0.9]);
    m_grid('linestyle','-','xtick',[-135:5:-110],'linewi',2);
    title('Demonstration of m\_vec')


    %% The following form, centered and without heads, is suitable for
    %% the square roots of the principle axes of variance ellipses,
    %% for example.  Color is given as a single RGB triplet; if it were
    %% given as an array of triplets, the two lines could be different
    %% colors; unfortunately, as of Matlab 5.1, if you do this the plot
    %% will be forced into Zbuffer mode with the warning:
    %%  "RGB color data not yet supported in Painter's mode"
    %% This example also illustrates specification of a black patch
    %% boundary, via direct patch property specifications following any
    %% vector parameter options.
    hpv1 = m_vec(100, [-133 -133], [49 49], [0 50], [100 0.0],...
	  [0.7 0.8 0.9],'centered','yes', ...
	  'shaftwidth', 5, 'headlength', 0,...
	  'EdgeColor','k');

    %% Here are three vectors, no color specified, all properties defaulted.
    hpv2 = m_vec(100, [-128 -128 -128], [46 46 46], ...
     [0 25*sqrt(2) 50], [50 25*sqrt(2) 0]);

    incs = (1:20)/20;
    vlat = 42 + incs*3;
    vlon = -127 - incs*2;
    uu = 50*sin(incs*2*pi);
    vv = 50*cos(incs*2*pi);

    %% This simulates a set of vectors such as currents measured
    %% with an ADCP along a cruise track.  Color is specified by letter,
    %% in this case. The arrows and heads are about as skinny as they
    %% can reasonably be.
    hpv3 = m_vec(100, vlon, vlat+2, uu, vv, 'm', ...
     'shaftwidth', 0.2, 'headlength', 2.5);

    %% Now let's make a similar vectors, but with colors based on the
    %% colormap, simulating SST, for example:
    vlon = vlon - 3;
    sst = 12 - incs * 4;
    hpv4 = m_vec(100, vlon-4, vlat, uu, vv, sst);
    colorbar(ha2);
    axes(ha2); xlabel('SST'); axes(ha1);
    %% These vectors will not show up well on the screen because they
    %% have no edges, but they will print adequately.  There does not
    %% seem to be any easy way to get around this; I would have to
    %% completely change the way the patches are specified.  Ordinarily,
    %% one would probably use fatter than default vectors with this
    %% method anyway, so that the colors are clear when printed; and in
    %% this case the screen display will look OK also.

    % Here I show the edgeclip property
    hpv4 = m_vec(100, vlon+4, vlat-5, uu, vv, sst,'edgeclip','on');


    %% Key: Note that it can be on or off the map. The vector and
    %% text are the same color by default. Handles to both the vector
    %% (patch) and the label (text) are returned so that you can
    %% change the font, color, etc.
    [hpv5, htv5] = m_vec(100, -115, 38, 50, 0, 'b', 'key', '50 cm s^{-1}');
    set(htv5,'FontSize',8);

    return

end











if isempty(MAP_PROJECTION),
  error('No Map Projection initialized - call M_PROJ first!');
end;



% Default arrow parameters:
centered = 0;
headlength = 5/72;
headwidth  = NaN;
headangle = 40;
shaftwidth = 1/72;
c = 'k';
key = '';

clip = 'on'; % except for the key
edgeclip = 'off';  % for arrows at the edge

if nargin < 5,   % Minimum: s,x,y,u,v
   help('mvec');
   error('not enough input arguments');
end

UVperIn = s;
x = x(:);
y = y(:);

%% Begin the slightly complicated parsing of arguments.  This
%% could probably be done more efficiently and elegantly.
% One cause of complexity is that the optional argument c
% could be either character or numeric.
keyvars = {};  % This will hold keyword/value pairs.
nvarargin = length(varargin);
% Find the index of the first string argument in varargin.
istr0 = 0;
for ii = 1:nvarargin
   if ischar(varargin{ii}), istr0 = ii; break, end
end

if istr0 > 0,    % There are strings.
   if rem(nvarargin - istr0, 2) == 0,  % an odd number
      c = varargin{istr0}; % should be a colorspec string
      istr0 = istr0 + 1;
   end
   % istr0: start of keyword/value pairs.
   n_numeric = istr0 - 1;  %Actually, numeric or colorspec.
   keyvars = varargin(istr0:nvarargin);
   finished = 0;
   while ~isempty(keyvars) & ~finished,
      kv = lower(keyvars{1});
      value = keyvars{2};
      if     strcmp(kv, 'headlength')
         headlength = value/72;
      elseif strcmp(kv, 'headwidth')
         headwidth = value/72;
      elseif strcmp(kv, 'headangle')
         headangle = value;

      elseif strcmp(kv, 'shaftwidth')
         shaftwidth = value/72;
      elseif strcmp(kv, 'centered')  & lower(value(1)) == 'y',
         centered = 1;
      elseif strcmp(kv, 'key')
         key = value;
         clip = 'off'; % Can put key outside the map.
      elseif strcmp(kv, 'edgeclip')
         edgeclip = value;
      else
         finished = 1;  % no match; break out
      end
      if ~finished,
         keyvars(1:2) = [];
      end
   end
else
   n_numeric = nvarargin;
end

% Calculate the headwidth if it is not given explicitly:
if isnan(headwidth) & headangle < 170 & headangle > 0
   headwidth = headlength * tan(headangle*pi/180);
end
headwidth = max([headwidth; shaftwidth]);

if n_numeric == 2 | n_numeric == 3,
   u = varargin{1}(:);
   v = varargin{2}(:);
   if n_numeric == 3,
      c = varargin{3};
   end
   z = zeros(size(u));
elseif n_numeric == 4,
   z = varargin{1}(:);
   u = varargin{2}(:);
   v = varargin{3}(:);
   c = varargin{4};
else
   help('mvec')
   error('not enough numeric arguments')
end

[nr,nc] = size(c);
if nr == 1 & nc == length(u) & (nc ~= 3 | (any(c<=1) | any(c>=0))),
   c = c(:);
end
% c could be a 1x3 colorspec

if (length(x) == 1 & length(y) == 1 & length(u) > 1)
   x = x(ones(size(u)));
   y = y(ones(size(u)));
end

%% End of input argument parsing.

OrigAxUnits = get(gca,'Units');
if OrigAxUnits(1:3) == 'nor'
   OrigPaUnits = get(gcf, 'paperunits');
   set(gcf, 'paperunits', 'inches');
   figposInches = get(gcf, 'paperposition');
   set(gcf, 'paperunits', OrigPaUnits);
   axposNor = get(gca, 'position');
   axWidLenInches = axposNor(3:4) .* figposInches(3:4);
else
   set(gca, 'units', 'inches');
   axposInches = get(gca, 'position');
   set(gca, 'units', OrigAxUnits);
   axWidLenInches = axposInches(3:4);
end

% Multiply inches by the following to get data units:
scX = diff(get(gca, 'XLim'))/axWidLenInches(1);
scY = diff(get(gca, 'YLim'))/axWidLenInches(2);
sc = max([scX;scY]);  %max selects the dimension limited by
                      % the plot box.

Width = shaftwidth*sc;
HeadWidth = headwidth*sc;
HeadLength = headlength*sc;

uvmag = abs(u + i*v);
% Arrow lengths in plot data units:
L = uvmag*sc/UVperIn;

% base of arrows. Don't plot if outside boundaries (except for keys for which clip is off)
[xs, ys] = m_ll2xy(x,y, 'clip', clip);

%[xsp, ysp] = m_ll2xy(x+0.1*u./uvmag, y+0.1*v.*cos(y*pi/180)./uvmag, 'clip', clip);
% Vector angles in the Cartesian data-unit system of m_map:
%Ang = angle( (xsp-xs) + i*(ysp-ys) );
% ABove angle calc could fail when arrows were near boundaries. Replace with this.
% Now we calculate angles. Use a small offset, and keep clipping off to prevent odd things
% happening.
%    - RP 7/Jun/06
[xsp, ysp] = m_ll2xy([x x+0.00001*u./uvmag]',[y y+0.00001*v.*cos(y*pi/180)./uvmag]' , 'clip', 'off');
Ang = angle( diff(xsp)' + i*diff(ysp)' );

if ~isempty(key), Ang = 0; end


nvec = length(L);
Zero = zeros(nvec,1);
One  = ones(nvec,1);

% Normal arrow dimensions:
HL = Zero+HeadLength;
HW = Zero+HeadWidth;
W =  Zero+Width;

% Distinguish zero-length vectors from non-zero:
mm = (L < 100*eps);
i_zero = find(mm);
i_nonzero = find(~mm);
% Don't plot if length is zero.
if ~isempty(i_zero)
   HL(i_zero) = NaN; HW(i_zero) = NaN; W(i_zero) = NaN;
end

if ~isempty(i_nonzero)
   ii = i_nonzero;
   if HeadLength == 0,   %% square end; no arrowhead
      W(ii)  = Width;
      HW(ii) = Width;
      HL(ii) = Zero(ii);  % Thanks for Pierre Jaccard for this fix
   else
      % If the arrow length is less than the headlength,
      % omit the arrow shaft and just plot a head scaled
      % to the length.
      i_short = ii( find(L(ii) < HeadLength) );
      W(i_short)  = 0;
      HL(i_short) = L(i_short);
      HW(i_short) = HL(i_short) * (HeadWidth/HeadLength);
   end
end


% Just a change of variable names for historical reasons:
Y = ys;
X = xs;
% It is not clear whether non-zero elevations will be useful;
% that depends on whether the entire m_map system will work
% with 3-D views, as in making a surface plot of topography.
% Then one might want to have a current profile represented
% as a stack of vectors at appropriate heights above the topog.
Z = z(:);


Xzero = Zero;
if centered,
   Xzero = -L/2;
end


nV = 7*nvec;  % number of Vertices: 7 per arrow.
Vert = zeros(nvec,3);

Vert(1:7:nV,:) =  [Xzero,  W/2, Z ];       % back corner of shaft
Vert(2:7:nV,:) = Vert(1:7:nV,:) + ...
                 [(L-HL), Zero,  Zero];    % shaft-head junction
Vert(3:7:nV,:) = Vert(2:7:nV,:) + ...
                 [Zero, (HW-W)/2,  Zero];  % point of barb
Vert(4:7:nV,:) = [Xzero + L, Zero, Z];     % tip of arrow

%% This could be done more efficiently with fancier indexing, but
%% it is probably not a bottleneck, hence not worth the trouble.

% Reflect the top half to get the bottom half.
% First replicate:
Vert(5:7:nV,:) = Vert(3:7:nV,:);
Vert(6:7:nV,:) = Vert(2:7:nV,:);
Vert(7:7:nV,:) = Vert(1:7:nV,:);
% Then negate y to reflect:
Vert(5:7:nV,2) = -Vert(5:7:nV,2);
Vert(6:7:nV,2) = -Vert(6:7:nV,2);
Vert(7:7:nV,2) = -Vert(7:7:nV,2);

% Make an index array for operating on all vertices of each vector:
ii = (1:nvec);
ii = ii(ones(7,1),:);
ii = ii(:);

%% Rotate:
Vxy = exp(i*Ang(ii)).*(Vert(:,1) + i*Vert(:,2));

%% Translate:
Vxy = Vxy + X(ii) + i*Y(ii);

%%
Vert(:,1) = real(Vxy);
Vert(:,2) = imag(Vxy);


Faces = [1:nV].';            %Top
Faces = reshape(Faces,7,nvec).';

% Extremely narrow patches don't show up on the screen (although they seem
% to be printed OK) when EdgeColor is 'none', so when the arrows are all
% the same color, set the EdgeColor to be the same as FaceColor.
% Set clip off here so arrows are complete - RP 7/Jun/06

% Request to clip arrows at the plot edge.

if strcmp(edgeclip,'off'),
 hp = patch('Faces', Faces, 'Vertices', Vert, 'tag', 'm_vec','clipping','off');
else,
  [LG,LN]=m_xy2ll(reshape(Vert(:,1),7,nvec),reshape(Vert(:,2),7,nvec)); % Converts vertices in 7 point lines (i.e. columns) in lat/long
  [X,Y]=m_ll2xy(LG,LN  ,'clip','patch');                                % Converts back to x/y, but does clipping on for each column
  hp = patch('Faces', Faces, 'Vertices', [X(:) Y(:)], 'tag', 'm_vec','clipping','off');
end;

if ischar(c) | (size(c,1) == 1 & size(c,2) == 3),
   set(hp, 'EdgeColor', c, 'FaceColor', c, 'LineWidth', 0.1);
else
   set(hp, 'EdgeColor', 'none', 'FaceColor','Flat', ...
     'FaceVertexCdata', c);
end
if ~isempty(keyvars)
   set(hp, keyvars{:});
end

if ~isempty(key)
   ht = text(X(1), Y(1)-0.5*HW, Z(1), key, ...
      'color', c, 'horiz','left','vert','top', 'tag','m_vec', ...
      'clipping','off');
   set(hp,'clipping','off')
else
   ht = [];
end

if nargout==0
  clear hp ht
end;