
load lagrangeterms_world

dLon=5;
dLat=2.5;

%Lon=min(lon_lag):dLon/2:max(lon_lag);
%Lat=min(lat_lag):dLat/2:max(lat_lag);
Lon=-180+dLon/2:dLon/2:180-dLon/2;
Lat=-72:dLat/2:80;

% indices are (Lon, Lat, lag, component)
Kbar=NaN*ones(length(Lat),length(Lon),size(K,2),4);
Upbar=NaN*ones(length(Lat),length(Lon),4); 
    % for convenience in looping; 
    %  Upbar(:,:,3)=Upbar(:,:,4)= <u' v'>
Tbar=Kbar;
NumObs=Kbar;

for i=1:length(Lat);
    fprintf('%d of %d\n',i,length(Lat));
    q=find(lat_lag>=Lat(i)-dLat/2 & lat_lag<=Lat(i)+dLat/2);
    
    for j=1:length(Lon);
        q2=find(lon_lag(q)>=Lon(j)-dLon/2 & lon_lag(q)<=Lon(j)+dLon/2);
        if (length(q2)>5)
            
            kxx=K(q(q2),3,1); qer=find(isfinite(kxx));kxx=kxx(qer);
            
            if (length(kxx)>20e3 & 0)
                plot(sort(kxx),'o'); 
                ttl=['mean: ',num2str(mean(kxx)), ...
                        '     median: ',num2str(median(kxx))];
                title(ttl);
                qer=log10(kxx+1e8); qer=median(qer); 
                qer=10^qer-1e8
                pause;
            end;
            
            Upbar(i,j,1)=nanmean(u_lag(q(q2)).^2);
            Upbar(i,j,2)=nanmean(v_lag(q(q2)).^2);
            Upbar(i,j,3)=nanmean(u_lag(q(q2)).*v_lag(q(q2)));
            Upbar(i,j,4)=Upbar(i,j,3);
            for k=1:size(K,2);  % loop over fixed lag values in T_lag
                for l=1:4;   % tensor components
                    Kbar(i,j,k,l)=nanmean(K(q(q2),k,l));
                    Tbar(i,j,k,l)=Kbar(i,j,k,l)/Upbar(i,j,l)/3600/24;
                    NumObs(i,j,k,l)=length(find(isfinite(K(q(q2),k,l))));
                end;
            end; 
        end;
    end;
end;

kxx=NaN*ones(length(Lat),length(Lon)); kxy=kxx;
kyx=kxx; kyy=kxx;
for i=1:length(Lat);
    for j=1:length(Lon);
        %kxx(i,j)=nanmedian(Kbar(i,j,2:length(T_lag),1)); 
        %kyy(i,j)=nanmedian(Kbar(i,j,2:length(T_lag),2)); 
        %kxy(i,j)=nanmedian(Kbar(i,j,2:length(T_lag),3)); 
        %kyx(i,j)=nanmedian(Kbar(i,j,2:length(T_lag),4));
        q=find(T_lag==20);
        kxx(i,j)=Kbar(i,j,q,1); 
        kyy(i,j)=Kbar(i,j,q,2); 
        kxy(i,j)=Kbar(i,j,q,3); 
        kyx(i,j)=Kbar(i,j,q,4);
    end;
end;

kxyS=(kxy+kyx)/2;
alpha=NaN*ones(size(kxx)); 
q=find(kxx>0 & kyy>0 & (kxx.*kyy)>(kxyS.^2));
alpha(q)=.5* atan2(2*kxyS(q),(kxx(q)-kyy(q)));

% k1: major principal component.  k2: minor principal component
%   Zhurbas and Oh (2003): minor component minimizes impact of shear
%   to yield an appropriate scalar lateral diffusivity.
k1=kxx.*(cos(alpha)).^2 + kxyS.*sin(2*alpha) + kyy.*(sin(alpha)).^2;
k2=kxx.*(sin(alpha)).^2 - kxyS.*sin(2*alpha) + kyy.*(cos(alpha)).^2;
 

uxx=Upbar(:,:,1); uyy=Upbar(:,:,2); uxy=Upbar(:,:,3);
term1=uxx + uyy;
term2=sqrt( (uxx-uyy).^2+4*uxy.*uxy);
u1=sqrt(.5*(term1+term2));
u2=sqrt(.5*(term1-term2));
q=find(term2>term1); u1(q)=NaN; u2(q)=NaN; 
T=k2./u2./u2/3600/24;
%T1=k1./u1./u1/3600/24;
L=k2./u2/1e3;

save meanlagrange Lon Lat NumObs L T Upbar alpha k1 k2 kx* ky* ux* uy*

