function [EOFm,PC,spec,RC]=EOF_hepta(Z,lat,thre);
% function [EOFm,PC,spec,RC]=EOF_hepta(Z,lat,thre)
%   
%   EOF analysis and synthesis of matrix Z, 
%   
%   feat.  @ area-weighting (sqrt(cos(lat))) of data matrix as in North et al 1982
%           @ PCs with unit variance . 
%           @ EOF  patterns scaled by the eigenvalue squared
%               ...  i.e. the total (area-weighted) variance they explain  . 
%            
%   in other words :  
%   EOF patterns have a physical meaning, the PCs have unit energy .
%      (this is somewhat unusual so make sure this is what you want)
%  
%     INPUTS :  - Z (nt,ny,nx) 
%      Z is assumed to contain NO NAN's . Get back to the drawing board if it does.   
%                    - lat : latitude vector  (length ny) in degrees
%                    - THRE : variance threshold for PCA reconstruction ( in %)   
%                                [default is 90%]   
%                     
%     OUTPUTS :
%                    - EOFm matrix, size nm*ny*nx (nm= number of modes retained)
%                     - PC matrix (nt*nm)
%                     - spec : (full) eigenvalue spectrum in terms of % variance explained
%                     - RC : PCA-reconstructed data matrix with nm modes 
%        
%    note : nm is such that sum(spec(1:nm)) >= THRE ;
%   
%   ref :
%   North, Gerald R., Bell, Thomas L., Cahalan, Robert F., Moeng, Fanthune J.
% "Sampling Errors in the Estimation of Empirical Orthogonal Functions"
%  Monthly Weather Review 1982 110: 699-706
%      
%   thanks to Yochanan Kushnir (LDEO) for precious advice.    
%   
%   Hepta Technologies, Aug 2006, Columbia University, New York.  
%   v1.1 : Bug in PC normalization fixed, Feb 2007, GaTech.
%   v1.2 : Now only retains nm modes (smaller matrices)     July, 25 2007, GaTech 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
if sum(isnan(Z)) >0
   disp('I told you not to put any NaNs in there, damnit !!')
end

%
if nargin<3
   thre=90;
end


% get sizes
[nt,ny,nx]=size(Z);
np=nx*ny; q=min(nt,np);
%%%%%% scale by area 
Zr=zeros(size(Z));
for jj=1:ny
   Zr(:,jj,:)=Z(:,jj,:)*sqrt(cos(pi/180*lat(jj)));
end

% reshape 
Zw=reshape(Zr,nt,np);
%  make sure avg is zero
for jp=1:np
   Zw(:,jp)=Zw(:,jp)-mean(Zw(:,jp))*ones(nt,1);
end


%%%%%%%%%%%% this is the intensive part %%%%%%%%%
% Determine singular value decomposition of the problem
[E,S,P]=svd(Zw',0); % works fast iff the SECOND dimension is small

% rewrite the large eigenvalue matrix to a vector and 
% apply the appropriate normalization
eigv=S.^2/(q-1);    % has q elements
%
spec=100*diag(eigv)/sum(diag(eigv)); % percentage of variance explained

% number of modes to retain
crit=abs(cumsum(spec)-thre);
r=find(crit==min(crit));  

% Normalize the EOFs
Et=E*sqrt(eigv(:,1:q)); % HERE LIES THE VARIANCE INFO, not in PC
EOFm=reshape(Et,ny,nx,q); 
EOFm=EOFm(:,:,1:r);

% Determine the EOF coefficients (Principal Components)
PC=P(:,1:r);
%
% Reconstruction using variance threshold

RC=E(:,1:r)*S(1:r,1:r)*P(:,1:r)';



% sanity check   
vz=sum(var(Zw,1));
vrc=sum(var(RC,1));
 %
RCt=reshape(RC',nt,ny,nx); 
RC=zeros(size(RCt));
%rescale by inverse of area, dim nt*ny*nx
for jj=1:ny; % north and south poles are ditched : bastard 
   if (cos(pi/180*lat(jj))==0) 
      RC(:,jj,:)=0;
   else
      RC(:,jj,:)=RCt(:,jj,:)/sqrt(cos(pi/180*lat(jj)));
   end
end

return