cd ~dmenemen/llc_2160/regions/latlon/drift
load global_v3/Release.mat

% extract trajectory for species 1 starting
% at 10.291S, 149.346E, on April 1, 2011, hour 13
ix = find( Release.Lon==149.346 & Release.Lat==-10.291 );
disp([Release.Lon(ix) Release.WetLon(ix) Release.Lat(ix) Release.WetLat(ix)])
startDte=datenum(2011,4,1,13,0,0);
dte=startDte:(1/24):(startDte+DriftLength);
Lat=zeros(length(dte),1);
Lon=zeros(length(dte),1);
for t=1:length(dte), mydisp(t)
    iDrift=find(Release.Tme<=dte(t) & Release.Tme>=(dte(t)-DriftLength));
    suf=['_' int2str(length(Release.WetLon)) 'x' int2str(length(iDrift)) '_' datestr(dte(t),30)];
    if exist(['global_v3/DriftLon' suf])
        idx=ix+(t-1)*length(Release.Lat);
        Lon(t)=readbin(['global_v3/DriftLon' suf],1,1,'real*4',idx-1);
        Lat(t)=readbin(['global_v3/DriftLat' suf],1,1,'real*4',idx-1);
    end
end
plot(Lon(find(Lon)),Lat(find(Lat)))

% extract all trajectories starting on March 1, 2012, hour 0
dt=1/24;                                   % model output time step (days)
eps=dt/100;                                % epsilon used for roundoff errors
startDte=datenum(2011,12,1,0,0,0);
dte=startDte:dt:(startDte+DriftLength);
Lat=zeros(length(Release.WetLon),length(dte));
Lon=zeros(length(Release.WetLon),length(dte));
for t=1:length(dte), mydisp(t)
    iDrift=find(Release.Tme<=dte(t)+eps & Release.Tme>=(dte(t)-DriftLength-eps));
    suf=['_' int2str(length(Release.WetLon)) 'x' int2str(length(iDrift)) '_' datestr(dte(t),30)];
    if exist(['global_v3/DriftLon' suf])
        it=find(Release.Tme(iDrift)==startDte);
        Lon(:,t)=readbin(['global_v3/DriftLon' suf],length(Release.WetLon),1,'real*4',it-1);
        Lat(:,t)=readbin(['global_v3/DriftLat' suf],length(Release.WetLon),1,'real*4',it-1);
    else
        error('choubichou')
    end
end
clf
plot(Lon',Lat','.','markersize',1)
plotland('k',12,-38:(1/2):322,-50:(1/12):50)
axis('equal')
axis([-38 322 -50 50])
eval(['print -dpsc figs/YearDrift' datestr(startDte,30)])
eval(['print -djpeg -r600 figs/YearDrift' datestr(startDte,30)])