eofrot.m

function [coef,Arot,varrotated]=eofrot(z,index)
  
% Algorithm for eof rotation with the varimax method
% Difiere de la eofrot de Antonio, pues le agregue en el "sort" los RPC
% (coef)
%
%Inputs:
%       z            Data Matrix
%       index        EOF to be rotated
%
%Outputs:
%       coef            Cofficient for the rotated EOFs
%       Arot            Rotated EOF in ascending order
%       varrotated      Variance explained by rotated EOFs
  
  
%  Time and space points
[npoints,ntime]=size(z);
  
% Compute unrotated EOF for variance calculation
[uneof,ss,vneof]=svd(z,0); diag(ss);
totvar = sum(diag(ss.^2));
  
lding = uneof(:,index);
sl    = diag(ss(index,index).^2);
% variance explained by the unrotated modes
varexpl = sum(sl)/totvar;  
  
[n,nf]=size(lding);
b=lding;
hl=lding*(diag(sl));
  
hjsq=diag(hl*hl');           
hj=sqrt(hjsq);                                  %  Normalize by the communalities
bh=lding./(hj*ones(1,nf));  
Vtemp=n*sum(sum(bh.^4))-sum(sum(bh.^2).^2);     % VARIMAX functional to be minimized
V0=Vtemp;
for it=1:10;                    % Number of iterations
for i=1:nf-1;                   %     Program cycles through 2 factors
  jl=i+1;                       %     at a time
  for j=jl:nf;
      xj=lding(:,i)./hj;        % notation here closely
      yj=lding(:,j)./hj;        % follows harman
      uj=xj.*xj-yj.*yj;
      vj=2*xj.*yj;
      A=sum(uj);
      B=sum(vj);
      C=uj'*uj-vj'*vj;
      D=2*uj'*vj;
      num=D-2*A*B/n;
      den=C-(A^2-B^2)/n;
      tan4p=num/den;
      phi=atan2(num,den)/4;
angle=phi*180/pi;
      [i j it angle];
      if abs(phi)> eps;
          Xj=cos(phi)*xj+sin(phi)*yj;
          Yj=-sin(phi)*xj+cos(phi)*yj;
          bj1=Xj.*hj;
          bj2=Yj.*hj;
          b(:,i)=bj1;
          b(:,j)=bj2;
          lding(:,i)=b(:,i);
          lding(:,j)=b(:,j);
      end
  end
end;
lding=b;
  
bh=lding./(hj*ones(1,nf));
Vtemp=n*sum(sum(bh.^4))-sum(sum(bh.^2).^2);
V=Vtemp;
fprintf('  Converging  V  V0 \n');
fprintf(' %f\n',[V V0])
  
if abs(V-V0)<.0001;break;else V0=V;end;
end;
  
% Reflect vectors with negative sums
  
for i = 1:nf
    if sum(lding(:,i)) < 0
        lding(:,i) = -lding(:,i);
    end
end
  
Arot=lding ;                        % rotated eof
coef=z'*Arot(:,1:nf);               % time series for rotated eof
whos Arot coef
  
for i=1:nf
varex(i) = sum(var(coef(:,i)*Arot(:,i)')*(ntime-1));
end
  
varexpl
varexplrot = sum(varex)/totvar
zvar=sum(var(z')*(ntime-1))
totvar
  
% Sort in decreasing order of variance
whos varex
varex   
[varex,I]=sort(varex);
Arot=Arot(:,I);
Arot = fliplr(Arot);
coef=coef(:,I);
coef=fliplr(coef);
varex = flipud(varex');
varunrotated = sl/totvar
varrotated   = varex/totvar
return
PS：我从别的地方copy的{B.Aires, 20-24/02/06 - Centro de Investigaciones del Mar y la Atmosfera, Department of Atmospheric and Oceanic Sciences (UBA)
E.Scoccimarro, A.F.Carril }