function [eof,pc,eig,s] = get_eof(varargin);





%[eof,lambda,tv,ev_percent,pc] = get_eof(x,do,time,n); %2d data structure
%[eof,lambda,tv,ev_percent,pc] = get_eof(x,dox,doy,time,n); %3d data
%[eof,lambda,tv,ev_percent,pc] = get_eof(x,dox,doy,time,n, 'sphere'); %3d data sphere


%% do = 1d domain
%% x = data matrix (2d)
%% n = number of eofs to retain

if nargchk(4, 4, nargin);
  error('Check input arguments - must have between four and four')
end

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

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

if nargin == 4 
    time = varargin{3};
    eofn = varargin{4};
    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, checwhoik 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
    lt='EOF 1';
    for n=2:eofn
        lt=[lt; ['EOF ', num2str(n)]];
    end
    
    for n=1:eofn
        Seof(:,n)=X*v(:,n)/Nsamples;
    end
    
    
    plotmulti(dox,Seof(:,1:eofn)', lt ,1); 
    xlabel2('structure domain'); ylabel2('structure value associated with 1 standard deviation change in pc')
    title2('EOFS of data')
    set2;
    %plot the pcs
    
    pc=s*v';
    lt='PC 1';
    for n=2:eofn
        lt=[lt; ['PC ', num2str(n)]];
    end
    plotmulti(time,pc(1:eofn,:), lt ,2);    %%%%%   v(:,1:eofn))'
    xlabel2('sampling domain'); ylabel2('unnormalized pc value-note')
    title2('pcs of data')
    set2;
    pc=pc(1:eofn,:);
    eof=u(:,1:eofn);
                        

end
   
%Estimate the degrees of freedom in the data set
  
for n = 1:Nstructures
     cor(n,:)=auto(X(n,:),10);
end
corm=mean(cor);

%plot the autocorrelation for each structure
if nargin == 4
    [LAG, XX]=meshgrid(0:10,dox);
    figure(3); clf;
    subplot(2,1,1); pcolor(LAG,XX, cor); colorbar;
    xlabel2('time lag (number of points)')
    ylabel2('structure domain value')
    title2('autocorrelation for each structure')
    set2
else
    figure(3);clf;
    subplot(2,1,1); plot(0:10,corm)
    xlabel2('time lag (number of points)')
    ylabel2('mean autocorrelation')
    title2('mean autocorrelation')
    set2
end




% %Now, estimate Nstar using the method of Brethereton
Ns=Nsamples*((1-corm(1,2)^2)/(1+corm(1,2)^2));

%calculate the error bars for each eigenvalue
S=(diag(s).^2)/Nsamples;  
SE=S*(2/Ns)^.5; eig=SE;
PS=S/sum(S); %percent of variance explained
PSE=PS*(2/Ns)^.5;  %error in percent variance explained

figure(3); subplot(2,1,2); hold on;
plot(1:1:eofn, PS(1:eofn), 'r*')
errorbar(1:1:eofn, PS(1:eofn), PSE(1:eofn), '--');
xlabel2('EOF PC pair number'); ylabel2('percent variance explained')
set2;


    
    










