% Find capture indices and times
cd ~dmenemen/llc_2160/regions/latlon/drift
load global_v3/Release.mat
N=length(Release.Lon);
pin='~dmenemen/llc_2160/regions/latlon/';  % model output files location
eval(['load ' pin 'grid/coastline'])       % load model coastline indices
iCoast=find(coastline);                    % coastline location indices (relative to XC/YC)
nx=8640; ny=4320;                          % model grid dimensions
suf=['_' int2str(nx) 'x' int2str(ny)];
for fld={'XG','YG'}                        % load model grid info
  fnm=[pin 'grid/' fld{1} suf];
  eval([fld{1} '=readbin(fnm,[nx ny]);'])
end
xg=XG(:,500);                              % longitude of grid corners
yg=YG(6000,:);                             % latitude of grid corners
xg(find(xg<xg(1)))=xg(find(xg<xg(1)))+360;
xg=[xg; 360];
clear *G c*

% create one file per release time that document all capture events
CaptureIdx{N}=[];                          % Capture index (relative to XC/YC)
CaptureTme{N}=[];                          % Capture time (decimal days from release time)
for r=1:length(Release.Tme)
    disp(datestr(Release.Tme(r)))
    flat=['global/DriftLat_n4331_dt3600_' datestr(Release.Tme(r),30)];
    flon=['global/DriftLon_n4331_dt3600_' datestr(Release.Tme(r),30)];
    D=dir(flon);
    Lon=readbin(flon,[N D.bytes/4/N]);     % propagule trajectories longitude
    Lat=readbin(flat,[N D.bytes/4/N]);     % propagule trajectories latitude
    tme=Release.Tme(r):(1/24):(Release.Tme(r)+DriftLength);
    [n idx]=histc(Lon(:),xg);
    [n idy]=histc(Lat(:),yg);
    idx=reshape(idx,size(Lon));            % longitude index (relative to xg)
    idy=reshape(idy,size(Lat));            % latitude index  (relative to yg)
    dx=diff(idx,1,2);
    dy=diff(idy,1,2);
    for i=1:N
        it=find(dx(i,:)|dy(i,:));
        IT=sub2ind([nx ny],idx(i,it+1),idy(i,it+1));
        LIA=ismember(IT,iCoast);
        CaptureIdx{i}=uint32(IT(find(LIA)));
        CaptureTme{i}=single(tme(it(find(LIA))+1)-Release.Tme(r));
    end
    eval(['save global/Capture' myint2str(r,4) ' Capture*'])
end

% create and save one array that only contains first-capture events
FirstCaptureIdx{N,length(Release.Tme)}=[]; % Capture index (relative to XC/YC)
FirstCaptureTme{N,length(Release.Tme)}=[]; % Capture time (decimal days from release time)   
for r=1:length(Release.Tme), mydisp(r)
    eval(['load global/Capture' myint2str(r,4)])
    for i=1:N
        Idx=CaptureIdx{i};
        Tme=CaptureTme{i};
        n=1;
        while n<length(Idx)
            in=find(Idx==Idx(n));
            Idx(in(2:end))=[];
            Tme(in(2:end))=[];
            n=n+1;
        end
        FirstCaptureIdx{i,r}=Idx;
        FirstCaptureTme{i,r}=Tme;
    end
end
save -v7.3 global/FirstCapture FirstCapture*
save -v7.3 global/FirstCaptureIdx FirstCaptureIdx
save -v7.3 global/FirstCaptureTme FirstCaptureTme

% save by species
ReleaseTme=Release.Tme;
for i=1:length(Release.SpeciesLocation), mydisp(i)
    SpeciesLocation=Release.SpeciesLocation{i};
    ReleaseLat=Release.Lat(SpeciesLocation);
    ReleaseLon=Release.Lon(SpeciesLocation);
    ReleaseWetLat=Release.WetLat(SpeciesLocation);
    ReleaseWetLon=Release.WetLon(SpeciesLocation);
    clear PotStrandIdx PotStrandTme
    PotStrandIdx{length(SpeciesLocation),N}=[];
    PotStrandTme{length(SpeciesLocation),N}=[];
    for s=1:length(SpeciesLocation)
        for t=1:length(Release.Tme)
            PotStrandIdx{s,t}=FirstCaptureIdx{s,t};
            PotStrandTme{s,t}=FirstCaptureTme{s,t};
        end
    end
    eval(['save global/Species' myint2str(i) ...
          ' Dr* Po* ReleaseL* ReleaseT* ReleaseW* Sp*'])
end