function [Seof,EOF,pc,eig] = get_eof3sph(x,dox,doy,time,eofn);

% x=shiftdim(SLP(dow,:,:),1);
% dox=lon;
%doy=lat;
 %time=dec_year(dow);
 %eofn=4;


%takes grided 2d data and does the eof anaylsis ON THE SPHERE WE CALL EARTH
%ALSO COSINE WEIGHTS THE DATA MATRIX


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

%eliminate NaN from the anaylsis
Nstructures=nx*ny;

[XX, YY]=meshgrid(dox,doy); %standard grid of domain
XXc=g2c(XX); %column form
YYc=g2c(YY);


XXi=zeros(0); %define null matrix to compose good data points
YYi=zeros(0); % y domain null matrix
Xi=zeros(0); %2d data matrix null matrix
badstructures=zeros(0);

for n=1:Nstructures
    if length(find(isnan(X(n,:)) == 1)) > 0
        %XXi=XXi;  %these steps do nothing but may help see what the loop
        %YYi=YYi;  % does/ if there is a NaN at any time, eliminate the
        %Xi=Xi;   %data from the matrix- keep track with XXi and YYi
        
        display(['eliminating data at location [', num2str(XXc(n)),',',num2str(YYc(n)),'] which has ', num2str(length(isnan(X(n,:)))), ' NaN entries out of ', num2str(nt),' samples']);
        badstructures=[badstructures n]; %add that structure number to the list
    else
        XXi=[XXi; XXc(n)];
        YYi=[YYi; YYc(n)];
        Xi=[Xi; X(n,:).*cos(YYc(n)*pi/180)]; % Cosine Weighting
    end

end
display([ num2str(length(badstructures)), ' bad structures out of ', num2str(Nstructures), ' eliminated from the data matrix']);
    
     
        
%weight by the area



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


%create the EOFs and PCs, putting NaN back in at the appropriate spots

    
for n=1:eofn
    badnum=1;
    offset=0;
    seof=((Xi*v(:,n)).*(1./cos(pi*YYi/180)))/(nt^.5); %scaled eof - Hartman's notes are off on this one AD 1/28/07
   
    
    for nn=1:Nstructures
        
        if length(badstructures) > 0
    
        
        
            if badstructures(badnum) == nn
        
                Seof(nn,n)=NaN;
                EOF(nn,n)=NaN;
                offset=offset+1;
                badnum=badnum+1;
                if badnum > length(badstructures)
                    badnum=length(badnum); %once you've found the last bad structure, quit looking for bad structures
                end

            else
                Seof(nn,n)=seof(nn-offset);
                EOF(nn,n)=u(nn-offset,n);
            end %if bad structure number
   
    
        
        else  %if- not bad structures
           Seof(nn,n)=seof(nn-offset);
           EOF(nn,n)=u(nn-offset,n); 
        end %no bad structures case

    end %structures
    
end%eof number



    %plot the first eofn eof
  
    
    for n=1:eofn
        
        figure(n); clf;    
        sphcolor2(c2g(Seof(:,n),ny,nx),doy, dox, n); 
        title(['EOF ', num2str(n), '  (with magnitude of 1 standard deviation)'])
        eof(1:ny,1:nx,n)=c2g(EOF(:,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 ,eofn+1);    %%%%%   v(:,1:eofn))'
    xlabel('sampling domain'); ylabel('unnormalized pc value-note, the absolute magnitude is not that physically meaningful as it scales as 1 over the root of the structure dimension ')
    title('pcs of data')
    pc=pc(1:eofn,:);
    
   
    
                        

   
%Estimate the degrees of freedom in the data set
  


for n = 1:Nstructures-length(badstructures)
     cor(n,:)=auto(Xi(n,:),10);
     %warning off last
end
corm=mean(cor);


figure(eofn+2);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(eofn+2); 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')

%make the EOFs grided
SEOF=Seof;
EoF=EOF;
clear Seof EOF;

for n=1:eofn;
    Seof(n,1:ny,1:nx)=c2g(SEOF(:,n),ny,nx);
    EOF(n,1:ny,1:nx)=c2g(EoF(:,n),ny,nx);
end
    