function [ev,ew,amps,rec]=getmssa(data,llag,ldim);
% function [ev,ew,amps,rec]=getmssa(data,llag,[ldim]);
%
% calculate (multichannel) singular spectrum analysis MSSA
% for SSA simply give a vector as data-array
% the method is similar to extended EOFs (EEOF)
%
% NaN's are possible and internally replaced by random values
% they should not exceed 10% of the series
%
% input  : data   		- the data-array (max 2-dim)
%	   llag			- length of lag (should not exceed one third
%				  of the length of the series)
%	   [ldim]	[1]	- dimension to be lagged
%
% output : ev     		- 'vertical' eigen-vectors of lagged 'data'
%				  normalized to 1
%				  MATLAB 5: dim1=lag dim2=series dim3=EOF#
%				  MATLAB 4: dim1=lag+series dim3=EOF#
%          ew     		- corresponding eigen-values
%				  ew/sum(ew) gives explained variance
%	   amps   		- amplitudes of ev
%				  sorted in columns (like 'ev')
%	   rec			- reconstructed fields (ev*amps)
%				  only available for MATLAB 5 or higher
%				  Warning: this is very slow !
%
% version 0.1.1		last change 15.09.1997

% Gerd Krahmann, LODYC Paris, Sep 1997
%
% Ref.: Y.Sezginer Unal & M.Ghil, Clim. Dyn., 1995, 11:255-278

% added 'rec'			G.Krahmann, LODYC/Paris Sep 1997, 0.1.0-->0.1.1


% check if NaN's exist in data
% if yes, replace them by white-noise random values
if isnan(sum(data(:)))
  dummy = sum(data);
  dummy = find(isnan(dummy));
  disp('WARNING: filling gaps with white noise')
  for n=1:length(dummy)
    bad = find(isnan(data(:,dummy(n))));
    va = nstd(data(:,dummy(n)));
    me = nmean(data(:,dummy(n)));
    if length(bad)>length(data)/10
      disp('WARNING: more than 10% gaps in series')
    end
    data(bad,dummy(n)) = randn(length(bad),1)*va + me;
  end
end


% check for length of lag in relation to length of time series
if llag>size(data,1)/3
  disp('WARNING: lag is longer than 1/3 of time series')
end


% remove mean of time series
data = data - repmat(mean(data),[size(data,1),1]);


% reorder data to lagged series
ind = zeros(llag,size(data,1)-llag+1,size(data,2));
dummy = [1:llag]';
ind = ind +...
	repmat(dummy,[1,size(data,1)-llag+1,size(data,2)]);
dummy = [0:size(data,1)-llag];
ind = ind +...
	repmat(dummy,[llag,1,size(data,2)]);
clear dummy
dummy(1,1,[1:size(data,2)]) = [0:size(data,2)-1]*size(data,1);
ind = ind +...
	repmat(dummy,[llag,size(data,1)-llag+1,1]);
ind = permute(ind,[2,1,3]);
ind = ind(:,:)';
ndata = data(ind);
sn = size(ndata);


% transpose the data matrix if sn(1)>sn(2)
% to get increased speed
if sn(1)>sn(2)
  ndata=ndata';
  transp=1;
else
  transp=0;
end


% calculate covariance matrix
% and divide by number of realizations
covmat=ndata*ndata'/size(ndata,2);


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


% normalize Eigenvectors to rms=1
rev=rms(ev);
ev=ev./(ones(size(ev,1),1)*rev);


% project to get amplitudes
amps=ndata'/ev';


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


% reshape results 
if intvers>4
  ev = reshape(ev,llag,size(ev,1)/llag,size(ev,2));
end

% reconstructed fields
if intvers>4 & nargout>3
  rec=zeros(size(data,1),size(data,2),size(ev,3));
  count=zeros(size(data,1),size(data,2),size(ev,3));
  for n=1:size(ev,3)
    for m=1:size(ev,1)
      for k=1:size(ev,3)
        rec(n+[0:size(ev,1)-1],:,k)=rec(n+[0:size(ev,1)-1],:,k)+...
		ev(:,:,n)*amps(n,k);
        count(n+[0:size(ev,1)-1],:,k)=count(n+[0:size(ev,1)-1],:,k)+1;
      end
    end
  end
  rec=rec./count;
end
