function rotation = ComputeRotation(output,k, coordinates, method, weights);
%Rotates the CA solution generated by CAnalysis.m
%
% INPUT
%       output      The 'output' structure generated by CAnalysis.m
%
%       k           Number of dimensions to retain
%
%       coordinates 1 -> Symmetrical coordinates
%                   2 -> Row principal coordinates, and column standard coordinates
%                   3 -> Column principal coordinates, and row standard coordinates
%
%       method      1 -> Orthogonal rotation of symmetrical coordinates
%                   2 -> Oblique rotation of symmetrical coordinates
%                   3 -> Orthogonal rotation of principal coordinates
%                   4 -> Orthogonal rotation of standard coordinates
%                   5 -> Oblique rotation of principal coordinates
%
%       weights     1 -> Row-wise normalization of coordinates
%                   2 -> To rescale coordinates using masses
%                   3 -> None. No weights scheme is used.
%
% OUTPUT
%
% 		rotation.coordinates    Coordinates model (1, 2, or 3, as expecified in the input)
% 		rotation.method         Rotation method computed (1, 2, 3, 4, or 5 as expecified in the input)
% 		rotation.weights        Weights scheme used in the rotation (1, 2, or 3 as expecified in the input)
% 		rotation.g              Function value of the rotation criterion
% 		rotation.g              Function value of the rotation criterion
% 		rotation.T              Tranfosmation matrix
% 		rotation.phi            Inter-dimension corralation matrix
% 		rotation.X_p            Rotated row coordiantes (Pattern) 
% 		rotation.Y_p            Rotated column coordiantes (Pattern)
% 		rotation.X_s            Rotated row coordiantes (Structure)
% 		rotation.Y_s            Rotated column coordiantes (Structure)
% 		rotation.X_ps           Rotated row coordiantes (Pseudo-Structure; only in oblique rotation of symmetrical coordinates)
% 		rotation.Y_ps           Rotated column coordiantes (Pseudo-Structure; only in oblique rotation of symmetrical coordinates)
%
%       rotation.X_sup_p        Rotated supplementary row coordiantes (Pattern)
%       rotation.X_sup_s        Rotated supplementary row coordiantes (Structure)
%       rotation.X_sup_ps       Rotated supplementary row coordiantes (Pseudo-Structure; only in oblique rotation of symmetrical coordinates)
%
%       rotation.Y_sup_p        Rotated supplementary column coordiantes (Pattern)
%       rotation.Y_sup_s        Rotated supplementary column coordiantes (Structure)
%       rotation.Y_sup_ps       Rotated supplementary column coordiantes (Pseudo-Structure; only in oblique rotation of symmetrical coordinates)
%
%	    rotation.bentler.r_b    Bentler's simplicity index before rotation for rows coordinates
%	    rotation.bentler.c_b    Bentler's simplicity index before rotation for columns coordinates
%	    rotation.bentler.r_a    Bentler's simplicity index after rotation for rows coordinates
%	    rotation.bentler.c_a    Bentler's simplicity index after rotation for columns coordinates

%       rotation.APCT           Absolute contributions: Contributions of points to axes (Principal coordinates)
%       rotation.ASCT           Absolute contributions: Contributions of points to axes (Standart coordinates)
%       rotation.ARSCT          Absolute contributions: Contributions of points to axes (Row symmetrical coordinates)
%       rotation.ACSCT          Absolute contributions: Contributions of points to axes (Column symmetrical coordinates)
%       rotation.RPCT           Relative contributions: Contributions of axes to points (Principal coordinates)
%       rotation.IN             Inertia of principal coordinates [raw percentage]
%       rotation.QLT            Measure of quality of display for principal coordinates in the reduced k-dimensional map. The qualities are equivalent to the communalities in PCA 
%       rotation.B              Correlation of principal coordinates with axes. Similar to component loadings in PCA
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rotation.coordinates = coordinates;
rotation.method      = method;
rotation.weights     = weights;

K=output.K;

%Checks for supplementary rows
SupRow=0; if isfield(output,'Fsup'),  SupRow=1; end;

%Checks for supplementary columns
SupCol=0; if isfield(output,'Gsup'),  SupCol=1; end;
    
switch method
% 1 -> Orthogonal rotation of symmetrical coordinates
case 1
    A1 = output.H(:,1:k);
    A2 = output.Z(:,1:k);
    Asr = output.H;
    Asc = output.Z;
    if SupRow, As_sr= output.Hsup(:,1:k);end;
    if SupCol, As_sc= output.Zsup(:,1:k);end;
    
    Dr = output.Dr;
    Dc = output.Dc;

    %weights scheme
    switch weights
    %row wise normalization
    case 1
        B1 = diag(diag(A1*A1'))^(-1/2) * A1;
        B2 = diag(diag(A2*A2'))^(-1/2) * A2;
    %weight for masses
    case 2
        B1 = Dr^(1/2) * A1;
        B2 = Dc^(1/2) * A2;
    %raw
    otherwise
        B1=A1;
        B2=A2;
    end;

    %rotation (varimax)
    [P1, P2, T1] = varimax2(B1,B2);

	rotation.T    = T1;

    %phi
    phi = T1'*T1;
	rotation.phi  = phi;
    
    %pattern 
	rotation.X_p     = A1*T1;
	rotation.Y_p     = A2*T1;
    if SupRow, rotation.X_sup_p = As_sr*T1;end;
    if SupCol, rotation.Y_sup_p = As_sc*T1;end;

	rotation.X_s     = rotation.X_p;
	rotation.Y_s     = rotation.Y_p;
    
    %bentler index before (b) rotation and after (a) rotation for rows (r) and columns (c) coordinates
	rotation.bentler.r_b = bentler(A1);
	rotation.bentler.c_b = bentler(A2);
	rotation.bentler.r_a = bentler(rotation.X_p);
	rotation.bentler.c_a = bentler(rotation.Y_p);

    % Diagnostics symmetrical case:
    Tt=[T1 zeros(k,K-k);zeros(K-k,k) eye(K-k)];
    L=diag(output.eig(1:k,1));
    varT=diag(diag(T1'*L*T1));
    L=diag(output.eig(:,1));
    varTt=diag(diag(Tt'*L*Tt));
    varTt4=varTt^0.25;
    PXTt=Asr*Tt*varTt4;
    PYTt=Asc*Tt*sqrt(sqrt(varTt));
    PXT=PXTt(:,1:k);
    PYT=PYTt(:,1:k);
    INr=sum((Dr*(PXTt.*PXTt))');
    INc=sum((Dc*(PYTt.*PYTt))');
    INr_p = INr / sum(output.eig(:,1));
    INc_p = INc / sum(output.eig(:,1));
    rotation.ARSCT=Dr*(PXT.*PXT)*inv(varT);
    rotation.ACSCT=Dc*(PYT.*PYT)*inv(varT);
    rotation.INSR = [INr' INr_p']; 
    rotation.INSC = [INc' INc_p'];
      
% 2 -> Oblique rotation of symmetrical coordinates
case 2
    
    A1 = output.H(:,1:k);
    A2 = output.Z(:,1:k);
    if SupRow, As_sr= output.Hsup(:,1:k);end;
    if SupCol, As_sc= output.Zsup(:,1:k);end;

    Dr = output.Dr;
    Dc = output.Dc;

    %weights scheme
    switch weights
    %row wise normalization
    case 1
        B1 = diag(diag(A1*A1'))^(-1/2) * A1;
        B2 = diag(diag(A2*A2'))^(-1/2) * A2;
    %weight for masses
    case 2
        B1 = Dr^(1/2) * A1;
        B2 = Dc^(1/2) * A2;
    %raw
    otherwise
        B1=A1;
        B2=A2;
    end;

    %rotation (quartimin)
    [T1, g] = Qmin(B1,B2);

    rotation.g    = g;
	rotation.T    = T1;

    %phi
    phi = T1'*T1;
	rotation.phi  = phi;
    
    %pattern and pseudo-structures
	rotation.X_p  = A1*inv(T1');
	rotation.Y_p  = A2*inv(T1');
	rotation.X_ps = A1*inv(T1')*phi^(1/2);
	rotation.Y_ps = A2*inv(T1')*phi^(1/2);
    if SupRow,
        rotation.X_sup_p  = As_sr*inv(T1');
        rotation.X_sup_ps = As_sr*inv(T1')*phi^(1/2);
    end;
    if SupCol,
        rotation.Y_sup_p  = As_sc*inv(T1');
        rotation.Y_sup_ps = As_sc*inv(T1')*phi^(1/2);
    end;
    
    %bentler index before (b) rotation and after (a) rotation for rows (r) and columns (c) coordinates
	rotation.bentler.r_b = bentler(A1);
	rotation.bentler.c_b = bentler(A2);
	rotation.bentler.r_a = bentler(rotation.X_p);
	rotation.bentler.c_a = bentler(rotation.Y_p);

% 3 -> Orthogonal rotation of principal coordinates
case 3
    %coordinates
    switch coordinates
    %Row principal coordinates, and column standard coordinates
    case 2
        A1 = output.F(:,1:k);
        A2 = output.Y(:,1:k);
        D   = output.Dr;
        Ap  = output.F;
        Dp  = output.Dr;
        Ast  = output.Y;
        Dst  = output.Dc;
        if SupRow, As_sr= output.Fsup(:,1:k);end;
        if SupCol, As_sc= output.Ysup(:,1:k);end;

    %Column principal coordinates, and row standard coordinates
    case 3
        A1 = output.G(:,1:k);
        A2 = output.X(:,1:k);
        D   = output.Dc;
        Ap  = output.G;
        Dp  = output.Dc;
        Ast = output.X;
        Dst = output.Dr;
        if SupRow, As_sr= output.Xsup(:,1:k);end;
        if SupCol, As_sc= output.Gsup(:,1:k);end;
    %error
    otherwise
        error('Coordinates are not coherent with the rotation method');
        return;
    end;
        
    %weights scheme
    switch weights
    %row wise normalization
    case 1
        B1 = diag(diag(A1*A1'))^(-1/2) * A1;
    %weight for masses
    case 2
        B1 = D^(1/2) * A1;
    %raw
    otherwise
        B1=A1;
    end;

    %rotation (varimax)
    [T1,g] = varim(B1,0.00001,10);

    rotation.g    = g;
	rotation.T    = T1;
    
    %phi
    phi = T1'*T1;
	rotation.phi  = phi;
    
    %pattern
    switch coordinates
    %Row principal coordinates, and column standard coordinates
    case 2
    	rotation.X_p  = A1*T1;
    	rotation.Y_p  = A2*T1;
    %Column principal coordinates, and row standard coordinates
    case 3
    	rotation.X_p  = A2*T1;
    	rotation.Y_p  = A1*T1;
    end;
    if SupRow, rotation.X_sup_p  = As_sr*T1; end;
    if SupCol, rotation.Y_sup_p  = As_sc*T1; end;
	rotation.X_s     = rotation.X_p;
	rotation.Y_s     = rotation.Y_p;
        
    %bentler index before (b) rotation and after (a) rotation for rows (r) and columns (c) coordinates
	rotation.bentler.r_b = bentler(A1);
	rotation.bentler.c_b = bentler(A2);
	rotation.bentler.r_a = bentler(rotation.X_p);
	rotation.bentler.c_a = bentler(rotation.Y_p);    

    %diagnostics
    Tt=[T1 zeros(k,K-k);zeros(K-k,k) eye(K-k)];
    L=diag(output.eig(1:k,1));
    varT=diag(diag(T1'*L*T1));
    %inertiaT=diag((1/output.ti)*varT);
    PT=Ap*Tt;
    ST=Ast*Tt;
	IN=sum((Dp*(PT.*PT))');
    IN_p = IN / sum(output.eig(:,1));
	rotation.APCT=Dp*((PT(:,1:size(varT,1)).*(PT(:,1:size(varT,1)))))*inv(varT);
    rotation.ASCT=Dst*((ST(:,1:size(varT,1))).*(ST(:,1:size(varT,1))));
	rotation.RPCT=inv(diag(IN))*Dp*(PT(:,1:size(varT,1)).*PT(:,1:size(varT,1)));
    rotation.IN = [IN' IN_p'];                  %Inertia [raw percentage]
    QCOR = rotation.RPCT;
    rotation.QLT = QCOR*ones(k,1);
    rotation.B = sign(PT(:,1:k)).*sqrt(QCOR);
    
% 4 -> Orthogonal rotation of standard coordinates
case 4
    %coordinates
    switch coordinates
    %Row principal coordinates, and column standard coordinates
    case 2
        A1 = output.F(:,1:k);
        A2 = output.Y(:,1:k);
        D  = output.Dc;
        Ap  = output.F;
        Dp  = output.Dr;
        Ast = output.Y;
        Dst = output.Dc;
        
        Dd = output.Dr;
        
        if SupRow, As_sr= output.Fsup(:,1:k);end;
        if SupCol, As_sc= output.Ysup(:,1:k);end;
    %Column principal coordinates, and row standard coordinates
    case 3
        A1 = output.G(:,1:k);
        A2 = output.X(:,1:k);
        D  = output.Dr;
        Ap = output.G;
        Dp = output.Dc;        
        Ast= output.X;
        Dst= output.Dr;
        
        Dd = output.Dc;
        
        if SupRow, As_sr= output.Xsup(:,1:k);end;
        if SupCol, As_sc= output.Gsup(:,1:k);end;
    %error
    otherwise
        error('Coordinates are not coherent with the rotation method');
        return;
    end;
    
    %weights scheme
    switch weights
    %row wise normalization
    case 1
        B2 = diag(diag(A2*A2'))^(-1/2) * A2;
    %weight for masses
    case 2
        B2 = D^(1/2) * A2;
    %raw
    otherwise
        B2=A2;
    end;

    %rotation (varimax)
    [T1,g] = varim(B2,0.00001,10);

    rotation.g    = g;
	rotation.T    = T1;

    %phi
    phi = T1'*T1;
	rotation.phi  = phi;

    %pattern
    switch coordinates
    %Row principal coordinates, and column standard coordinates
    case 2
    	rotation.X_p  = A1*T1;
    	rotation.Y_p  = A2*T1;
    %Column principal coordinates, and row standard coordinates
    case 3
    	rotation.X_p  = A2*T1;
    	rotation.Y_p  = A1*T1;
    end;
	rotation.X_s     = rotation.X_p;
	rotation.Y_s     = rotation.Y_p;

    if SupRow, rotation.X_sup_p  = As_sr*T1; end;
    if SupCol, rotation.Y_sup_p  = As_sc*T1; end;
    
    %bentler index before (b) rotation and after (a) rotation for rows (r) and columns (c) coordinates
	rotation.bentler.r_b = bentler(A1);
	rotation.bentler.c_b = bentler(A2);
	rotation.bentler.r_a = bentler(rotation.X_p);
	rotation.bentler.c_a = bentler(rotation.Y_p);    

    %diagnostics
    Tt=[T1 zeros(k,K-k);zeros(K-k,k) eye(K-k)];
    L=diag(output.eig(1:k,1));
    varT=diag(diag(T1'*L*T1));
    %inertiaT=diag((1/output.ti)*varT);
    PT=Ap*Tt;
    ST=Ast*Tt;
	rotation.APCT=Dp*((PT(:,1:size(varT,1)).*(PT(:,1:size(varT,1)))))*inv(varT);
    rotation.ASCT=Dst*((ST(:,1:size(varT,1))).*(ST(:,1:size(varT,1))));
	IN=sum((Dp*(PT.*PT))');
    IN_p = IN / sum(output.eig(:,1));
    rotation.IN = [IN' IN_p'];                  %Inertia [raw percentage]

% 5 -> Oblique rotation of principal coordinates
case 5
    
    %coordinates
    switch coordinates
    %Row principal coordinates, and column standard coordinates
    case 2
        A1 = output.F(:,1:k);
        A2 = output.Y(:,1:k);
        D  = output.Dr;
        if SupRow, As_sr= output.Fsup(:,1:k);end;        
        if SupCol, As_sc= output.Ysup(:,1:k);end;        
    %Column principal coordinates, and row standard coordinates
    case 3
        A1 = output.G(:,1:k);
        A2 = output.X(:,1:k);
        D  = output.Dc;
        if SupRow, As_sr= output.Xsup(:,1:k);end;
        if SupCol, As_sc= output.Gsup(:,1:k);end;
    %error
    otherwise
        error('Coordinates are not coherent with the rotation method');
        return;
    end;
    
    %weights scheme
    switch weights
    %row wise normalization
    case 1
        B1 = diag(diag(A1*A1'))^(-1/2) * A1;
    %weight for masses
    case 2
        B1 = D^(1/2) * A1;
    %raw
    otherwise
        B1=A1;
    end;

    %rotation (quartimin)
    [Pr, phi, g ] = oblimin (B1, 0,0.00001,10);
    U1 = (B1'*B1)^(-1)*B1'*Pr;
    T1=(U1')^(-1);

    rotation.g    = g;
	rotation.T    = T1;

    %phi
	rotation.phi  = phi;
    
    %pattern and structures
    switch coordinates
    %Row principal coordinates, and column standard coordinates
    case 2
		rotation.X_p  = A1*inv(T1');
		rotation.Y_p  = A2*inv(T1');
		rotation.Y_s  = A2*inv(T1')*phi;
        if SupRow,
			rotation.X_sup_p  = As_sr*inv(T1');
        end;
        if SupCol,
			rotation.Y_sup_p  = As_sc*inv(T1');
			rotation.Y_sup_s  = As_sc*inv(T1')*phi;
        end;
    %Column principal coordinates, and row standard coordinates
    case 3
		rotation.Y_p  = A1*inv(T1');
		rotation.X_p  = A2*inv(T1');
		rotation.X_s  = A2*inv(T1')*phi;
        if SupRow,
			rotation.X_sup_p  = As_sr*inv(T1');
			rotation.X_sup_s  = As_sr*inv(T1')*phi;
        end;
        if SupCol,
			rotation.Y_sup_p  = As_sc*inv(T1');
        end;
    end;
    
    %bentler index before (b) rotation and after (a) rotation for rows (c) and columns (c) coordinates
	rotation.bentler.r_b = bentler(A1);
	rotation.bentler.c_b = bentler(A2);
	rotation.bentler.r_a = bentler(rotation.X_p);
	rotation.bentler.c_a = bentler(rotation.Y_p);
    
end;

%-----------------------------------------------------------------------------------------------
function [APCT, ASCT, RPCT, RSCT]= RotatedDiagnostics(P,S,T,varT,l,Dp,Ds)

		ST=S*T;
		PT=P*T;
		varTt=diag(diag(T'*(l*l)*T));
		
		totp=sum((Dp*(PT.*PT))');
		tots=sum(((Ds*(ST.*ST))*varTt)');
		
		APCT=Dp*((PT(:,1:size(varT,1)).*(PT(:,1:size(varT,1)))))*inv(varT);
		ASCT=Ds*((ST(:,1:size(varT,1))).*(ST(:,1:size(varT,1))));
		
		RPCT=inv(diag(totp))*Dp*(PT(:,1:size(varT,1)).*PT(:,1:size(varT,1)));
		RSCT=inv(diag(tots))*Ds*(ST(:,1:size(varT,1)).*ST(:,1:size(varT,1)))*varT;
return;