 function [deof,pev,dev0,dev] = deof_2sets_analysis(eof0, eof, ev0, ev, lat)
%
% evaluating EOFs of two data sets against each other to find the 
% DEOF pattern that are more important in terms of explained variance
% in data set2 than in data set1. See Bayr and Dommenget (2012) for details.
% 
% Input:
%
% eof0(neof,xdim,ydim): eofs of data set1 (this eofs will not be rotated)
%  eof(neof,xdim,ydim): eofs of data set2 (this eofs will be rotated)
%           (undefined values should be set to zero in both eof sets)
%                  ev0: explained variances of eof0
%                   ev: explained variances of eof
%                  lat: latitudes as vector for lat weighting sqrt(cos(lat))
%                        
%
% Output:   
%
% deof:  DEOF pattern of the data-matrix that have a higher 
%         explained variance in data set2 than in data set1     
%  pev:  projection of eigenvalues of data set1 on data set2 
%        (compare with ev0; eq.[1] and [2] of Bayr and Dommenget (2012)) 
% dev0:  exp. variances of DEOFs in data set1        
%  dev:  exp. variances of DEOFs in data set2
%
% Author: Tobias Bayr (tbayr@geomar.de) based on a Dietmar Dommenget script
% date: 18.09.2012 V1.0

% parameters & initialization
[neof,xdim,ydim] = size(eof);
xydim = xdim*ydim;           
deof  = zeros(neof,xdim,ydim);

% latitudes weight (apply same lat weight in eof analysis
% to ensure that pattern are orthogonal!)
wg = ones(xdim,1) * sqrt(cosd(lat))';
wgv = reshape(wg,xydim,1);

% reshape data(neof,xdim,ydim) -> xydata(xydim,neof)
eofv = reshape(eof,neof,xydim)';
eof0v = reshape(eof0,neof,xydim)';

% latitude weight data
eofv = eofv .* (wgv * ones(1,neof));
eof0v = eof0v .* (wgv * ones(1,neof));

% normalize eofs
for k=1:neof
    eofv(:,k) = eofv(:,k)./sqrt(eofv(:,k)'*eofv(:,k));
    eof0v(:,k) = eof0v(:,k)./sqrt(eof0v(:,k)'*eof0v(:,k));
end

% warning if eofs are not orthogonal
pcor1 = eofv(:,1)' * eofv(:,2); 
pcor2 = eof0v(:,1)' * eof0v(:,2);
if abs(pcor1) > 1e-12 || abs(pcor2) > 1e-12
    warning('Warn:Unorthogonal',...
'EOF pattern are not orthogonal to each other. Please check lat weighting!');
end

% project null-eof on eof (eq. [1] and [2] Bayr and Dommenget 2012)
r = (eofv'*eof0v); pev = ev0*(r.*r)';  

% find DEOF (rotation of eof around eof0)
[deofv, dev, dev0]=deof_rotation(eofv,ev,eof0v,ev0);

% reshape & scale fields
for n=1:neof
    deof(n,:,:)=sqrt(sum(ev0)*dev(n)/100)*reshape(deofv(:,n)./wgv,xdim,ydim);
end

return

% ====================================================================
function [reofv, rev, rev0] = deof_rotation(eofv,ev,eof0v,ev0)
% estimate maximum variance diff. by pairwise rotation 
% ====================================================================
neof = length(eofv(1,:)); rev=ev; reofv=eofv; 
da   = 0.05;                 % step width in rotation    [fractions of pi]
amax = pi;                   % maximum angle of rotation
a=0.:da*pi:amax; 
iter = 3;                    % number of interations of rotation
xdif = zeros(1,neof);

% rotation
for i0=1:neof
    for k=1:iter                                                        
        for i2=i0+1:neof
            x1=cos(a)'*reofv(:,i0)'+sin(a)'*reofv(:,i2)';
            x2=sin(a)'*reofv(:,i0)'-cos(a)'*reofv(:,i2)';
            xev1=ev*((eofv'*x1').^2);  xev2=ev*((eofv'*x2').^2);
            evbx=ev0*((eof0v'*x1').^2);
            [xd ii]=max(xev1-evbx);         % eq. [10] of (Dommenget 2007)
            reofv(:,i0)=x1(ii,:); rev(i0)=xev1(ii); 
            reofv(:,i2)=x2(ii,:); rev(i2)=xev2(ii);
        end
        xdif(i0)=xd;
    end
end
rev0=rev-xdif;
return
