%initialize
pin='~dmenemen/llc_2160/regions/latlon/'; % model output files location
eval(['cd ' pin 'drift'])
N=4331;                                   % total number of release locations
nx=8640; ny=4320;                         % model grid dimensions
suf=['_' int2str(nx) 'x' int2str(ny)];
for fld={'XC','YC','hFacC'}               % load model grid info arrays
  fnm=[pin 'grid/' fld{1} suf];
  eval([fld{1} '=readbin(fnm,[nx ny]);'])
end
XC(find(XC<XC(1)))=XC(find(XC<XC(1)))+360;% unwrap longitude array
xc=XC(:,500);                             % longitude of grid centers
yc=YC(6000,:);                            % latitude of grid corners
land=abs(hFacC-1);                        % array of ones for land and zeros for water

% load species information file
s=1;                                      % species index (1:94)
eval(['load global/Species' myint2str(s)])

% load trajectories
r=1;                                      % release time index (1:8785)
flat=['global/DriftLat_n4331_dt3600_' datestr(ReleaseTme(r),30)];
flon=['global/DriftLon_n4331_dt3600_' datestr(ReleaseTme(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
Lon=Lon(SpeciesLocation,:);
Lat=Lat(SpeciesLocation,:);

% plot trajectory and potential stranding events
t=1;                                      % trajectory index (1:length(SpeciesLocation))
figure(1), clf
ix=find(xc>=min(Lon(t,:)-.2)&xc<=max(Lon(t,:)+.2));
iy=find(yc>=min(Lat(t,:)-.2)&yc<=max(Lat(t,:)+.2));
shadeland(xc(ix),yc(iy),land(ix,iy)',.8*[1 1 1]);
hold on
plot(Lon(t,:),Lat(t,:))
plot(Lon(t,1),Lat(t,1),'k*','linewidth',2,'markersize',16)
plot(Lon(t,end),Lat(t,end),'r*','linewidth',2,'markersize',16)
plot(XC(PotStrandIdx{t}),YC(PotStrandIdx{t}),'go')
title(['Species ' int2str(s) ', trajectory ' int2str(t) ...
       ', released on ' datestr(ReleaseTme(r))])

% plot potential stranding events distance and time
figure(2), clf
Distance=dist2([ReleaseLat(t) YC(PotStrandIdx{t})], ...
               [ReleaseLon(t) XC(PotStrandIdx{t})]);
plot(PotStrandTme{t},Distance,'.')
xlabel('Time from release (days)')
ylabel('Distance from release (km)')
title(['Potential stranding events for species ' int2str(s) ...
      ', trajectory ' int2str(t) ', released on ' datestr(ReleaseTme(r))])