function [EOFs,expl,PCs]=eof(xdata);
%A program to do EOF analysis, and scale the EOFs.
% This is an example of scaling EOFs
% usage: function [EOFs,expl,PCs]=eof(xdata)
% xdata[m,n], m: spatial dimension, n:temporal dimension
% output: EOFs, expl, PCs
[M,N]=size(xdata);
xm=xdata;
% calculate covariance matrix
C=xm*xm'/N;
%C
% Do eigenanalysis
[EOFs,D]=eig(C);
PCs=EOFs'*xm;
%calculate autocorrelation
for j=1:M
alpha(j,:)=auto(xm(j,:),1);
end
% the temporal autocorrelation is different for each spatial grid point,
% a problem. Let's kluge this by using the mean of the mean and the maximum
% autocorrelation. We'll take the absolute value before averaging, so that
% random negative correlations can't cancel out positive ones. This is
% conservative.
ralph=(mean(abs(alpha(:,2)))+max(alpha(:,2)))/2.;
tau=-1.0/log(ralph);
Nstar=N/(2*tau);
factor=sqrt(2./Nstar);
% Plot Eigenvalues
L=diag(D,0);
LE=L*factor;
Lexp=L(:);
wkp=100*Lexp/sum(Lexp);
expla=flipud(wkp);
[ntmp1,ntmp2]=size(expla);
stmp=0.0;
for ikp=1:ntmp1
stmp=stmp+expla(ikp);
wkexpl(ikp)=stmp;
end
expl=[expla(:) wkexpl(:)]; % [m,2]
KKa= flipud(L);KKb=flipud(LE);
errorbar(1:10,KKa(1:10),KKb(1:10));
title(['Eigenvalue Spectrum, first 10']);
xlabel('Number');
EOFs=fliplr(EOFs);% [m,m] now the most important five are EOFs[:,1:5]
% Plot eigenvectors
%plot(E(:,M-2:M));
%MP=M+1;
%axis([0 MP -0.60 0.60]);
%set(gca,'DataAspectRatio',[1 22 1])
%title(['First Three Eigenvectors: Norm= ' ans]);
%xlabel('Structure dimension');
PCs=flipud(PCs);% ?[m,n] now the most important five are PCs[1:5,:]
PCs=PCs';%[n,m] now the most important five are PCs[:,1:5]
% Project eigenvectors onto original data to get PCs
%Z=E'*x;
%plot(Z(M-2:M,:)');
%title(['First Three Principle Components: Norm= ' ans]);
%xlabel('Sampling dimension');
% clear stmp ikp ntmp1 ntmp2 wkexpl expla wkp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%程序里的auto()没有呀
%把auto改成autocorr就可以了
