function []=vecplt3(w,wll,wind,windll,runoff,runoffll,jd)
% VECPLT3 animates vector currents on a map 
%  USAGE vecplt3(w,wll,wind,windll,runoff,runoffll,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(:,m+2)=runoff;
w=replace(w,NaN,0.);  %replace NaN with zeros to increase speed of plotting

wll(m+1,:)=windll;
wll(m+2,:)=runoffll;
%
%  get coastline and 50 m isobath
load coast
load bath50
%
%  set up plot
lat_mid=42.3;
lat_max=42.6;
lon_mid=-70.4;
dlon=1.4;
lat_fac=cos(lat_mid*pi/180);
dlat=dlon*lat_fac;
%
% factor for arrows
%
 xfac=.01;
 yfac=xfac*lat_fac;
%
[n,m]=size(w);
lon=wll(:,1);
lat=wll(:,2);
%  Do movie until user says no!
axis('square');
s='y';
while(s~='n')
 iclear=24*28;   % clear screen every 2 weeks
 greg1=gregorian(jd(1));
 for i=1:nsub:n
   if (rem(i-1,iclear)==0),  % clear screen every iclear hours
     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
     plot(bath50(:,2),bath50(:,1),'r.')
     u=zeros(1,m);
     v=zeros(1,m);
   end
   arrows(lon,lat,u,v,xfac,yfac,'i-');
   u=real(w(i,:));
   v=imag(w(i,:));
   arrows(lon,lat,u,v,xfac,yfac,'g-');
 end
 s=input('Do it again? [y or n]','s');
end
