function [V,D100,alpha,D] = EOF(X,k,l)
%用法：[VCN05,DCN05,alphaCN05] = EOF(CN05JJAlz',28)
%EOF 計算時空場X的EOF前10個模態(tài).
% V = eof(X) 計算X的前10個空間型，即V的10列
% [V,D100]=eof(X) 計算出前10個空間型以及每個空間型的方差解釋率
% [V,D100,alpha]=eof(X) 附加計算每個空間型對應的時間系數(shù)
% [V,D100,alpha,D]=eof(X) 附加計算每個空間型(特征向量)的特征值
% [......] = eof(X,k) 計算前k個模態(tài)
%
% X —— [m,n]矩陣，m為空間格點數(shù)帽借，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;
