function [V,D100,alpha,D] = EOF(X,k,l)
%用法：[VCN05,DCN05,alphaCN05] = EOF(CN05JJAlz',28)
%EOF 计算时空场X的EOF前10个模态.
% V = eof(X) 计算X的前10个空间型，即V的10列
% [V,D100]=eof(X) 计算出前10个空间型以及每个空间型的方差解释率
% [V,D100,alpha]=eof(X) 附加计算每个空间型对应的时间系数
% [V,D100,alpha,D]=eof(X) 附加计算每个空间型(特征向量)的特征值
% [......] = eof(X,k) 计算前k个模态
%
% X —— [m,n]矩阵，m为空间格点数，n为时间长度
% 如果m>n，则利用时空转换

if nargin<1
error('MATLAB:eof:NotEnoughInputs', 'Not enough input arguments.');
end
if nargin<2
k = 10;
end
[m,n] = size(X);
if k>=m
error('MATLAB:eof:ErrorInputs', 'Error Inputs: "K".');
end

% X矩阵距平化或者标准化
if l == 1
X = X-repmat(mean(X,2),1,n); %mean(A,2)give the mean value along the row
elseif l == 2
X = (X-repmat(mean(X,2),1,n))./repmat(std(X')',1,n);
end

% 资料预处理，剔除0方差格点
%[X,inX] = EOF_PRE(X);
M = m;
m = m;
if m<n % 不转换时空
% 求解 XX'/n 的特征值和特征向量
[V,D] = eigs(XX'/(n-1),k);
D = diag(D); %get the diagonal value of the matrix
D = D(:);
alpha = V'X;
D100 = D100./trace(XX'/(m-1));
else % 时空转换
S = X'X/(m-1);
% 求解 S 的特征值和特征向量
[alpha,D] = eigs(S,k);
V = Xalpha;% 特征向量
D = diag(D);
D = D(:);
V = V./repmat(sqrt(D'(m-1)),m,1);
alpha= V'X;
D100 = D100./trace(X'*X/(m-1));
end
tmp = zeros(M,k);
tmp(:,:) = V;
V = tmp;