function []=vecplt6(w,wll,wind,windll,jd)
% VECPLT6 steps through vector currents on a map 
%         using mouse buttons as controls
%
%  Does hard copy
%      right mouse button advances in time
%      middle mouse quits
%      left mouse button goes back in time
%      "p" prints out plot
%
%  click farther right on the plot to take greater time steps
%
% plots surface, mid and deep as different colors
% 0-10 m red
% 10-30 m green
% 30-60 m blue
% > 60 m  yellow
%
%
%  USAGE vecpltp(w,wll,wind,windll,jd)
%        w  matrix where each column is a complex velocity time series
%        wll two-column matrix [lon lat] for velocity stations
%        w  complex wind vector
%        run real runoff vector
%
%        if it's too slow, try increasing nsub
%
nsub=1;
[n,m]=size(w);
w(:,m+1)=wind;
w=replace(w,NaN,0.);  %replace NaN with zeros to increase speed of plotting

wll(m+1,:)=windll;
%
%  get coastline and 40 m isobath
load coast
load bath40
%[nbath,foo]=size(bath40);
%bath40=bath40(1:4:nbath,:);
%
%  set up plot
lat_mid=42.22;
lon_mid=-70.55;
dlon=1.1;
lat_fac=cos(lat_mid*pi/180);
dlat=dlon*lat_fac;
%
% factor for arrows
%
 yfac=24*3600/1.e5/111;   % length particles would travel in one day :
 xfac=yfac/lat_fac;
%
% 111 km/ deg lat;
[n,m]=size(w);
lon=wll(:,1);
lat=wll(:,2);
%  Do movie until user says no!
axis('square');
s='y';
iw=m;   %wind
is=[1 3 6 12 14 16 19 21 24 ];   %surface 4-8 m
im=[2 4 7 13 15 17 20 22 25 ];   % mid 18-28 m
inb=[5 8 9  23 ];        %  33-60 m        % near bot
ib=[10 11 18  ];           % > 60
ibb=[2  5  9 11 13 15  20  23 ]; % 10 m or less from  bottom
while(s~='n')
 iclear=24*28;   % clear screen every 2 weeks
 greg1=gregorian(jd(1));
 i=nsub;
 mb=1
 k=0;
 while(mb~=0),
   k=k+1;
   [x,y,mb]=ginput(1,'sc');
   if(mb==112),  % print out plots if user hits "p"
     clg
hold off
axis('square')
xmin=lon_mid-.5*dlon;
xmax=lon_mid+.5*dlon;
ymin=lat_mid-.5*dlat;
ymax=lat_mid+.5*dlat;
xx = [xmin xmax xmax xmin; xmax xmax xmin xmin];  % x coords of box
yy = [ymin ymin ymax ymax; ymin ymax ymax ymin];  % y coords of box
axis([xmin xmax ymin ymax])          % set axis limits
hold on
     !rm brad.met
%    for ii=1:16;    % how many plots to output?
     for ii=1:1;
     greg=gregorian(jd(i));
     tit=greg2str(greg);
     title(tit(1:15))
     hold on
     sca=-12;
     u=real(w(i,:));
     v=imag(w(i,:));
     i=i+1;
     arrow3(lon(is),lat(is),u(is),v(is),sca,-2,1,xfac,yfac,'r-')
     arrow3(lon(im),lat(im),u(im),v(im),sca,-2,1,xfac,yfac,'g-')
     arrow3(lon(ibb),lat(ibb),u(ibb),v(ibb),sca,-2,1,xfac,yfac,'b-');
     arrow3(lon(iw),lat(iw),u(iw),v(iw),sca,-2,1,xfac,yfac,'-');
     meta brad
     end
     !matlab2lpr brad
   else
     mb=mb-2;
     if(mb==0),return,end;
     i=i+50*x*mb*nsub;
     i=rem(i,n);
     if(i<0),i=i+n,end;
   if (rem(k-1,iclear)==0 | mb==97),  % clear screen every iclear button hits
     hold off
     axis([lon_mid-.5*dlon lon_mid+.5*dlon lat_mid-.5*dlat lat_mid+.5*dlat]);
     plot(coast(:,1),coast(:,2));
     greg2=gregorian(jd(min(i+iclear,n)));
     title(date_lbl(greg1,greg2));
     greg1=greg2;
     hold on
     mplot(bath40,'r.')
     u=zeros(1,m);
     v=zeros(1,m);
   end
   end
   greg2=gregorian(jd(i));
   str1=greg2str(greg2)
   arrows(lon(is),lat(is),u(is),v(is),xfac,yfac,'i-');
   arrows(lon(im),lat(im),u(im),v(im),xfac,yfac,'i-');
   arrows(lon(inb),lat(inb),u(inb),v(inb),xfac,yfac,'i-');
   arrows(lon(ib),lat(ib),u(ib),v(ib),xfac,yfac,'i-');
   arrows(lon(iw),lat(iw),u(iw),v(iw),xfac,yfac,'i-');
   u=real(w(i,:));
   v=imag(w(i,:));
   arrows(lon(is),lat(is),u(is),v(is),xfac,yfac,'c2-');
   arrows(lon(im),lat(im),u(im),v(im),xfac,yfac,'c3-');
   arrows(lon(inb),lat(inb),u(inb),v(inb),xfac,yfac,'c4-');
   arrows(lon(ib),lat(ib),u(ib),v(ib),xfac,yfac,'c5-');
   arrows(lon(iw),lat(iw),u(iw),v(iw),xfac,yfac,'c6-');
 end
end
