function [B1,B2,T,f]=orthmax2(A1,A2,gam1,gam2,conv);
% function [B1,B2,T,f]=orthmax2(A1,A2,gam1,gam2,conv);
% produces simultaneous orthomax rotation of two matrices A and B 
% maximizes f = ssq (A1.^2) - 1/m1 * gam1 sum ((sum(A1.^2))^2 )
%             + ssq (A2.^2) - 1/m2 * gam2 sum ((sum(A2.^2))^2 )
% based on Jennrich 1970, explicit expression for coefficients from Clarkson & Jennrich
%
% input: A1,A2    (matrices to be rotated)
%	 gam1, gam2 (orthomax parameters)	
%	 conv (convergence criterion)
% output: T (rotation matrix)
%	  f (orthomax function value)


if nargin<5,conv=1e-6;opt=1;end;	% default convergence value

[m1,r]=size(A1);
[m2,r2]=size(A2);
if r2~=r, disp('column orders of A and B unequal');beep;return;end;

T=eye(r);
B1=A1;
B2=A2;

f=sum(B1(:).^4)-gam1/m1*sum((sum(B1.^2)).^2)+sum(B2(:).^4)-gam2/m2*sum((sum(B2.^2)).^2);
fold=f-2*conv*abs(f);
if f==0, fold=-conv; end;
iter=0;

while f-fold>abs(f)*conv 
  fold=f;iter=iter+1;
  for i=1:r
    for j=i+1:r
    
 	% Jennrich & Clarkson

        xx=T(:,i);yy=T(:,j);
        a=0;b=0;

        % for A1

        x=B1(:,i);y=B1(:,j);
        x2=x.^2;y2=y.^2;
        a=a+(-gam1/m1)*(.25*(sum(x2-y2))^2 - (sum(x.*y))^2) + .25*sum(x2.^2 + y2.^2 - 6*x2.*y2);
        b=b+(-gam1/m1)*sum(x.*y)*sum(x2-y2) + sum((x.^3).*y - x.*(y.^3));  

        % for A2

        x=B2(:,i);y=B2(:,j);
        x2=x.^2;y2=y.^2;
        a=a+(-gam2/m2)*(.25*(sum(x2-y2))^2 - (sum(x.*y))^2) + .25*sum(x2.^2 + y2.^2 - 6*x2.*y2);
        b=b+(-gam2/m2)*sum(x.*y)*sum(x2-y2) + sum((x.^3).*y - x.*(y.^3));  
        
        theta=0.25*atan2(b,a);
        cs=cos(theta);sn=sin(theta);
        x=B1(:,i);y=B1(:,j);
        B1(:,i)=cs*x+sn*y;      
        B1(:,j)=cs*y-sn*x;
        x=B2(:,i);y=B2(:,j);
        B2(:,i)=cs*x+sn*y;      
        B2(:,j)=cs*y-sn*x;

        T(:,i)=cs*xx+sn*yy;       
        T(:,j)=cs*yy-sn*xx; 
     end;
  end;  
  f=sum(B1(:).^4)-gam1/m1*sum((sum(B1.^2)).^2)+sum(B2(:).^4)-gam2/m2*sum((sum(B2.^2)).^2);
%  fprintf('f after %g iterations is %12.8f \n',iter,f);
end;