function [eof,pc,eig,s] = get_eof2(x,dox,doy,time,eofn);

%takes grided 2d data and does the eof anaylsis


if length(size(x)) < 3 | length(size(x)) > 3
    
    error('data matrix must be 3d for this format, try get_eof for 2d matrix')
    
end

[ny, nx, nt]= size(x);

if ((ny == length(doy)) == 0)
    
    error('data matrix and y domain dimensional mismatch - y position should change moving down a column (y is 1st dimension in data)')
    
end

if ((nx == length(dox)) == 0)
    
   error('data matrix and x domain dimensional mismatch - x position should change moving across a row (x is 2nd dimension in data)')
   
end

if ((nt == length(time)) == 0)
    
   error('data matrix and time domain dimensional mismatch - time should be third dimension in array, try shiftdim if you dont dig it')
   
end

if eofn>8
    error('thats a lot of eofs, I doubt all of those are significant')
end
    
display([' number of time steps is ', num2str(length(time)), ' and grid of data is ', num2str(length(dox)), ' in x and ', num2str(length(doy)),' in y'])

%compose the data matrix from the 3d matrix

for t=1:nt
    X(1:nx*ny, t)=g2c(x(:,:,t));
end

[u,s,v]=svd(X, 0);  %zero inout is economy size
    

%create a mesh of the spatial domain

[XX, YY]=meshgrid(dox,doy);

    %plot the first eofn eof
    figure(1); clf;
    
    for n=1:eofn
        hold on;
        subplot(eofn,1,n);
        pcolor(XX, YY, c2g(u(:,n),ny,nx)); colorbar;
        title(['EOF ', num2str(n), '  (normalized)'])
        
        eof(1:ny,1:nx,n)=c2g(u(:,n),ny,nx);
        
    end
    
    %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))'
    xlabel('sampling domain'); ylabel('unnormalized pc value-note, the absolute magnitude is not that physically meaningful as it scales as 1 of the root of the structure dimension ')
    title('pcs of data')
    pc=pc(1:eofn,:);
                        

   
%Estimate the degrees of freedom in the data set
  
Nstructures=nx*ny;

for n = 1:Nstructures
     cor(n,:)=auto(X(n,:),10);
end
corm=mean(cor);


figure(3);clf;
subplot(2,1,1); plot(0:10,corm)
xlabel('time lag (number of points)')
ylabel('mean autocorrelation')
title('mean autocorrelation (averaged over all structures')
    




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

%calculate the error bars for each eigenvalue
S=(diag(s).^2)/nt;  
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), '--');
xlabel('EOF PC pair number'); ylabel('percent variance explained')