% calc_transit_matrix
% 
% calculate the transit matrix for float locations
% dT days later.   Do not normalize: i.e., there will be fewer
% floats after dT days due to them dying from various causes.
% This can be normalized later, if the tracer should be
% conserved.  
%
% This version ignores seasonality.

 
clear all;

if (0)  % load yomaha file, make evenly interpolated version
    worldfile=matfile('../yomaha_deep.mat');
    lon=worldfile.lon;
    lat=worldfile.lat;
    t=worldfile.t;
    id=worldfile.id;
    q=find(lon<0); lon(q)=lon(q)+360;
    clear q ind worldfile
    % loop through, interpolating to daily intervals
    % (this is to ensure that all bins are occupied by trajectories that
    %  otherwise can jump across several bins in a 10 day interval)
    ID=unique(id);
    lon2=[]; lat2=[]; id2=[]; t2=[];
    for i=1:length(ID);
        fprintf('%d of %d\n',i,length(ID));
        q=find(id==ID(i));
        if (q>1)  % drop if there is only one point
            if (median(diff(t(q)))<30)  % only retain and interpolate regularly-sampling floats
                t22=(min(t(q)):1:max(t(q)))';
                lon2=[lon2;rewraplon(interp1(t(q),unwraplon_nodespiker(lon(q)),t22,'linear'))];
                lat2=[lat2;interp1(t(q),lat(q),t22,'linear')];
                id2=[id2;id(q(1))*ones(length(t22),1)];
                t2=[t2;t22];
            end;
        end;
    end; 
    clear t22 ID i
    lon=lon2; lat=lat2; id=id2; t=t2; clear *2
    q=find(lon<0); lon(q)=lon(q)+360;
    clear q 
    
    lon2=3/2:3:360-3/2;
    lat2=floor(min(lat))+3/2:3:ceil(max(lat))-3/2;
    NumObs=NaN*ones(length(lat2),length(lon2));
    for i=1:length(lat2);
        q=find(lat>=lat2(i)-3/2 & lat<lat2(i)+3/2);
        for j=1:length(lon2);
            q2=find(lon(q)>=lon2(j)-3/2 & lon(q)<lon2(j)+3/2);
            NumObs(i,j)=length(q2);
        end;
    end;
    
    save('eraseme.mat','-v7.3');
else
    worldfile=matfile('eraseme.mat');
end;


dT=30;  % number of days later
 

load /data/lumpkin/surfdrift/tracer/OceanTracerTransitMatrix lon lat
q=find(lon<0); lon(q)=lon(q)+360;

% subset for the Pacific (TPOS exercise)
%q=find(lon<100 | lon>290); lon(q)=[];
%q=find(lat<-30 | lat>30); lat(q)=[];

dLon=abs(median(diff(lon)));
dLat=abs(median(diff(lat)));
latmat=ones(length(lon),1)*lat;
lonmat=lon'*ones(1,length(lat));
latmat=latmat(isfinite(latmat));
lonmat=lonmat(isfinite(lonmat));
  
p=sparse(length(latmat),length(latmat)); 

loadflag=1; 
worldfile=matfile('eraseme.mat');
latB=10;  % load data in band +-latB
latBflag=0;
 
for i=1:length(latmat);  
    
    if (i/100==floor(i/100))
        fprintf('%d of %d\n',i/100,round(length(latmat)/100));
    end;
    
    if (loadflag)  % load the data, truncate 
        lon=worldfile.lon;
        lat=worldfile.lat;
        t=worldfile.t;
        id=worldfile.id; 
        % by construction, looping through indices of lonmat, latmat
        % will loop slowly through latitudes, quickly through longitudes
        %   i.e., latmat(1:360)=-74
        ind=find(lat>latmat(i)-latB & lat<latmat(i)+latB);
        lon=lon(ind); lat=lat(ind); t=t(ind); id=id(ind);  
        loadflag=0;
    end;
     
    q=find(lat>=latmat(i)-dLat/2 & lat<latmat(i)+dLat/2);
    q=q(lon(q)>=lonmat(i)-dLon/2 & lon(q)<lonmat(i)+dLon/2);  
    if (~isempty(q))
        % find all IDs going through the box
        ID=unique(id(q));
        % loop through IDs
        for j=1:length(ID);
            q2=find(id==ID(j));
            % find closest approach to center of box
            del=(lon(q2)-lonmat(i)).^2+(lat(q2)-latmat(i)).^2;
            t2=t(q2);
            % define t2=0 at closest approach
            mini=find(del==min(del));
            mini=max(mini);
            t2=t2-t2(mini);
            
            % look for data at t2=dT-5 to dT+5 
            qq=find(t2>=dT-5 & t2<=dT+5); 
            if (~isempty(qq));    % this won't be true if float dies
                                % before dT days after closest to
                                % lonmat(i),latmat(i), or has gap
                %qp=q2(mini:qq(1));
                %plot(lon(qp),lat(qp),'k','linewidth',.5);
                %plot(lon(qp(length(qp))),lat(qp(length(qp))),'ok','markersize',4,'markerfacecolor','w');
                ind=sqrt((latmat-lat(q2(qq(1)))).^2+(lonmat-lon(q2(qq(1)))).^2); 
                c0=ind(find(ind<1)); c0=1-c0; c0=c0/sum(c0)/length(ID);
                ind=find(ind<1);
                p(i,ind)=p(i,ind)+c0';
                
                if (lat(q2(qq(1)))>=max(lat)-latB/2 && latmat(i)<80)
                    loadflag=1;  
                    if (lat(q2(qq(1)))>=max(lat)-1)
                        disp('WARNING: hitting edge of data') 
                        latBflag=1;
                    end;
                end;
                
            end;
            %  p2=reshape(full(p(i,:))',[length(unique(lonmat)) length(unique(latmat))])';
            

        end; % loop over IDs
        if (latBflag)
            latB=latB+1
            latBflag=0;
        end;
   
    end; % if length(q)  ( if floats are in box)
     
    
end;  % loop over lonmat, latmat 
 

latest_data=datestr(max(t)+722815);
created_on=datestr(now);
 
save transit_matrix lonmat dT latmat p latest_data created_on


 