function handle = tdvec(x, b, t, d, xylimits, uscale, option1, option2, option3);
% function tdvec(x, b, t, d, xlimits, ylimits, option1, option2, option3)
%
% Plots a complex matrix as a field of arrows.
%
%  x*b      = u + i*v; row -> depth, column -> time
%  t        = row vector with time of each column in x*b
%  d        = column vector with depth of each row in x*b
%             (above 3 arrays are automatically transposed if necessary)
%  xylimits = [xmin xmax ymin ymax] for plot x- and y-axes
%  uscale   = velocity units per x-axis unit;
%             e.g.,  (m/s) per day
%  options:
%     (option1 to option3 can be any one of the following)
%     arrowlength = length of arrowhead in velocity units
%                   Note: set arrowlength = 0 to omit arrowheads;
%                         this greatly reduces memory requirements.
%     key         = [key_x key_y key_u];
%                   key_x, key_y are location of key in x,y units;
%                   key_u is length of key in velocity units
%     pos_kind    = [pos kind];
%                   pos  = 0  default, uses t as x-axis for stick plots
%                        > 1  no unit for x-axis and it uses 'pos'
%                             as its x-axis position;
%                   kind = 0  default, plot sticks
%                        = 1  for mean and trend (at the end_time)
%                        = 2  plot ellipse for inertial, semidiurnal
%                             and diurnal
%                        = 3  plot phase and direction
%
%           Eric Firing, Sat  08-19-1989
%
%  Modified the option part and added option pos_kind to plot ellipse,
%  phase and direction of inertial wave, semidiurnal and diurnal tides.
%
%            Willa Zhu,   05/30/91
% ---------------------------------------------------------------------------
%  EF: Fri  93/11/12; mod to take advantage of Matlab4
%  aspect ratio control.

axis_ratio =1/ 0.714;
xrange = abs(xylimits(2)-xylimits(1));
yrange = abs(xylimits(4)-xylimits(3));
data_ratio =  axis_ratio * (yrange/xrange);
uv_aspect  = data_ratio;
   % to mult by (v/uscale) to get depth units

sticks       = 0;
mean_trend   = 1;
plot_ellipse = 2;
phase_direct = 3;

uvc = x*b;

if (nargin < 6 | nargin > 9)
   error('ERROR: number of arguments should be 6 to 9'); end
arrowlength = 0;
key = [0 0 0];
pos_kind = [0 0];

if nargin >= 7
   if (length(option1) == 1)
      arrowlength = option1;
   elseif (length(option1) == 2)
      pos_kind = option1;
   elseif (length(option1) == 3)
      key = option1;
   end
end

if nargin >= 8
   if (length(option2) == 1)
      arrowlength = option2;
   elseif (length(option2) == 2)
      pos_kind = option2;
   elseif (length(option2) == 3)
      key = option2;
   end
end

if nargin == 9
   if (length(option3) == 1)
      arrowlength = option3;
   elseif (length(option3) == 2)
      pos_kind = option3;
   elseif (length(option3) == 3)
      key = option3;
   end
end

[nn,nt] = size(t);
if nn>nt, t=t';, nt=nn;, end;

[nd,nn] = size(d);
if nn>nd, d=d';, nd=nn;, end;

[m, n] = size(uvc);     % m is number of depth bins, n is number of times
if (m ~= nd | n ~= nt)
   if (m == nt & n == nd)
      uvc = uvc.';
      disp('input matrix has been transposed');
   else
      disp('ERROR: uvc dimensions should be length(d) rows b length(t) cols')
      return
   end
end

% Define corners of horizontal arrowhead when tip is at the origin.
if arrowlength > 0
   arrow1 = (-1 + i * 0.2) * arrowlength / uscale;
   arrow2 = (-1 - i * 0.2) * arrowlength / uscale;
end



if pos_kind(2) ~= 2 % ellipse
  if pos_kind(1) == 0
    tails = (ones(nd,1) * t)  + i * (d * ones(1,nt));
  else
    tails = (ones(nd,nt) * pos_kind(1))  + i * (d * ones(1,nt));
  end

  if pos_kind(2) == phase_direct
    alpha = atan2(imag(uvc(:,1)),real(uvc(:,1)));
    uvc(:,1) = 0.1*cos(alpha).*ones(size(uvc(:,1))).*real(uvc(:,1))./real(uvc(:,1)) ...
         + i * 0.1*sin(alpha).*ones(size(uvc(:,1))).*imag(uvc(:,1))./imag(uvc(:,1));
    signs = sign(abs(b(1,:))-abs(b(2,:)));
  end

  heads = (real(tails) + real(uvc) / uscale) ...
          + i * (imag(tails) + (imag(uvc) / uscale) * uv_aspect);
  directions = uvc ./ abs(uvc);

  if ((arrowlength > 0) & (pos_kind(2) ~= phase_direct))
    arrows1 =         real(arrow1 * directions)  ...
                + i * imag(arrow1 * directions) * uv_aspect;
    arrows2 =         real(arrow2 * directions)  ...
                + i * imag(arrow2 * directions) * uv_aspect;

    vectors = [tails(:) heads(:) ...
               (heads(:) + arrows1(:)) ...
               (heads(:) + arrows2(:)) heads(:) ...
               NaN * tails(:)].';
  else
    vectors = [tails(:) heads(:) NaN * tails(:)].';
  end
end

if pos_kind(2) == mean_trend
  end_ind = vectors(:,((nt-1)*nd+1):(nt*nd));
  handle = plot(end_ind(:));
elseif pos_kind(2) == sticks
  handle = plot(vectors(:));
elseif pos_kind(2) == plot_ellipse
% jr-+: replaced below with ellipse calculated as
% u + iv = b(1,:) * exp(i * omega * t) + b(2,:) * exp(-i * omega * t)
% where omega * t goes from 0 to 2*pi; this allows us to plot a smoother ellipse
% using a line rather than dots
%%  ii=length(vectors(2,:))/nt*[1:nt]-nd;
%%  c=[];
%%  for j=1:nd
%%   c =[c; vectors(2,ii+j)];
%%  end
%  c = reshape(vectors(2,:), nd, nt); % replaces above for loop
%  ellipse= conj(c');
  th = 2*pi*[(0:0.01:1) ].';
  ellipse = exp(i * th) * b(1,:) + exp(-i * th) * b(2,:);
  ellipse = [i * ones(length(th),1) * d' + pos_kind(1) + ...
             real(ellipse) / uscale + i * imag(ellipse) / uscale * uv_aspect; ...
             NaN * ellipse(1,:)];
  clear th
  handle = plot(ellipse(:));
elseif pos_kind(2) == phase_direct
  vectors = vectors(:,1:nd);
  handle = plot(vectors(:));
  hold on
  phase = signs'.*sin(alpha)*0.025 + i * (signs'.*cos(alpha)*0.025);
  direct_tail = vectors(2,:)';
  direct_head = (real(direct_tail) - real(phase)/uscale) ...
          + i * (imag(direct_tail) - imag(phase)/uscale*uv_aspect);

  if arrowlength > 0
     directions = phase ./ abs(phase);

     arrows1 =   -     real(arrow1 * directions(:,1) )  ...
                 - i * imag((arrow1) * directions(:,1)) * uv_aspect;
     arrows2 =   -     real(arrow2 * directions(:,1) )  ...
                 - i * imag((arrow2) * directions(:,1)) * uv_aspect;

     direct = [direct_tail(:) direct_head(:) ...
                 (direct_head(:) + arrows1(:)) ...
                 (direct_head(:) + arrows2(:)) direct_head(:) ...
                 NaN * direct_tail(:) ]';
  else
     direct = [direct_tail(:) direct_head(:) NaN * direct_tail(:)]';
  end
  handle = plot(direct(:));
else
  return
end

if key(3) ~= 0
   % Adjust key length to have only 1 digit after decimal.
   key(3) = (fix( key(3)*10.0 ) ) / 10.0;
   hold on
   plot([key(1) (key(1)+key(3)/uscale)], [key(2) key(2)]);
%   text(key(1)+key(3)/uscale+0.02*xrange, key(2)-0.025*yrange, ...
   text(key(1)+key(3)/uscale+0.01*xrange, key(2), ...
               sprintf('%3.1f m/s', key(3)));
   hold off
end

%axis([xylimits]);
set(gca, 'Xlim', xylimits(1:2), 'Ylim', xylimits(3:4));
if (v4)
   set(gca, 'AspectRatio', [axis_ratio, data_ratio]);
else
   set(gca, 'DataAspectRatio', [1 data_ratio 1]);
   set(gca, 'PlotBoxAspectRatio', [1 axis_ratio 1]);
end

