function [Tr,fr,ir] = varim(A1,conve,nrs)
% VARIM       : Produces varimax function value and rotation matrix T.
%               Based on Nevels (1986). See also Pascal source.
%
% SYNTAX      : [Tr,fr,ir] = varim(A,conv,nrs)
%
% ARGUMENTS   : A    matrix to be rotated.
%               conv convergence criterion.
%               nrs  number of random starts.
%
% RETURN VALUE: T    rotation matrix.
%               f    varimax function value.
%               iter number of iterations.
%
% DATE WRITTEN: ? ???? ????. Henk Kiers.
% LAST UPDATE : 8 june 1994. Ren‚ van der Heijden.

A=A1;

[m,r] = size(A);
for nos=1:nrs+1,
   
  if nos == 1, T = eye(r);
  else         T = orth(rand(r,r) - .5); end;
      
  B=A*T;
  f=ssq((B.*B)-ones(m,1)*mean(B.*B));
  fold=f-2*conve*f;
  if f==0, fold=-conv; end;
  iter=0;
  while f-fold>f*conve
    fold=f;iter=iter+1;
    for i=1:r
      for j=i+1:r
        x=B(:,i);y=B(:,j);
        xx=T(:,i);yy=T(:,j);
        u=x.^2-y.^2;v=2*x.*y;
        u=u-ones(m,1)*mean(u);v=v-ones(m,1)*mean(v);
        a=2*sum(u.*v);b=sum(u.^2)-sum(v.^2);c=(a^2+b^2)^.5;
        e = (u'*u) - (v'*v);
        if a>=0; sign=1; end;
        if a<0; sign=-1; end;
        if c<.00000000001
          disp(' No rotation anymore');
           T = eye(3);
           return; 
          cos=1;sin=0;
          fr = 0; ir = 0; Tr = 0;
        end;
        if c>=.00000000001
          vvv=-sign*((b+c)/(2*c))^.5;
          sin=(.5-.5*vvv)^.5;cos=(.5+.5*vvv)^.5;
        end;
        v=cos*x-sin*y;w=cos*y+sin*x;
        vv=cos*xx-sin*yy;ww=cos*yy+sin*xx;
        if vvv>=0       % prevent permutation of columns
          B(:,i)=v;B(:,j)=w;T(:,i)=vv;T(:,j)=ww;
        end;
        if vvv<0
          B(:,j)=v;B(:,i)=w;T(:,j)=vv;T(:,i)=ww;
        end;
      end;
    end;
    f = ssq((B.*B)-ones(m,1)*mean(B.*B));
  end;
  if nos == 1
    fr = f; ir = iter; Tr = T;
  elseif f > fr
    fr = f; ir = iter; Tr = T;
  end;
end;
