
% calculate individual values of effective diffusivity tensor 
%   and eddy speed.  Kxx=K(:,:,1).  Kyy=K(:,:,2).
%   Kxy=K(:,:,3).  Kyx=K(:,:,4).
% These calculations follow the approach of Zhurbas and Oh (2004)
%   JGR 109, C05015


clear all;

filename='world';
eval(['load ',filename,' lon lat t id drogue']);
eval(['load eddyterms_',filename,' up vp'])

if length(id) ~= length(up)
    disp('warning! up length not matching id length')
    keyboard;
end;

q=find(drogue<1);
truncatedrift;
clear drogue;

medt=median(diff(t)); 

T_lag=[5 10 15 20];
K=NaN*ones(length(up),length(T_lag),4);


% loop through trajectory segments of length max(T_lag)/2
ind=2:(max(T_lag)/medt) +2;

while(max(ind)<length(t));
    fprintf('%3.2f percent\n',100*max(ind)/length(t))
    
    if (~diff(id(ind)));
        
        tt=t(ind);
        tt=tt-tt(1);
        uu=up(ind);
        vv=vp(ind);
        x=cumsum(uu).*[diff(tt);last(diff(tt))]*3600*24;
        y=cumsum(vv).*[diff(tt);last(diff(tt))]*3600*24; 
        for j=1:length(T_lag);
            q=find(tt==T_lag(j));
            if (~isempty(q))
                %  k_jk(t)=v'(t_o)d'(-t).
                K(ind(q),j,1)=-uu(q)*(x(1)-x(q));
                K(ind(q),j,2)=-vv(q)*(y(1)-y(q));
                K(ind(q),j,3)=-uu(q)*(y(1)-y(q));
                K(ind(q),j,4)=-vv(q)*(x(1)-x(q));
            end;
        end;
        
    end;
    
    ind=ind+max(T_lag)/medt /2;
    
end;

q=[];
for j=1:length(T_lag);
    q=[q;find(isfinite(K(:,j,1)))];
end;
q=unique(q);
lon_lag=lon(q);
lat_lag=lat(q);
t_lag=t(q);
u_lag=up(q);
v_lag=vp(q);
id_lag=id(q);
K=K(q,:,:);

eval(['save lagrangeterms_',filename,' *_lag K T_lag']);

