function [eof, eof0, deof, ev, pev0, dev, dev0] = deof_analysis(data,dlon,dlat,lat0,ndgf)
%
% evaluating EOFs against a fitted isotropic diffusion process as discussed
% in section 4 of Dommenget (2006)
% 
% Input: data(tdim,xdim,ydim): data matrix of anomalies
%                              (undefined values should be set to zero)
%                   dlon,dlat: step width in logitudes,latidues
%                        lat0: latidues of first point
%                        ndgf: estimated number degree of freedom of the time series of data  
%
% Output:   eof = EOFs of the data-matrix         (E_obs)
%            ev = Eigenvalues of the data-matrix  (e_obs) 
%          eof0 = EOFs of fitted process          (E_null) 
%          pev0 = eq.[9] of Dommenget (2006)      (e_obsnull) 
%          deof = DEOFs of the data-matrix        (D_obs)
%           dev = exp. variances of DEOFs          
%          dev0 = exp. variances of DEOFs under the null hypothesis          
%
% Author: Dietmar Dommenget (dietmar.dommenget(a)monash.edu)
%

% parameters & initialization
neof =20;                                     % number of EOFS computed
xdim   =length(data(1,:,1));   ydim =length(data(1,1,:));
tdim   =length(data(:,1,1));  xydim =xdim*ydim;           
xydata =zeros(tdim,xydim);       wg =zeros(xydim,1);    
eof    =zeros(neof,xdim,ydim);   ev =zeros(1,neof);
eof0   =zeros(neof,xdim,ydim);  ev0 =zeros(1,neof);
deof   =zeros(neof,xdim,ydim);  icv =zeros(xydim,xydim);
jcv    =zeros(xydim,xydim); 

% latitudes weight
for i=1:ydim; wg(xdim*(i-1)+(1:xdim))=cos(6.28*(dlat*(i-1)+lat0)/360); end 

% reshape data(tdim,xdim,ydim) -> xydata(tdim,xydim)
for t=1:tdim;  xydata(t,:)=wg'.*reshape(data(t,:,:),1,xdim*ydim); end

% standard deviation field
sig=squeeze(std(xydata,1));

% eofs
if ( tdim >= xydim ); cv=(xydata'*xydata)/tdim; end
if ( tdim <  xydim ); cv=(xydata*xydata')/xydim; end
[U,S] = svds(cv,neof);
for k=1:neof; ev(k)=100*S(k,k)/sum(trace(cv)); end
if ( tdim >=  xydim ); eofv=U; end
if ( tdim <  xydim )
    for k=1:neof
        eofv(:,k)=(xydata'*U(:,k))'./sqrt((xydata'*U(:,k))'*(xydata'*U(:,k)) );
    end   
end
neff=1/sum((ev/100.).^2);             % effective spatial number degrees of freedom (N_eff)

% create the null-hypo (isotropic diffusion process)
joff=lat0/dlat;
ix=1:xdim; for j=1:ydim; jcv(xdim*(j-1)+ix,:) =j; end
jx=1:ydim; for i=1:xdim; icv(xdim*(jx-1)+i,:) =i; end
d=dsphere(icv,jcv+joff,icv',jcv'+joff,dlon,dlat); 
d0=sqrt((xdim*ydim*dlat*dlon))/neff;  % first guess decorrelation length [degree]
dsmin=0.5/sqrt(ndgf);                 % min of diff. in neff eof vs eof0 [fraction of sdgf]
ds=1.0;
while abs(ds) > dsmin 
    cv0=(sig'*sig).*exp(-d/d0);
    [eof0v,S0] = svds(cv0,neof);      % eofs of null hypo
    for k=1:neof; ev0(k)=100*S0(k,k)/sum(trace(cv0)); end
    neff0=1/sum((ev0/100.).^2); ds=neff0/neff-1;  d0=(ds+1)*d0;
end

% project null-eof on eof (eq. [8] and [9] Dommenget 2006)
r=(eofv'*eof0v); pev0=ev0*(r.*r)';

% find DEOF
[deofv, dev, dev0]=deof_rotation(eof0v,ev0,eofv,ev);

% reshape & scale fields
for n=1:neof
    eof(n,:,:) =sqrt( S(n,n))*reshape( eofv(:,n)./wg,xdim,ydim);
    eof0(n,:,:)=sqrt(S0(n,n))*reshape(eof0v(:,n)./wg,xdim,ydim);
    deof(n,:,:)=sqrt(sum(trace(cv))*dev(n)/100)*reshape(deofv(:,n)./wg,xdim,ydim);
end
return

% ====================================================================
function [reofv, rev, rev0] = deof_rotation(eof0v,ev0,eofv,ev)
% 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 2006)
            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

% ====================================================================
function d = dsphere(i,j,i0,j0,dlon,dlat)
% distance of two points on a sphere
% ====================================================================
dx=cos(pi/180*dlon*i).*cos(pi/180*dlat*j)-cos(pi/180*dlon*i0).*cos(pi/180*dlat*j0);
dy=sin(pi/180*dlon*i).*cos(pi/180*dlat*j)-sin(pi/180*dlon*i0).*cos(pi/180*dlat*j0);
dz=sin(pi/180*dlat*j)-sin(pi/180*dlat*j0);
dxyz=sqrt(dx.^2+dy.^2+dz.^2);
d=360/pi*asin(dxyz/2);
return




