% This procedure follows algorithm as spelled out in
% Harman (1960) in Chapter 14, section 4.  To run the
% program - the loadings are put in an array called
% lding.  Type return to continue processing.
% The notation follows Harman.  The routine vfunct.m is
% called to compute the variance of the loadings
% squared.
disp('Set a variable called lding equal to your loadings')
disp('Type "return" to resume Varimax Rotation')
keyboard
b=lding;
[n,nf]=size(lding)
disp('PAUSE, hit any key to continue');pause
hjsq=diag(lding*lding')   % communalities
hj=sqrt(hjsq);
disp('PAUSE, hit any key to continue');pause
vfunct                    % function to compute
                          % variances of loadings^2
V0=Vtemp
disp('PAUSE, hit any key to continue');pause
for it=1:10; % Never seems to need very many 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)>.00001;
          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;
vfunct;
V=Vtemp;
if abs(V-V0)<.0001;break;else V0=V;end;
end;
lding
V
disp('PAUSE, hit any key to continue');pause