%function [pcs,eofs,vxpl,conf,lambda,weights]=eof(data,area);
function [pcs,eofs,vxpl,lambda,weights]=eof(data,area);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
% Computes the EOFs from the covariance matrix by SVD
%
% Input:  data    - data matrix (time,space)                    
%
% Output: pcs  	- principal components (columns)
%         eofs    - empirical orthogonal functions (columns)
%         lambda  - the eigenvalues of the covariance/correlation matrix   
%         vxpl    - variance represented by the EOFs
%
% Example: [pcs,eofs,v,l,w]=eof(mydata,myarea);
%
%
% Ingo Bethke 2004 (based on is_sdv by Ivar Seierstad,2003) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

warning off MATLAB:divideByZero

% Checking input arguments
if nargin<1
help eof; return; end

% Remove mean from time series
data=detrend(data,'constant');
disp('Data are detrended');

% Remove zero-variance points
[d1,d2old]=size(data);
subind=find(sum(abs(data),1)~=0);
data=data(:,subind);
[d1,d2]=size(data);
area=area(subind);

% Scale area 
area=area/max(area(:));

% Keep a copy of the data matrix prior to normalization/weighting
dataold=data;  

% Area weighting 
weight=sqrt(area);
data=data.*meshgrid(weight,1:d1);
disp('Data are weighted by the sqrt of the cell area');   
weights=weight; % return the weights
     
% EOFs calculated by SVD
if d1>d2
  [U,L,A]=svd(data,0);
  lambda=L.^2;         % the singular values are the square roots of the eigenvalues
  U=U*L;               % amplitude information in PCs      
  A=A';                % put the EOFs as columns 
else           
  [A,L,U]=svd(data',0);
  lambda=L.^2;         % the singular values are the square roots of the eigenvalues
  U=L*U;               % amplitude information in the PCs 
  U=U';                % put the PCs as columns
end;
disp('The eigenvectors are calculated by singular value decomposition')
vxpl=100*abs(diag(lambda)/sum(diag(lambda))); % Variance represented by the PCs


% Scale PCs
disp('PCs have st.dev 1');
pcs=U./sqrt(meshgrid(diag(lambda)/(d1-1),1:size(U,1))); % PCs normalized


% Data regressed on principal components, so as to represent
% the amplitude of the data
eofssub=((pcs*dataold)/(d1-1))'; 
disp('EOFs scaled to represent the amplitude of the data');  

% Fill gaps with zeros
%eofs=nans(d2old,d1);eofs(subind,:)=eofssub;
eofs=nan(d2old,d1);eofs(subind,:)=eofssub;  % cguo

% Estimate confidence for the explained variance  
disp('Estimate confidence for explained variance');
%c=autocorr(pcs');
%for s=1:size(c,2)
% ind=find(c(:,s)<0.3679);
% T(s)=mean([ind(1) max([1 (ind(1)-1)])]);
%end
%T=mean(T);
%N=size(c,1);
%conf=vxpl*2*sqrt(T/N);
