function [evector,ew,amps,fac,covmat]=geteof(data,col,filtpar,area);
% function [evector,ew,amps,fac,covmat]=geteof(data,[col],[filtpar],[area]);
%
% calculates the eigen-vectors and -values of the covariance matrix of 'data'
% the EOFs 'ev' are the solution of (data*data')/(length(dim2)-1) * ev = ev * ew
% the 'amps' solve  data = ( amps' * ev' )'
%
% NaN's are possible but the amplitudes are not any longer orthogonal!
%
% input  : data   		- the data set (no data : NaN)
%				  assumed to have 1 temporal and 1 or more
%				  spatial dimensions
%	   col		[1]	- column of the temporal dimension
%	   filtpar	[0]	- -1 for no tranposing
%				  -2 for no normalization
%				  >0 for filtering covariance matrix
%				  filtpar is the 3rd argument for filt2d
%				  DON'T DO THIS !!!
%
% the next is experimental !!!
%	   area	        [1]	- one weight for EACH grid point (in time
%                                 and space !)
%
% output : evector     		- 'vertical' (temporal) eigen-vectors of 'data'
%				  normalized to 1
%          ew     		- corresponding eigen-values
%				  ew/sum(ew) gives explained variance
%	   amps   		- (spatial) amplitudes of evector
%				  sorted in columns (as 'evector')
%          fac    		- number of realisations of pairs
%	   covmat 		- covariance matrix
%
% uses :	nans.m  filt2d.m
%
% version 1.3.1		last change 30.07.2002

% Gerd Krahmann, IfM Kiel, Feb 1993
%
% Ref.: U.Send's  EOF notes
% and 
% http://www.dkrz.de/forschung/reports/report11-1.1/pingo-45.html
% for the area stuff

% added filtering mode 		G.Krahmann, Sep 1994
% changed filtering mode	G.Krahmann, Jan 1995
% changed output order		G.Krahmann, Sep 1995
% transponing for speed		G.Krahmann, Feb 1996
% removed problem with amps	
% and introduced 'gap_par'	G.Krahmann, Mar 1996
% added compatibility to MATLAB 5	G.Krahmann, LODYC Paris, Aug 1997	
% changed multi-dim handling		G.Krahmann, LODYC Paris, Aug 1997
% added area handling		
% and removed comp. to MATLAB 4 G.K., LDEO Jul 2000
% bugfix			G.K., LDEO Jul 2002	1.3.0-->1.3.1

% add weights, if demanded
if nargin>3
  disp('ask Gerd, if this fails')
  data = data.*sqrt(area);
end

% check for MATLAB version and reshape/reorder if necessary
nd=ndims(data);
if nargin>1
  if isempty(col)
    col=1;
  end
else
  col=1;
end
dimorder=[1:nd];
dimorder(col)=nan;
dimorder=dimorder(find(~isnan(dimorder)));
dimorder=[col,dimorder];
if nd>2
  data=permute(data,dimorder);  
  oldsize=size(data);
  data=data(:,:);
end


% transpose the data matrix if s(1)>s(2)
% to get increased speed
sd=size(data);
if sd(1)>sd(2)
  data=data';
  transp=1;
  if nargin>2
    if filtpar==-1
      data=data';
      transp=0;
    else
      disp('Matrix transposed for increased speed!')
    end
  else
    disp('Matrix transposed for increased speed!')
  end
else
  transp=0;
end


% check if NaN's exist in data
% if yes, calculate the number of data-pairs for calculating the
% covariance matrix separately for each pair, else use the
% length of the series as number of pairs
h1=~isnan(data);
if (length(data(:))~=sum(h1(:)))
    nan_exist=1;
    s=size(data);
    fac=zeros(s(1));
    for i=1:s(1)
      for k=1:s(1)
        fac(i,k)=length( find( h1(i,:).*h1(k,:) ) ) ;
      end
    end
    allnan=find(sum(h1)==0);
    fac = fac-1;
    if ~isempty(find(fac<=0))
      disp('ask Gerd')
    end
else
    nan_exist=0;
    s=size(data);
    fac=s(2);
    allnan=[];
end


% replace NaNs with 0, that gives no contribution when calculating the
% covariance matrix 
[data,nan_ind]=nans(data,0,nan);


% calculate covariance matrix
% and divide by number of realizations
covmat=data*data';
covmat=covmat./fac;

h=find(isnan(covmat));
if length(h) > 0
   covmat(h)=zeros(1,length(h));
end


% filter covariance matrix, if wanted
if (nargin>2)
  if filtpar<1 & filtpar~=-1 & filtpar~=-2 & filtpar>0
    covmat=filt2d(covmat,2,filtpar);
  end
end


% calculate Singular Value Decomposition
[u,s,v]=svd(covmat);
evector=u;
ew=s;
ew=diag(ew);


% normalize Eigenvectors to rms=1
if transp==0
  norms=0;
  if exist('filtpar')
    if filtpar==-2
      norms=1;
    end
  end
  if norms==0
    rev=rms(evector);
    evector=evector./(ones(size(evector,1),1)*rev);
  end
end


% project to get amplitudes
amps=data'/evector';


% if the data matrix was transposed, reorder output arguments
% and correct for rms(evector)=1
% at the same time check for time series full of NaN's
if transp==1
  ew=ew/sd(2)*sd(1);
  evn=amps(:,1:length(ew));
  ampsn=evector(:,1:length(ew));
  evector=evn;
  amps=ampsn;
  if ~isempty(allnan)
    evector(allnan,:)=evector(allnan,:)*nan;
  end
  rev=rms(evector);
  evector=evector./(ones(size(evector,1),1)*rev);
  amps=amps.*(ones(size(amps,1),1)*rev);
else
  if ~isempty(allnan)
    amps(allnan,:)=amps(allnan,:)*nan;
  end
  ew=ew(1:sd(1));
end


% check for MATLAB version and reshape if necessary
if exist('oldsize')
  if length(oldsize)>2
    if transp==0
      amps=reshape(amps',[size(amps,2),oldsize(2:length(oldsize))]);
    else
      amps=reshape(amps',[size(amps,1),oldsize(2:length(oldsize))]);
    end
    [dummy,ind]=sort(dimorder);
    amps=permute(amps,ind);
  end
end

% add weights, if demanded
if nargin>3
  disp('ask Gerd, if this fails')
  if ndims(amps)==2
    amps = amps./sqrt(area');
  else
    amps = amps./sqrt(area);
  end
end

