% script to do a very simple "worm" movie, based on reading
% sub-tidal velocity fields form the model, and then by tracking
% drifters at a fixed depth. 
%
% the tracking algorithm works as follows:
%  1. obtain 2D velocity fields from closest time steps bracketing the 
%    current time.   These should be subtidal velocity fields. 
%  2. linearly interpolate to get the "now" velocity
%  3. interpolate the "now" velocity to the drifter positions using 
%    griddata. 
%  4. move the drifters forward in time by a small time step
%  5. back to step 1.

cdf='/snafu/d2/rps/kagom/ecom/charlie/ecom3d.cdf';
dt=2/24;   % interpolate velocity field every hour 
isub=3;   % start particle at every 3rd i and j point in model domain
dlong=4;  % drifter length in days (length of worm)
zlev=-5;  % level to track drifters
start=[1998 3 20 0 0 0];   % start tracking from this date
subframe=6;  % every 3rd step, update frame
i=sqrt(-1);
[d,x,y]=kslice(cdf,'depth');
[lon,lat]=merc2ll(x,y);
m_proj('transverse mercator','lon',[-72 -61],'lat',[38 47]);
earth_radius=6371.e3;
[x,y]=m_ll2xy(lon,lat);
x=x*earth_radius;   % convert to (approximate) meters
y=y*earth_radius;   % convert to (aprroximate) meters

igood=find(finite(d));
xgood=x(igood);
ygood=y(igood);

%id=1:length(igood);
%data=[id(:) xgood(:) ygood(:)];
%saveascii('test.node',data,'%d %10.1f %10.1f\n');

tri=delaunay(xgood,ygood);
%trisurf(tri,xgood,ygood,d(igood));view(2);shading('flat');dasp;

jd=ecomtime(cdf);
jstart=tind(jd,start);

dsub=d(1:isub:end,1:isub:end);
xsub=x(1:isub:end,1:isub:end);
ysub=y(1:isub:end,1:isub:end);
iwater=find(finite(dsub));
nframe=0;
t0=clock;
nt=round(jd(2)-jd(1))/dt;
dt=(jd(2)-jd(1))/nt;
dts=dt*24*3600;

nlength=round(dlong/dt);

i=sqrt(-1);
% initial particle positions

% set up the size of the drifter array
X=nan*ones(length(iwater),nlength);

X0=xsub(iwater)+i*ysub(iwater);
X(:,end-1)=X0;
set(gcf,'pos',[10 10 800 600]);
pslice(x,y,d,[0 300]);
c=jet(128);
colormap(c(20:100,:));
h=line(real(X)',imag(X)','color','black');

% old velocity
w1=zsliceuv(cdf,jstart,zlev);
for j=jstart:length(jd);
% new velocity
  w2=zsliceuv(cdf,j,zlev);
  for k=1:nt
     fac=(k-1)/nt;
     wnow=(1-fac)*w1+fac*w2;   % velocity at this time
% velocity at particle positions
     wX=triterp(tri,xgood,ygood,wnow(igood),...
         real(X(:,end-1)),imag(X(:,end-1)));  
% step forward in time
     X(:,end)=X(:,end-1)+dts*wX;
% if on land, put back to original position
     inan=find(isnan(X(:,end)));
     X(inan,end)=X0(inan);
     g=greg2str(gregorian(jd(j)+k*dt));
     title(['Time = ' g]);
% draw line segments
% set the xdata and ydata for each drifter 
     for jj=1:size(h);
      set(h(jj),'xdata',real(X(jj,:))','ydata',imag(X(jj,:))');
     end
     nframe=nframe+1;
     if (mod(nframe,subframe)==0),
       anim_frame('drifter',nframe);
     end
     X(:,1:end-1)=X(:,2:end);
     etime(clock,t0)
     t0=clock;
  end
  w1=w2;
end
anim_make