function [Xeofs,Xpcs,Var,Xrecon]=eof_analysis(X,neof)
% function [EOFs,PCs,Var,Xrecon]=eof_analysis(X,neof)
% Wrapper function to perform PCA of a field X 
% with TWO spatial dimensions. (This code will also 
% work with a SINGLE spatial dimension. But it might 
% be easier to directly call 'principal_component_analysis'.)
% This function basically transforms the data matrix into 
% a standard 2-d data matrix, taking into account NaNs.
% Input:
%   X: (x,y,t) or (x,t). 
%   neof: number of EOF/PC to return
% Output:
%   EOFs: 3-d (x,y,e) or 2-d (x,e) matrix with spatial 
%         patterns.
%   PCs:  2-d matrix (t,e) with principal components (scores) 
%         in the columns
%   Var:  variance of each principal component
%   Xrecon: 3-d (x,y,t) or 2-d (x,t) matrix with reconstructed X
%           WITHOUT adding back the mean
% Note: Xpc is the projection of X onto the corresponding EOF.
% That is, Xpc(it,ie) = nsum(nsum(X(:,:,it).*Xeofs(:,:,ie)))
% Use this to project a physical field onto EOF space.

% If X only has 2 dimensions, assume second dimension is time.
% Insert 'y' dimension to conform with rest of function which 
% assumes X represents multiple realizations of a 2-D field.

% Samar Khatiwala (spk@ldeo.columbia.edu)

% SPK 3/9/2006: minor changes to use less memory

if ndims(X)==2
  [nx,nt]=size(X);
  X=reshape(X,[nx 1 nt]);
end

% Flatten 2-d fields into single vector
[nx,ny,nt]=size(X);
Xd=reshape(X,[nx*ny nt]); % (space,time)
% remove NaNs and transpose
inn=find(~isnan(Xd(:,1)));
Xd=Xd(inn,:)'; % (time,space)  

% PCA
if nargout>3
  [EOFs,Xpcs,Var,XR]=principal_component_analysis(Xd,neof);
else
  [EOFs,Xpcs,Var]=principal_component_analysis(Xd,neof);
end  
% Xpcs (t,e)

% Reshape EOFs to physical space
Xeofs=repmat(NaN,[nx*ny neof]);
Xeofs(inn,:)=EOFs;
Xeofs=squeeze(reshape(Xeofs,[nx ny neof]));

% Xrecon=repmat(NaN,[nx ny nt]);
% for ie=1:neof
%   for it=1:nt
%     Xrecon(:,:,it)=Xrecon(:,:,it)+Xpcs(it,ie)*Xeofs(:,:,ie);
%   end
% end
  
if nargout>3
  XR=XR'; % (x,t)
  Xrecon=repmat(NaN,[nx*ny nt]);
  Xrecon(inn,:)=XR;
  Xrecon=squeeze(reshape(Xrecon,[nx ny nt]));
end