addpath code
clear all
close all

%--------------------
%PARAMETERS
%--------------------

dirOut	='scratch/';
netcdfDir='jplroms/forecast/';
gridFile='jplroms/mb_das_2012012403.nc';

fSaveGrid=0;
fSaveMat=1;
fSaveTxt=1;
fSavePng=1;
fSaveFig=1;
fPlot=1;

dt		= 0.25/24;	%timestep for advection (days)
zpark		= -1e3;		%max parking depth (meters)
parkTime	= 17*3600;	%time to park underwater
telemetryTimes	= [0 14]*60;	%telemetry periods	
profileSpeed	= 0.08;		%profile speed (m/s)
deploymentLength= 3;		%total length of the deployment	(days)

startTime=datenum('2/1/2012 00:00');
xy0=[-122 36.8;-122.4 36];
x0=xy0(:,1);y0=xy0(:,2);

%--------------------
%load grid
%--------------------

%run ./start_make_grid

load('jplroms/grid.mat','X','Y','Z','zcoord','wc','belowBottomIdx');
%for some reason interp2 only likes trnasposed matrices
Xt=X';Yt=Y';Zt=Z';

%initial parking depth of a drifter
z0=interp2(Xt,Yt,Zt,x0,y0);
%z0=max(zpark,z0);
z0=zeros(size(x0));

%--------------------
%predict drifter location
%--------------------

inTLine=startTime+[0:dt:deploymentLength]';

[U,V,t]=readJplUV(netcdfDir,inTLine,[size(X), length(zcoord)]);
U(isnan(U))=0;
V(isnan(V))=0;

%loop through floats
xadv=[];yadv=[];zadv=[];tadv=inTLine;
for ifloat=1:length(x0)
  [xadv(ifloat,:) yadv(ifloat,:) zadv(ifloat,:) action]=...
      advectParticles3dfloat2(Xt,Yt,Zt,zcoord,U,V,inTLine,x0(ifloat),y0(ifloat),z0(ifloat),...
	zpark,parkTime,telemetryTimes,profileSpeed);
%    plot3(xadv(ifloat,:), yadv(ifloat,:), zadv(ifloat,:),'.-');grid on
end

%--------------------
%save results 
%--------------------

fnameOut=['_' datestr(min(double(tadv)),'yyyymmdd') '-' datestr(max(double(tadv)),'yyyymmdd')];

if fSaveMat
  save(fullfile(dirOut ,['traj' fnameOut '.mat']),'xadv', 'yadv', 'zadv', 'tadv')
end

if fSaveTxt
  for ifloat=1:size(xadv,1)
    tmpOut=double([tadv xadv(ifloat,:)' yadv(ifloat,:)' zadv(ifloat,:)']);
    save(fullfile(dirOut,['traj' fnameOut '_float' sprintf('%03d',ifloat) '.txt']),'-ascii','tmpOut')
  end
end

%--------------------
%plot results
%--------------------
if fPlot
  for ifloat=1:size(xadv,1)
      figure(1);clf
      contour(X,Y,Z,[-100 -200 -400 -1000],'k')
      colormap jet
      hold on
      for ifloat=1:size(xadv,1)
        plot(xadv(ifloat,1:end),yadv(ifloat,1:end),'k')
        scatter(xadv(ifloat,1:end),yadv(ifloat,1:end),5,tadv(1:end)-tadv(1),'filled')
      end
      plot(wc.LON,wc.LAT,'k')
      hold off
      colorbar
      axis([min(X(:)) max(X(:)) min(Y(:)) max(Y(:)) ])
      caxis([0 tadv(end)-tadv(1)])
      axis equal
      title('Number of days in the domain')
      
      fn =fullfile(dirOut ,['traj_' fnameOut '_float' sprintf('%03d',ifloat)]);
      if fSavePng, saveas(1,[fn '.png']);end
      if fSaveFig, saveas(1,[fn '.fig']);end

      figure(2);clf
      plot3(xadv(ifloat,:), yadv(ifloat,:), zadv(ifloat,:),'.-');grid on
      fn =fullfile(dirOut ,['traj3d_' fnameOut '_float' sprintf('%03d',ifloat)]);
      if fSavePng, saveas(2,[fn '.png']);end
      if fSaveFig, saveas(2,[fn '.fig']);end    
      
  end
end




