function [eof,lambda,tv,ev_percent,pc] = get_eof(varargin);
%[eof,lambda,tv,ev_percent,pc] = get_eof(x,do,n); 
%% do = 1d domain
%% x = data matrix (2d)
%% n = number of eofs to retain

if nargchk(3, 6, nargin);
  error('Check input arguments - must have between three and six')
end

X = varargin{1};
dox = varargin{2};

%determine if the data matrix is 3d or 2d

if nargin == 3 
    eofn = varargin{3};
    if eofn>6
        error('thats a lot of eofs, I doubt all of those are significant')
    end
    
    if ((length(size(X)) == 2) == 0)
        error('3 inputs requires a 2d data matrix, check input format')
    end
    
    
    Nsamples=length(X(1,:)); %data is organized with a time series at a single
    Nstructures=length(X(:,1)); %location as row; column represent all structures at one time
    display(['sampling dimension (most likely number of time steps) is ', num2str(Nsamples), ' and number of locations (structures) is ', num2str(Nstructures), ' if the reverse is intended, tranpose the input matrix'])
    [u,s,v]=svd(X, 0);  %zero inout is economy size
    
    
    %plot the first  eofn eo
    figure(1); clf;
    lt='EOF 1'
    for n=2:eofn
        lt=[lt; ['EOF ', num2str(n)]];
    end
    plotmulti(dox,v(:,1:eofn), lt ,1); 
    xlabel('structure domain'); ylabel('normalized structure value')
end
   
%    %plot the first 4 pcs
%    figure(2); clf; subplot(2,1,1);
%    plot(1:512,u(1,:).*s(1,1),'r',1:512, u(2,:).*s(2,2),'b', 1:512, u(3,:).*s(3,3),'g', 1:512,u(4,:).*s(4,4),'k' ); legend('PC 1','PC 2','PC 3', 'PC 4');
%    xlabel('time'); ylabel('principal component')
%    figure(2); hold on; subplot(2,1,2);
%    plot(ub(:,1:4));hold on; plot(1:1:512, -mean(x552b')/(32), 'b*')
%    legend('PC 1','PC 2','PC 3', 'PC 4', 'negative structure mean');
%    xlabel('time'); ylabel('principal component')
%    
%    figure(3); clf; subplot(2,1,1);
%    plot(1:1:32,(diag(s).^2)/N, 'b*');
%    xlabel('eof #'); ylabel('eigenvalue')
%    figure(3); hold on; subplot(2,1,2);
%    plot(1:1:32,(diag(sb).^2)/N, 'b*');
%    xlabel('eof #'); ylabel('eigenvalue')
%    
%   %Estimate the degrees of freedom in the data set
%   x=x552a; xb=x552b;
%   
% 
% %Begin by testing the auto.m function auto(x,m), where x is a data vector and m is 
% %the number of lags over which you want to compute the lag correlation. Here we 
% %ask to calculate the lag correlation of the first COLUMN of x out to 10 lags.
% auto(x(:,1),10)
% %now do the same for all the columns of x
%  for n = 1:32
%      cor(n,:)=auto(x(:,n),10);
%      corb(n,:)=auto(xb(:,n),10);
%  end
% 
% %Next plot the transpose  hof the cor matrix. This gives you a plot of the temporal 
% %autocorrelation functions of each of the 32 spatial grid points. Notice that 
% %they are all a little different because of sampling noise.
% figure(4); clf; subplot(2,1,1);
% plot(cor'); xlabel('time lag')
% ylabel('autocorrelation of each data structure')
% 
% figure(5); clf; subplot(2,1,1);
% plot(corb'); xlabel('time lag')
% ylabel('autocorrelation of each data structure')
% %Next we can set up a loop to calculate the average autocorrelation function for 
% %all 32 spatial grid points. First set the average array to zero.
% %Why would we do it this way, just take the mean of the columns of cor
% corm=mean(cor); cormb=mean(corb);
% figure(4); hold on; subplot(2,1,2)
% plot (corm); xlabel('time lag')
% ylabel('average autocorrelation over the 32 structures')
% figure(5); hold on; subplot(2,1,2)
% plot (cormb); xlabel('time lag')
% ylabel('average autocorrelation over the 32 structures')
% 
% %Now, estimate Nstar using the method of Brethereton
% Ns=N*((1-corm(1,2)^2)/(1+corm(1,2)^2));
% Nsb=N*((1-cormb(1,2)^2)/(1+cormb(1,2)^2));
% %calculate the error bars for each eigenvalue
% S=(diag(s).^2)/N; Sb=(diag(sb).^2)/N
% SE=S*(2/Ns)^.5; SEb=Sb*(2/Nsb)^.5;
% figure(3); subplot(2,1,1); hold on;
% errorbar(1:1:32, S, SE);
% figure(3); subplot(2,1,2); hold on;
% errorbar(1:1:32, Sb, SEb);



    
    
% %      [eof,d,tv,ev_percent,pc] = get_eof(x,n);
% %
% %    Calculates the EOF-s for a given field
% %    Actually:  x = x_cos
% %    eof - the eigenvectors (EOF)
% %    d - the eigenvalues 
% %    lambda -- array of the eigenvalues in descending order
% %    tv - the total variance
% %    ev_percent - the explained variance for the first n EOF-s in %
% %    pc - the principal components
% %    n - the number of modes

% % % Calculates the decomposition into singular values
% % % d - diagonal matrix with the eigenvalues on the diagonal
% % % eof - contains the eigenvectors (EOF) as its columns

%  % disp(sprintf('Calculating the first %g EOF-s and the eigenvalues',n))
%    [eof,d] = eig(cov);

%    [eof,d] = sorteig(eof,d);	% sorting in a descending order
   
% take only the number of eofs requested.
eof = eof(:,1:n);

%  disp('Calculating the total variance')
% calculate total variance
tv = trace(cov)/cols;
  
%  disp(sprintf('Calculating the explained variance for the %g EOF-s',n))
   z = sum(lambda);
   ev_percent = 100*lambda(1:n)/z;		% in percent ('100=%')
  
  disp('Calculating the principal components')
   pc = zeros(rows,n);
     for i = 1:n
       pc(:,i) = x*eof(:,i);
     end
    
%  disp('         Done!')









