function Xpcs=eof_project(Xanom,Xeofs,numd)
% function Xpcs=eof_analysis(Xanom,Xeofs,numd)
% Function to project physical field (1 or 2 spatial 
% dimensions) onto EOF.
% Input:
%   Xanom: (x,y,t) or (x,t) - physical ANOMALY field
%   Xeofs: 3-d (x,y,e) or 2-d (x,e) matrix with spatial 
%          patterns.
%   numd: if Xanom is (x,t), MUST set ndims=1
% Output:
%   Xpcs: 2-d matrix (t,e) with projection coefficients 
%         (principal components) in the columns. Number 
%         of projection coefficients returned is equal 
%         to the number of EOFs passed in the input argument.

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

if nargin<3 % default is 2-d
  numd=2;
end

if ndims(Xanom)==2 & numd==1  % data is (x,t)
  Xpcs=Xeofs'*Xanom;
  Xpcs=Xpcs';
else
%   nt=size(Xanom,3)
%   neofs=size(Xeofs,3)
%   Xpcs=repmat(NaN,[nt neofs]);
%   for it=1:nt
%     Xpcs(it,:)=squeeze(nansum(nansum(Xeofs.*repmat(Xanom(:,:,it),[1 1 neofs]))))';
%   end
%   for ie=1:neofs
%     Xpcs(:,ie)=squeeze(nansum(nansum(Xanom.*repmat(Xeofs(:,:,ie),[1 1 nt]))));
%   end  
% Flatten 2-d fields into single vector
  [nx,ny,nt]=size(Xanom);
  neofs=size(Xeofs,3);  
  Xanom=reshape(Xanom,[nx*ny nt]); % (x,t)
  Xeofs=reshape(Xeofs,[nx*ny neofs]); % (x,e)  
% remove NaNs
  inn=find(~isnan(Xanom(:,1)));  
  Xanom=Xanom(inn,:); % (x,t)
  Xeofs=Xeofs(inn,:); % (x,e);
  Xpcs=Xeofs'*Xanom;
  Xpcs=Xpcs';
end